Contents

A&A 443, 1033-1046 (2005)
DOI: 10.1051/0004-6361:20042272

Dissipation of Alfvén waves in complex 3D coronal force-free structures[*]

F. Malara - M. F. De Franceschis - P. Veltri

Dipartimento di Fisica, Università della Calabria, via P. Bucci, 87036 Rende (CS), Italy

Received 28 October 2004 / Accepted 19 May 2005

Abstract
The propagation and dissipation of Alfvén waves in a three-dimensional inhomogeneous force-free equilibrium structure is considered. Assuming small wavelengths the wave propagation is described within the framework of a WKB theory. The equilibrium structure is a complex magnetic field, resulting from a superposition of several harmonics, which represents a model of the coronal magnetic field above a quiet-Sun region. The properties of wave dissipation are related to the topological complexity of the background magnetic field, and the typical scale law of fast dissipation is recovered. A partial reflection of waves is included in the model by calculating reflection and transmission coefficients of both Alfvén and magnetosonic waves. An energy balance is derived between the wave energy dissipated inside the coronal structure and the energy lost through the lower boundary. We find that even for large values of the Reynolds number, a non-negligible fraction of the Alfvén wave energy is dissipated inside the corona.

Key words: magnetohydrodynamics (MHD) - waves - Sun: magnetic fields

1 Introduction

Among the mechanisms invoked to explain the nonradiative heating of the solar corona, dissipation of Alfvén waves represents one of the most studied from the theoretical point of view. The ubiquitous presence of a strong magnetic field, which represents the backbone of coronal structures, allows Alfvénic fluctuations to propagate from distant part, carrying both mechanical and electromagnetic energy. Alfvénic fluctuations have never been directly observed in the corona. However, they represent the main component of solar wind fluctuation spectrum (Belcher & Davis 1971) and it is likely they are also present in the corona, from where the solar wind is emanated. Due to the very low dissipative coefficients (resistivity and relevant viscosity) of the coronal plasma, an important point of wave-based theories is how to dissipate fluctuations before they leave the corona: dissipation can take place only if the wave energy is efficiently transferred to small scales. In highly structured magnetic fields as in the corona, interaction between Alfvén waves and inhomogeneities of the background structure naturally leads to generation of small scales. In this framework, we are not interested in specific dissipation mechanisms active at small scales. Rather, we focus on the dynamical processes that generate small-scale fluctuations.

In a 2D inhomogeneous background, where the Alfvén velocity $c_{\rm A}$ varies in a direction perpendicular to the magnetic field, two mechanisms have been studied in detail: (1) phase-mixing (Heyvaerts & Priest 1983; Nocera et al. 1986; Parker 1991; Nakariakov et al. 1997; Moortel et al. 2000; Botha et al. 2000; Tsiklauri et al. 2001, 2002b, 2003; Tsiklauri & Nakariakov 2002a), in which differences in phase velocity at different locations progressively bend wavefronts; and (2) resonant absorption which concentrates the wave energy in a narrow layer where the wave frequency locally matches a characteristic frequency (Alfvén or cusp). These processes have been studied both investigating normal modes of the inhomogeneous structure (Kappraff & Tataronis 1977; Mok & Einaudi 1985; Steinolfson 1985; Davila 1987, 1991; Califano et al. 1990, 1992; Hollweg 1987a,b), and by considering the time evolution of an initial disturbance (Lee & Roberts 1986; Malara et al. 1992, 1996). Effects of density stratification and magnetic line divergence (Ruderman et al. 1998), as well as linear (Tsiklauri & Nakariakov 2002a; Tsiklauri et al. 2003) and nonlinear (Nakariakov et al. 1997, 1998) coupling with compressive modes have also been considered. Recently, Tsiklauri et al. (2005) have studied the phase-mixing of Alfvén waves in a kinetic description, showing that the scale law of dissipation derived within the MHD is also valid in that context.

Similon & Sudan (1989) considered Alfvén wave propagation in three-dimensional magnetic structures. In this case, there can be regions where magnetic lines are chaotic: in these regions initially nearby magnetic lines exponentially move apart with increasing distance along the lines (Zimbardo et al. 1984, 1995). Then, Alfvénic disturbances that locally follow magnetic lines are stretched exponentially in time, leading to a fast formation of small scales. As a consequence, the typical dissipation time $t_{\rm d}$ of fluctuations follows the scaling law:

 \begin{displaymath}
t_{\rm d} \propto \log S
\end{displaymath} (1)

where S is the relevant viscous and/or resistive Reynolds number. Thus, the dissipation time remains relatively small even for large values of the Reynolds number, as in the coronal plasma. The mechanism proposed by Similon & Sudan has been investigated by Petkaki et al. (1998) and by Malara et al. (2000), who built a model based on a WKB expansion of incompressible magnetohydrodinamic (MHD) equations. In this model propagation and dissipation of Alfvénic wavepackets in 3D magnetic fields is described. From this model they found that the scaling law (1) corresponding to fast dissipation of wavepackets is recovered in regions of chaotic magnetic lines. This effect can co-exist with a slower dissipation related to phase-mixing, in regions where the topology of magnetic lines is quasi-2D. More recently, Malara et al. (2003) extended the above model to the case of a cold plasma, which is more suitable than the incompressible approximation to describe the low-$\beta$ coronal plasma.

The equations derived by Malara et al. (2003) describing the evolution of Alfvénic wavepackets are quite general and they can be used independently of the specific form of the background structure. However, the background magnetic fields considered by Petkaki et al. (1998) and by Malara et al. (2000, 2003) do not represent a realistic model of coronal magnetic fields. The purpose of the present paper is to apply the above theory to a more realistic description of coronal magnetic fields. In the present paper we will consider a force-free equilibrium magnetic field which can be used as a model of the coronal magnetic field above a quiet-Sun region. Propagation and dissipation of Alfvénic wavepackets within such a structure will be described by means of the equations derived by Malara et al. (2003). This equilibrium magnetic field was first proposed by Nakagawa & Raadu (1972) within the problem of deriving the magnetic field in the corona using measurements in the photosphere. The mathematical expression is quite general and it can describe a complex magnetic structure with several regions of the two polarities. The complex 3D topology of magnetic lines results in a fast increase of the wavevector of Alfvénic packets and fast dissipation. This magnetic field represents a closed structure because magnetic lines mostly start and end at the coronal base. Alfvénic packets which propagate along fieldlines are partially trapped inside this closed structure and, in order to describe their dynamics, it is necessary to represent the partial reflection process of packets at the coronal base. This will be done with a simple model which takes into account the geometry of the equilibrium magnetic field. Then, the wave energy deposition inside the corona will result as a balance between dissipation and energy losses through the partially reflecting boundary at the coronal base.

2 Alfvén waves in an inhomogeneous magnetic field

The coronal plasma is dominated by the magnetic force and it is then characterized by a low value of the plasma $\beta$ ($\beta$ being the gas pressure to magnetic pressure ratio). In the corona, typically $\beta \sim 10^{-2}$. In this case, a usual representation of large scale phenomena is that given by the cold plasma MHD equations, in which pressure and gravity terms are neglected with respect to the magnetic force in the momentum equation. These equations are written in the following dimensionless form:

 \begin{displaymath}
{\partial \rho \over \partial t}+v_j{\partial \rho \over \partial x_j}+
\rho{\partial v_j \over \partial x_j} = 0
\end{displaymath} (2)


 \begin{displaymath}
\rho\left[ {\partial v_i \over \partial t}+
v_j{\partial v_i...
...\partial \over \partial x_i}
{\partial v_j\over \partial x_j}
\end{displaymath} (3)


 \begin{displaymath}
{\partial B_i \over \partial t}+v_j{\partial B_i \over \part...
...lambda_N}{\partial^2 B_i \over \partial x_j\partial x_j} \cdot
\end{displaymath} (4)

In these equations, lengths have been normalized to a characteristic length L; the magnetic field is normalized to a characteristic magnetic field  ${\hat B}_0$, the density to a characteristic density ${\hat \rho}_0$, the velocity to the corresponding Alfvén velocity ${\hat c}_{\rm A0}={\hat B}_0/(4\pi {\hat \rho}_0)^{1/2}$, and the time to the Alfvén time $t_{\rm A}=L/{\hat c}_{\rm A0}$. The dimensionless dissipative coefficients are $\eta_N =\eta /{\hat \rho}_0 L {\hat c}_{\rm A0}$, $\zeta_N =\zeta /{\hat \rho}_0 L {\hat c}_{\rm A0}$, and $\lambda_N =\lambda / L {\hat c}_{\rm A0}$, where $\eta$ and $\zeta$ are the viscosity coefficients and $\lambda$ is the magnetic diffusivity, which are assumed to be constant. Summation over dummy indices is hereafter assumed.

We consider the propagation and the dissipation of Alfvén waves in a 3D inhomogeneous magnetic structure, representing a model of the magnetic field in the solar corona, above a quiet Sun region. A full description of wave propagation in a complex 3D magnetic field represents a formidable task, since different wave modes (Alfvén and magnetosonic) are coupled among them and with the inhomogeneity of the background structure. A possible approach to such a problem is represented by direct numerical simulations, in which the full set of MHD equations is solved. These methods have been used to study the wave evolution in simple 2D equilibria (e.g., Malara et al. 1996). For a 3D configuration, the limited number of numerical degrees of freedom makes inaccessible the high Reynolds/Lundquist numbers which are typical of the coronal plasma. This fact limits the applicability of direct numerical simulations to model wave dissipation in a complex 3D magnetic structure.

2.1 Alfvén wave evolution in the WKB limit

An alternative approach to the problem of Alfvén wave propagation and dissipation in 3D equilibria was proposed by Similon & Sudan (1989), and subsequently investigated by Petkaki et al. (1998) and by Malara et al. (2000, 2003). Here we use this same method, which is briefly summarized in the following. Some assumptions are required: (i) Alfvénic perturbations are supposed to have a small amplitude: $\delta v/c_{\rm A0} \ll 1$, where $\delta v$ is the typical velocity fluctuation and $c_{\rm A0}$ is the background Alfvén velocity; (ii) perturbations have a small wavelength $\lambda /L \ll 1$, where $\lambda$ is the wavelength and L is the scale of variation of the equilibrium magnetic field. Under the above assumptions the MHD equations are linearized and a WKB expansion is applied. The Alfvénic perturbation is divided into a large number of packets, each occupying an infinitesimal volume; thus, each packet is characterized by a position $\vec{x}=\vec{x}(t)$, a wavevector $\vec{k}=\vec{k}(t)$ and an energy e=e(t). When the packet propagates inside the equilibrium structure, these quantities evolve in time. The time evolution of each packet is calculated by solving the equations:

 \begin{displaymath}
{{\rm d} x_i^{\pm\rm A} \over {\rm d}t}= \pm c_{\rm A0} b_{\rm0i}
\end{displaymath} (5)


 \begin{displaymath}
{{\rm d}k_i^{\pm\rm A} \over {\rm d}t}=
\mp {\partial (c_{\rm A0} b_{0n}) \over \partial x_i} k_n^{\pm\rm A}
\end{displaymath} (6)


 \begin{displaymath}
{{\rm d}{\rm e}^{\pm\rm A} \over {\rm d}t}= -{2{k^{\pm\rm A}}^2 \over S} {\rm e}^{\pm\rm A}
\end{displaymath} (7)

where $\vec{b}_0=\vec{B}_0/B_0$ is the unit vector in the direction of the equilibrium magnetic field $\vec{B}_0$; $\rho_0$ is the background density; $c_{\rm A0}=B_0/(4\pi \rho_0)^{1/2}$ is the Alfvén velocity; all these quantities depends on the position $\vec{x}$. The Reynolds number S is defined by $S^{-1}=(\eta_N / {\rho}_0 + \lambda_N)/2$. The sign in the upper index identifies the propagation direction of the packet: +A in the direction of $\vec{B}_0$ and -A in the opposite direction. The same normalization of Eqs. (2)-(4) has been used.

Equations (5)-(7) have been derived by Petkaki et al. (1998) for an incompressible plasma. Recently, Malara et al. (2003) starting from Eqs. (2)-(4) showed that Eqs. (5)-(7) describe the evolution of Alfvénic packets also in a compressible cold plasma. Equations (5)-(7) indicate that: (i) each Alfvénic packet propagates at the Alfvén velocity $\vec{c}_{\rm A0}$ along the magnetic field lines; (ii) wavevector changes are due to spatial variations of $\vec{c}_{\rm A0}$, which distort the packet; (iii) energy dissipation takes place with a rate proportional to k2 and to 1/S.

Considering simple 2D configurations of the equilibrium structure, where the Alfvén velocity $\vec{c}_{\rm A0}$ is directed along a given direction (say, z) and varies along another perpendicular direction (say, x), Eqs. (5)-(7) can be exactly solved (Petkaki et al. 1998): in this case it is seen that the wavevector k asymptotically increases proportional to time: $k \propto t$, while the typical dissipation time $t_{\rm d}$ scales according to the law $t_{\rm d} \propto S^{1/3}$. These features are typical of phase-mixing of the Alfvén waves (Heyvaerts & Priest 1983). Since in the corona the dissipative coefficients (relevant viscosity and resistivity) are very low, the dissipation time due to phase-mixing is exceedingly long; i.e., Alfvén waves would be dissipated after traveling over lengths much larger than the size of coronal structures.

For more complex 3D equilibria, where the Alfvén velocity $\vec{c}_{\rm A0}$ depends on the three spatial coordinates, regions may exist where magnetic field lines are chaotic, i.e., nearby magnetic lines exponentially move apart. Equation (5) indicates that Alfvénic perturbations propagate locally following magnetic lines; then, in a chaotic region Alfvénic packets are rapidly stretched, and the wavevector exponentially increases in time. Equations (5)-(7) describe this behaviour which leads to the scaling law (1) of dissipation time, typical of fast dissipation (Petkaki et al. 1998). Thus, even for very low values of the dissipative coefficients (viscosity and/or resistivity), relatively short dissipation times are achieved, contrary to the phase-mixing-driven dissipation. The two phenomena (phase-mixing and 3D fast dissipation) can co-exist in different regions within the same equilibrium structure (Malara et al. 2000).

2.2 Equilibrium magnetic field

In the above-cited papers the main interest was to investigate the properties of Alfvén wave propagation and dissipation in complex 3D magnetic field. Thus, the explicit form used for the background structure was not particularly representative of the magnetic field in the corona. In the present paper, we want to study Alfvén wave dissipation in an equilibrium magnetic field which is a more realistic representation of the coronal magnetic field. In particular, we will use a force-free magnetic field which was first introduced by Nakagawa & Raadu (1972) within the problem of extrapolating the magnetic field in the corona starting from measurements in the photosphere. For completeness, in the following we briefly describe the derivation of this magnetic field, which requires some assumptions:

(i)
the equilibrium magnetic field $\vec{B}_0(\vec{x})$ satisfies the force-free condition

 \begin{displaymath}
\nabla \times \vec{B}_0 = \alpha \vec{B}_0
\end{displaymath} (8)

with $\alpha$ constant. The force-free condition is a natural requirement, since it is well satisfied by an equilibrium magnetic field in a low-$\beta$ plasma, like in corona. The condition $\alpha = {\rm const.}$ is a much stronger assumption; on the other hand, it allows us to simplify the derivation of the equilibrium magnetic field;

(ii)
the curvature due to the spherical geometry of the corona is neglected. This assumption is valid if relatively small portions of the corona are considered. This allows us to use a Cartesian reference frame, in which the xy plane represents the base of the corona, while the z axis is in the upward vertical direction;

(iii)
we assume statistical homogeneity along the horizontal directions x and y. This implies that no dominant structures (e.g., a magnetic dipole) are present. On the contrary, the magnetic field distribution is characterized by several regions of different polarities randomly distributed at various spatial scales. This situation can be found in quiet-Sun regions, and it is illustrated, for instance, in the magnetogram plotted in Fig. 1 of Schrijver et al. (1997). In order to describe this kind of configuration, we assume that along the horizontal x and y directions the equilibrium magnetic field has a characteristic scale of variation L0 and it is periodic with a periodicity length L. Statistical homogeneity is reproduced when $L_0 \ll L$ (Pommois et al. 1998); in particular, we used L = 4   L0. Comparing with the magnetogram of Schrijver at al. (1997) (their Fig. 1), the characteristic length for the magnetic field variation can be estimated as $L_0 \sim 2.5 {-} 3 \times 10^4$ km, thus, the unit length will be $L\sim 1{-}1.2 \times 10^5$ km. When using dimensionless variables, the periodicity length is equal to 1, while the magnetic field characteristic scale is l0 = L0/L = 1/4;

(iv)
following Nakagawa & Raadu (1972), we assume that the magnetic field vanishes in the limit $z \rightarrow +\infty$. This implies that the magnetic energy in a single periodicity box is finite (Alissandrakis 1981).
Using the above assumptions, a form for the equilibrium magnetic field can be derived. Using the periodicity condition, the magnetic field can be expanded in Fourier series in the x and y directions:

 \begin{displaymath}
B_{\rm0i}(x,y,z) = \sum_{\kappa_x,\kappa_y} {\hat B}_{\rm0i}...
...a_y,z)
\exp \left [ {\rm i}(\kappa_x x + \kappa_y y) \right ]
\end{displaymath} (9)

where the wavenumbers are $\kappa_x = (2\pi /L) n_y$ and $\kappa_y = (2\pi /L) n_y$, with nx and ny integers. The dependence on z is included in the Fourier coefficients ${\hat B}_{\rm0i}$. Inserting the expansion (9) into the force-free condition (8), we derive the equations:


 \begin{displaymath}
{{\rm d}^2 {\hat B}_\parallel \over {\rm d}z^2} =
\left( \kappa^2 - \alpha^2 \right) {\hat B}_\parallel
\end{displaymath} (10)


 \begin{displaymath}
{{\rm d} {\hat B}_\perp \over {\rm d}z} = - \alpha {\hat B}_\parallel
\end{displaymath} (11)


 \begin{displaymath}
{\hat B}_{0z}={\rm i \over \alpha} \left(
\kappa_x {\hat B}_y - \kappa_y {\hat B}_x \right)
\end{displaymath} (12)

where $\kappa =(\kappa_x^2 + \kappa_y^2)^{1/2}$, while ${\hat B}_\parallel$ and ${\hat B}_\perp$ are the components in the horizontal direction, parallel and perpendicular to the wavevector $(\kappa_x,\kappa_y)$, respectively. The only solution of Eq. (10) which satisfies the above assumption (iv) is in the form:

 \begin{displaymath}
{\hat B}_\parallel = b(\kappa_x,\kappa_y) \exp \left(-\sqrt{\kappa^2-\alpha^2} z\right)
\end{displaymath} (13)

where $\kappa > \alpha$. This condition implies that $\alpha$ is a lower bound for the horizontal wavevector $\kappa$. We will assume that $\kappa$ varies in the range $\left [ \kappa_0 , \kappa_{\rm max} \right ]$, where $\kappa_0=2\pi /l_0$ corresponds to the typical scale of $\vec{B}_0$, while $\kappa_{\rm max}$ corresponds to the smallest length in the equilibrium structure. Then, the largest length in the equilibrium magnetic field structure corresponds to the characteristic horizontal length l0. In particular, $\alpha < \kappa_0$.

Using the expression (13) the Eq. (11) is solved, and the solution is used to write the equilibrium magnetic field components in the form:

 
    $\displaystyle B_{0x}(x,y,z) = \sum_{\kappa = \kappa_0}^{\kappa_{\rm max}}
\left...
... \left[ {\rm i}(\kappa_x x + \kappa_y y) - \sqrt{\kappa^2 - \alpha^2}z \right ]$  
    $\displaystyle B_{0y}(x,y,z) = \sum_{\kappa = \kappa_0}^{\kappa_{\rm max}}
\left...
... \left[ {\rm i}(\kappa_x x + \kappa_y y) - \sqrt{\kappa^2 - \alpha^2}z \right ]$ (14)
    $\displaystyle B_{0z}(x,y,z) = \sum_{\kappa = \kappa_0}^{\kappa_{\rm max}}
{\rm ...
... \left[ {\rm i}(\kappa_x x + \kappa_y y) - \sqrt{\kappa^2 - \alpha^2}z \right].$  

The specific form of the magnetic field is determined by the complex Fourier amplitudes $b(\kappa_x ,\kappa_y)$. The coefficient $1/\sqrt{\kappa^2 - \alpha^2}$ in the expressions (14) diverges in the limit $\kappa \rightarrow \alpha$. Though $\kappa$ is strictly larger than $\alpha$, it is useful to eliminate such a divergence by simply re-defining the amplitude of the Fourier components. In particular, we define:

 \begin{displaymath}
a(\kappa_x,\kappa_y) = {\kappa \over \sqrt{\kappa^2 - \alpha^2}}
b(\kappa_x,\kappa_y)
\end{displaymath} (15)

where $ a(\kappa_x,\kappa_y)=\vert a (\kappa_x,\kappa_y)\vert
\left [ \cos \phi (\kappa_x,\kappa_y) +
{\rm i} \sin \phi (\kappa_x,\kappa_y) \right ]$ is the new complex amplitude. The reality condition $B_{\rm0i}(x,y,z) = B_{\rm0i}^*(x,y,z)$ (where the asterisk indicate complex conjugate) implies:

 \begin{displaymath}
a(-\kappa_x,-\kappa_y) = -a^*(\kappa_x,\kappa_y).
\end{displaymath} (16)

Using this condition, and the definition (15), after some algebra the equilibrium magnetic field components can be put in the following form:
 
    $\displaystyle B_{0x}(x,y,z) = 2 \sum_{\kappa_0 \le \kappa \le \kappa_{\rm max},...
... \phi(\kappa_x,\kappa_y) \right]
\exp \left(-\sqrt{\kappa^2 - \alpha^2}z\right)$  
    $\displaystyle B_{0y}(x,y,z) = 2 \sum_{\kappa_0 \le \kappa \le \kappa_{\rm max},...
... \phi(\kappa_x,\kappa_y) \right]
\exp \left(-\sqrt{\kappa^2 - \alpha^2}z\right)$ (17)
    $\displaystyle B_{0z}(x,y,z) = - 2 \sum_{\kappa_0 \le \kappa \le \kappa_{\rm max...
... \phi(\kappa_x,\kappa_y) \right]
\exp \left(-\sqrt{\kappa^2 - \alpha^2}z\right)$  

where only real quantities appear and the sums are extended only on half of the $\kappa_x \kappa_y$ plane. The expressions (17) indicate that the equilibrium magnetic field is obtained as a superposition of several harmonics, each characterized by an horizontal wavevector $(\kappa_x,\kappa_y)$ which defines a corresponding scale $2\pi /\kappa$. Each harmonic exponentially decreases with increasing the altitude z with a rate $\gamma_z=-\sqrt{\kappa^2-\alpha^2}$; thus, the contribution of harmonics with small horizontal wavelengths is limited to low altitudes, while at higher altitudes only the longest wavelength harmonics are present ( $\kappa \sim \alpha$).

The expression (17) depend on some parameters: namely, the coefficient $\alpha$; the minimum and maximum wavenumber $\kappa_0$ and $\kappa_{\rm max}$; the amplitudes $\vert a(\kappa_x,\kappa_y)\vert$ and the phases $\phi(\kappa_x,\kappa_y)$ of harmonics. In the context of extrapolating the coronal magnetic field from measures in the photosphere, the above parameters should be determined using appropriate boundary conditions at the coronal base z=0. In our approach we are not interested in representing any specific observed configuration; on the contrary, the equilibrium magnetic field should represent a generic magnetic structure above a quiet Sun region. In fact, we are mainly interested in studying the general properties of Alfvén wave dissipation, which should not depend on the details of the equilibrium structure. For this reason, we have chosen the above parameters in the following way.

(i)
We assume that harmonic amplitudes $\vert a(\kappa_x,\kappa_y)\vert$ depend only on the modulus $\kappa$ of the horizontal wavevector, according to a power law:

 \begin{displaymath}
\vert a(\kappa_x,\kappa_y)\vert = c \kappa^{-\mu} , \;\;\; {\rm for} \;\;
\kappa_0 \le \kappa \le \kappa_{\rm max}
\end{displaymath} (18)

where $\mu$ is the spectral index and c is a normalization constant. In order to choose a value for $\mu$, we consider the magnetic energy at the coronal base z=0 in a periodicity box, which can be expressed in terms of the magnetic field Fourier components by using the Parseval theorem:

 \begin{displaymath}
E_{\rm M}(z=0)=\int_0^L \int_0^L {B_0^2(x,y,z=0) \over 8 \pi...
...ppa_y}
\vert\vec{\hat {B}}_{0}(\kappa_x,\kappa_y,z=0)\vert^2.
\end{displaymath} (19)

In the limit of large L, the sum can be approximated by an integral, according to the prescription

\begin{displaymath}\sum_{\kappa_x,\kappa_y} \longrightarrow
\left({L\over 2\pi }\right)^2 \int {\rm d}\kappa_x {\rm d}\kappa_y.
\end{displaymath}

Thus, Eq. (19) gives

 \begin{displaymath}
E_{\rm M}(z=0) \simeq {L^4 \over 32 \pi^3} \int {\rm d}\kapp...
...rt{\vec{\hat B}}_{0}(\kappa,z=0)\vert^2
\kappa ~ {\rm d}\kappa
\end{displaymath} (20)

where we used the assumption that ${\vec{\hat B}}_{0}$ depends on the modulus $\kappa$ of the wavevector. If the magnetic field structure is the results of a turbulent process taking place beneath the coronal base, we can assume that

 \begin{displaymath}
E_{\rm M}(z=0) = \int \varepsilon (\kappa ) {\rm d}\kappa \;...
...rm with} \;\;\;
\varepsilon (\kappa ) \propto \kappa^{-\beta}.
\end{displaymath} (21)

Comparing Eqs. (20) and (21) we see that

 \begin{displaymath}
\vert a(\kappa_x,\kappa_y)\vert \propto \vert\vec{{\hat B}}_{0}(\kappa,z=0)\vert
\propto \kappa^{-{\beta+1 \over 2}}.
\end{displaymath} (22)

Using the Kolmogorov spectral index $\beta=5/3$ of the spectral energy density $\varepsilon (\kappa )$, we obtain the spectral index $\mu = (\beta + 1)/2 = 4/3$ of the amplitude $\vert a(\kappa)\vert$. The normalization constant c in Eq. (18) is determined requiring that $\langle B_0^2 \rangle_{xy} (z=0)=1$. Thus, the unit magnetic field in dimensionless units corresponds to the rms magnetic field at the base.

(ii)
Due to the form (18), the maximum of energy is at the largest wavelengths $l_0=2\pi / \kappa_0$. As discussed above, in order to reproduce statistical homogeneity we used l0=L0/L=1/4, corresponding to $\kappa_0 =4(2 \pi )$. The maximum wavenuber is $\kappa_{\rm max}=10(2\pi )$. With these values, each component of the equilibrium magnetic field is a superposition of 133 harmonics, which gives a complex magnetic field but computational times not exceedingly long.

(iii)
The parameter $\alpha$ determines the current density $\vec{j}_0=(c/4\pi ) \nabla \times \vec{B}_0$ to the equilibrium magnetic field $\vec{B}_0$ ratio. Moreover, in our model the value of $\alpha$ is bounded by the minimum horizontal wavevector: $\alpha < \kappa_0$. We considered two different values for the parameter $\alpha$: a) $\alpha =0$, corresponding to a potential magnetic field ( $\vec{j}_0=0$); b) $\alpha = 7\pi $, which is close to the upper bound $\kappa_0$ and gives a nonvanishing current $\vec{j}_0$.

(iv)
Since $\vec{B}_0$ does not represent any particular magnetic field, the phases $\left \{ \phi(\kappa_x,\kappa_y) \right \}$ of harmonics have been randomly chosen in the interval $[0, 2 \pi [$.

In Fig. 1 a contour plot of the vertical component B0z(x,y,0) at the base z=0 in the periodicity domain is represented, in the case $\alpha = 7\pi $. Darker or clearer tones correspond to positive or negative values of B0z, respectively. It can be seen that B0z is randomly distributed, with structures at different length scales. This is reminiscent of a portion of the magnetogram in a quiet Sun region displayed in Schrijver et al. (1997). In Fig. 1 several magnetic lines are represented; they have been obtained integrating the Eq. (5). Each magnetic line starts from a point randomly chosen at the base z=0 and it ends at another location at z=0, linking two points with a different magnetic polarity (sign of B0z). Since the magnetic field intensity decreases with the height z, the density of field lines must decrease, too: thus, most of lines remain at low altitudes and very few reach higher values of z. Magnetic lines which are limited to low altitudes are mostly in the form of an arc which crosses a neutral line at the base, defined by B0z(x,y,0)=0. Magnetic lines which extends to higher altitudes have a more complex behavior: they are more twisted and do not have a definite shape.


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{2272fig1.eps}
\end{figure} Figure 1: A contour plot of the vertical component B0z of the equilibrium magnetic field, along with the trajectories of 70 magnetic lines. The stating point of each magnetic line is randomly chosen inside the periodicity domain.

The overall structure of the magnetic field represented in Fig. 1 is quite complex. Dissipation properties of Alfvén waves are directly related to the topological complexity of the magnetic field (Petkaki et al. 1998; Malara et al. 2000). In order to illustrate some features of the topology of the equilibrium magnetic field we drew some magnetic surfaces in form of flux tubes. These magnetic surfaces are formed by several magnetic lines. Each line starts from a point at the coronal base located in a small circle; in other words, the initial section of the flux tube is circular, and it is located within a region with the same magnetic polarity. We calculated several of these flux tubes, finding that they can have two kinds of structures, which we call a "regular'' flux tube and a "singular'' flux tube, respectively. In Fig. 2 (upper panel) an example of a regular flux tube is represented: it is a sort of cylinder linking two regions of different polarities located at the base z=0. The tube section, initially circular, changes along the tube axis, being larger at higher altitudes, but it remains a closed line at any location along the tube. Magnetic lines wrap around the magnetic surface and their twist is related to the presence of a nonvanishing current. In Fig. 2 (lower panel) an example of a singular flux tube is represented. At variance with the previous situation, in this case the magnetic surface divides into several sheets, each following a different path: all lines forming this surface start from the same region with a given polarity, but they end in four different regions with the opposite polarity. The singularity is in the behaviour of magnetic lines around locations where the surface separates in different sheets: two initially infinitely close magnetic lines can either remain close to each other (if they belong to the same sheet) or they can follow two completely different trajectories (if they belong to two different sheets). This behaviour is found also into phase space trajectories of chaotic dynamical systems. Considering the overall 3D structure of the magnetic field, this separation takes place onto separatrix surfaces: magnetic lines which are on the two sides of the separatrix will follow diverging paths. Although we did not carry out any systematic study of these features in the magnetic structure, we expect that it contains several separatrices. A detailed study of the general properties of separatrices in 3D fields can be found in Priest & Titov (1996).


  \begin{figure}
\par\includegraphics[width=9.4cm,clip]{2272fig2.eps}\par\vspace*{2mm}
\includegraphics[width=8.9cm,clip]{2272fig3.eps}
\end{figure} Figure 2: A regular ( upper panel) and a singular ( lower panel) flux tube.

2.3 Single packet evolution

Divergence of nearby magnetic lines is related to the increase of the wavevector of Alfvénic disturbances, which propagate locally following magnetic lines. The volume of an Alfvénic packet is approximately constant (this is strictly verified for uniform $\rho_0$, but it holds also when $\rho_0$ varies on a sufficiently large scale). Thus, when a packet goes through locations where magnetic lines diverge it is stretched and squeezed. This results in an increase of the wavevector. Another independent mechanism which determines wavevector increase is related to the exponential dependence of $\vec{B}_0$ on the z coordinate. A simple configuration where the Alfvén velocity varies exponentially with the altitude was considered by Ruderman et al. (1998). These authors showed that in such a case the damping of an Alfvénic perturbation takes place exponentially with the distance along a field line. This is due to the fact that two points in the perturbation which propagate along two nearby magnetic lines move apart exponentially, thus resulting in an exponential increase of the wavevector. This same effect is present also in the more complex structure we are considering.


  \begin{figure}
\par\includegraphics[width=6cm,clip]{2272fig4.eps}\par\vspace*{2m...
...ar\vspace*{2mm}
\includegraphics[width=6.5cm,clip]{2272fig6.eps}
\end{figure} Figure 3: The trajectory of a packet ( upper panel); the wavevector k ( middle panel) and the energy e ( lower panel) of the packet, as functions of time t.

A description of the dynamics of a single packet can be obtained by integrating the Eqs. (5)-(7). As an example, the time evolution of the position $\vec{x}$, the wavevector $k=\vert\vec{k}\vert$ and the energy of a packet is shown in Figs. 3. The initial condition is the following:

(i)
the initial position $\vec{x}_0=(x_0,y_0,0)$ is located at the base;

(ii)
the initial wavevector is $\vec{k}_0=(20\pi ,20\pi ,20\pi )$, giving $k_0 \simeq 35\pi \gg \kappa_0$, i.e., large enough to satisfy the WKB approximation;

(iii)
the initial energy is e0=1; this value is not relevant to determine the time evolution: since the Eq. (7) is linear with respect to e, the physically relevant quantity is the ratio e(t)/e0. The sign in Eqs. (5) and (6) determining whether the packet propagates parallel or antiparallel to $\vec{B}_0$ is chosen so that the packet initially propagates in the upward direction. The value of the Reynolds number is $S=2 \times 10^5$.
Following a magnetic line, the trajectory of the packet starts and ends at the base z=0; the time evolution of Fig. 3 corresponds to this path. The wavevector k increases roughly exponentially in time (Fig. 3, middle panel), in accordance with the previous considerations; this characterizes the fast dissipation mechanism (Petkaki et al. 1998). The energy of the packet is fully dissipated (Fig. 3, lower panel). A time evolution similar to that shown in Fig. 3 has been found for many other packets with different initial positions, especially for those whose trajectory reaches large enough altitudes. In such cases there is enough time for the energy to dissipate before the packet reaches the base. The dissipation time $t_{\rm d}$ for a packet is defined by $e(t_{\rm d})/e_0=\exp (-1)$; $t_{\rm d}$ increases with increasing the Reynolds number S. In particular, $t_{\rm d}$ is proportional to $\log S$ if the wavevector increases exponentially in time (Petkaki et al. 1998). The dissipation time of the packet of Fig. 3 is represented in Fig. 4 for different values of S: the rough exponential growth of k (Fig. 3) gives $t_{\rm d}$ approximately proportional to $\log S$. This is typical of fast dissipation.


  \begin{figure}
\par\includegraphics[width=6.5cm,clip]{2272fig7.eps}
\end{figure} Figure 4: The dissipation time $t_{\rm d}$ of the same packet as in Figs. 6, as a function of the Reynolds number S (dots). The full line indicates a scale law $t_{\rm d} \propto \log S$.

3 Evolution of Alfvénic perturbations

Packets whose trajectory is confined at lower altitudes has a time evolution different from that displayed in Fig. 3. If the magnetic line where the packet propagates is too short, only a fraction of the initial energy is dissipated before the packet reaches the base. Moreover, when the Reynolds number is increased the dissipation rate becomes smaller; in such a case also packets reaching higher altitudes do not have enough time to dissipate all their initial energy. For these reasons it is important to establish what happens when a packet reaches the coronal base. Since large vertical gradients of the background quantities are present beneath the base of the corona, we will model the coronal base as a discontinuity surface. In this case, an Alfvén wave which reaches the base z=0 will give origin to transmitted and reflected waves, both Alfvénic and magnetosonic. In particular, the subsequent time evolution of the reflected Alfvénic packet can be described by using the same Eqs. (5)-(7). Reflection and transmission coefficients determine how the energy of the incident packet is redistributed into reflected and transmitted packets. A simple model of the discontinuity at the coronal base, as well as the detailed calculation of the corresponding transmission and reflection coefficients are given in the Appendix.

3.1 Evaluation of the reflected and transmitted packets parameters

The value of reflection $r_{\rm A}$, $r_{\rm M}$ and transmission $t_{\rm A}$, $t_{\rm M}$ coefficients of Alfvénic and magnetosonic waves calculated in the Appendix (Eqs. (A.103)-(A.106)) depends on some parameters:

(i)
the angle $\phi $ between the magnetic field $\vec{B}_0$ and the horizontal direction. This angle has different values at different positions at the base z=0;

(ii)
the ratio $r=(\rho_0^+/\rho_0^-)^{1/2}=c_{\rm A}^-/c_{\rm A}^+$, where the upper index + or - refers to the value just above or below the discontinuity, respectively. The density ratio between the corona and the chromosphere has been estimated as $\rho_0^+/\rho_0^- \sim 10^{-3}$ (see Mariska 1992). More recently, Ofman (2002) in a 1D model of Alfvén wave energy leakage from the corona used the value $ c_{\rm A}^-/c_{\rm A}^+=10^{-2}$. Then, in our model, we considered the values r=0.01 and r=0.032;

(iii)
the wavevector $\vec{k}^{\rm IA}$ of the incident Alfvén wave. It can be verified that the reflection and transmission coefficients do not depend on the modulus but only on the orientation of the wavevector $\vec{k}^{\rm IA}$. Moreover, they do not change under the transformation $\vec{k}^{\rm IA} \rightarrow -\vec{k}^{\rm IA}$.
Sometimes the reflected and/or the transmitted magnetosonic wave is evanescent (see the Appendix). In such a case, the corresponding reflection/transmission coefficient vanishes.

As an example, in Fig. 5 we plotted the reflection and transmission coefficients as functions of the angle $\phi $, for r=0.01 (upper panel) and r=0.032 (lower panel), with a given value of the wavevector $\vec{k}^{\rm IA}$. The value of $\vec{k}^{\rm IA}$ used in Figs. 5 has been chosen in the following way. We have verified that during the Alfvénic packet propagation the wavevector component parallel to $\vec{B}_0$ remains roughly constant, while the perpendicular components tend to increase in time. Then, we considered the parallel component $k_{\vert\vert}^{\rm IA}=k_Y \cos \phi + k_Z^{\rm IA} \sin \phi$ and the two perpendicular components $k_{1\perp}^{\rm IA}=k_X$ and $k_{2\perp}^{\rm IA}= k_Y \sin \phi - k_Z^{\rm IA} \cos \phi$. In Figs. 5 we chose the wavevector $\vec{k}^{\rm IA}$ so that $k_{\vert\vert}^{\rm IA}/ k_{1\perp }^{\rm IA}=0.1$ and $k_{\vert\vert}^{\rm IA}/ k_{2\perp }^{\rm IA}=0.1$. At r=0.01 the Alfvén wave reflection coefficient is close to unity: $r_{\rm A} \simeq 0.96$, and it is largely independent of the angle $\phi $; transmission coefficients $t_{\rm A}$ and $t_{\rm M}$ are much smaller, both being at most $\simeq$0.04; the reflected magnetosonic wave is evanescent ( $r_{\rm M}=0$) for any angle $\phi $. At r=0.032 the reflection coefficients is slightly lower: $r_{\rm A} \simeq 0.88$, and again mostly independent of $\phi $; both transmission coefficients are of the order $\sim $0.1; the maximum transmission of Alfvén waves is at $\phi \simeq 0$ and $\phi \simeq 90$ ( $t_{\rm A} \simeq 0.12$), while the maximum transmission of magnetosonic waves is at $\phi \simeq 60$ ( $t_{\rm M} \simeq 0.1$); the reflected magnetosonic wave is evanescent for any angle $\phi $. The wavevector $\vec{k}^{\rm IA}$ used in Figs. 5 is only an example (in the model, the actual value of $\vec{k}^{\rm IA}$ is determined by the packet evolution). However, these figures shows that for sufficiently small values of the parameter r a large fraction of the energy of the incident Alfvénic packet is transferred to the reflected Alfvénic packet, while the remaining part is essentially lost in waves (Alfvénic and magnetosonic) transmitted below the coronal base.


  \begin{figure}
\par\includegraphics[width=6.6cm,clip]{2272fig8.eps}\par\vspace*{4mm}
\includegraphics[width=6.6cm,clip]{2272fig9.eps}
\end{figure} Figure 5: The reflection coefficient $r_{\rm A}$ (thick line), the transmission coefficients $t_{\rm A}$ (dashed line) and $t_{\rm M}$ (thin line), as functions of the angle $\phi $, calculated for $k_{\vert\vert}^{\rm IA}/ k_{1\perp }^{\rm IA}=0.1$ and $k_{\vert\vert}^{\rm IA}/ k_{2\perp }^{\rm IA}=0.1$, and for r=0.01 ( upper panel) and r=0.032 ( lower panel).

Energy conservation requires that

 \begin{displaymath}
r_{\rm A} + t_{\rm A} + r_{\rm M} + t_{\rm M} = 1.
\end{displaymath} (23)

We checked that this condition is verified by our calculation of reflection and transmission coefficients.

The development described in the Appendix allowed us to determine all the parameters (wavevector and energy) of the reflected and transmitted wavepackets, as functions of the parameters of the incident Alfvénic packet. The procedure can be summarized as follows. We follow the time evolution of an Alfvénic packet integrating the Eqs. (5)-(7),

(i)
when the packet reaches the level z=0 coming from above, we consider its wavevector $\vec{k}^{\rm IA}=k_x \vec{e}_x + k_y \vec{e}_y + {k_z}^{\rm IA} \vec{e}_z $ and its energy $e^{\rm IA}$;

(ii)
the components of the wavevector are given with respect to the xyz reference frame. Thus, first we calculate the components of $\vec{k}^{\rm IA}=k_X \vec{e}_X + k_Y \vec{e}_Y + {k_Z}^{\rm IA} \vec{e}_Z $ with respect to the XYZ reference frame (see the Appendix, Fig. A.2); the direction of the X and Y axes and the $\phi $ angle are obtained calculating the equilibrium magnetic $\vec{B}_0$ at the packet location onto the plane z=0;

(iii)
the vertical (Z) components of the wavevectors of the reflected and transmitted Alfvénic packets are calculated using the Eqs. (A.57) and (A.60). The sign of the RHS of Eqs. (A.64) and (A.73) determines whether the magnetosonic packets are propagating or evanescent; in the former case, to calculate the components $\vec{k}_z^{\rm RM}$ and $\vec{k}_z^{\rm RM}$ Eqs. (A.67) and (A.75) are used, while in the latter case we use Eqs. (A.69) and (A.77). The horizontal (x and y) components of all reflected and transmitted packets are the same as those of the incident wave;

(iv)
the wave amplitude ratios are calculated by numerically solving the system of complex linear Eqs. (A.79)-(A.82). Such ratios are used in Eqs. (A.99) and (A.100) to determine the energy density ratios $E^{\rm RA}/E^{\rm IA}$, $E^{\rm TA}/E^{\rm IA}$, $E^{\rm RM}/E^{\rm IA}$, and $E^{\rm TM}/E^{\rm IA}$, which are used to calculate the reflection and transmission coefficients (Eqs. A.103)-(A.106));

(v)
the energy of reflected Alfvénic $e^{\rm RA}$, reflected magnetosonic $e^{\rm RM}$, transmitted Alfvénic $e^{\rm TA}$, and transmitted magnetosonic $e^{\rm TM}$ packets are determined by the equations

 \begin{displaymath}
e^{\rm RA}=r_{\rm A} e^{\rm IA} , \;\;\; e^{\rm RM}=r_{\rm M...
...t_{\rm A} e^{\rm IA} , \;\;\; e^{\rm TM}=t_{\rm M} e^{\rm IA}
\end{displaymath} (24)

where $e^{\rm IA}$ is the energy of the incident Alfvénic packet.

  
3.2 Energy balance of Alfvénic perturbations

We considered an Alfvénic perturbation formed by a large number of packets. At each time t the energy E(t) of the perturbation is given by the sum of the energy e(t) of each packet. The time evolution of each packet inside the equilibrium magnetic field $\vec{B}_0$ is determined by integrating the Eqs. (5)-(7), starting from the initial conditions specified in Sect. 2.3. In particular, the packet initial position $\vec{x}_0=(x_0,y_0,0)$ is located at the base, and it is randomly chosen within the periodicity domain. If the number $N_{\rm p}$ of packets is large enough, they can represent the evolution of the whole Alfvénic perturbation.

In particular we used $N_{\rm p}=1000$, and we verified that increasing $N_{\rm p}$ the results do not significantly change. Each time a packet reaches the base z=0 it is replaced by a reflected Alfvénic packet, which propagates back. The wavevector and the energy of the reflected packet are calculated as described in Sect. 3.2. At the same time, the energies of the magnetosonic reflected perturbation and Alfvénic and magnetosonic transmitted perturbations are stored. The evolution of each packet is followed until its energy e(t) drops under the value 10-6e0, where e0=e(t=0); at that time the initial energy e0 has been completely dissipated and/or converted in energy of transmitted or magnetosonic reflected perturbations. The run ends when all the packets have reached this condition.

We performed the above calculation for different values of the parameters of the problem. In particular, we considered two values for $\alpha$: namely, $\alpha =0$ (potential magnetic field) and $\alpha = 7\pi $ (nonvanishing equilibrium current). The Reynolds number varies from $S=2\times 10^4$ up to $S=2\times 10^{10}$. We considered two values for the ratio r: r=0.01 and r=0.032. The energy balance is represented by the dissipated energy $E_{\rm diss}$ of Alfvénic perturbations, the energies $ E_{\rm t}^{\rm A} $ and $ E_{\rm t}^{\rm M} $ lost in transmitted Alfvénic and magnetosonic perturbations, respectively, and the energy $E_{\rm r}^{\rm M} $ of reflected magnetosonic perturbations. These energies satisfy the balance equation

 \begin{displaymath}
E_{\rm diss}+ E_{\rm t}^{\rm A} + E_{\rm t}^{\rm M} + E_{\rm r}^{\rm M} =E_0
\end{displaymath} (25)

where E0 is the initial perturbation energy. In Figs. 6 we plotted the ratios $E_{\rm t}^{\rm A}/E_0$, $E_{\rm t}^{\rm M}/E_0$, and $E_{\rm r}^{\rm M}/E_0$, calculated by the above procedure, for different values of the parameters $\alpha$, S and r. The dissipated to initial energy ratio $E_{\rm diss}/E_0$ is derived by the Eq. (25) and it is also plotted.


  \begin{figure}
\par\includegraphics[width=8cm,clip]{2272fi10.eps}\par\vspace*{2mm}
\includegraphics[width=8cm,clip]{2272fi11.eps}
\end{figure} Figure 6: The energy ratios $E_{\rm diss}/E_0$ (full line), $E_{\rm t}^{\rm A}/E_0$ (dashed line), $E_{\rm t}^{\rm M}/E_0$ (dotted line), and $E_{\rm r}^{\rm M}/E_0$ (dash-dotted line) are shown as functions of the Reynolds number S, for $\alpha = 7\pi $ (no symbols) and for $\alpha =0$ (black dots). The upper and lower panels refer to the case r=0.01 and r=0.032, respectively.

The upper panel of Fig. 6 shows the results obtained for r=0.01. For low values of the Reynolds number $S=2\times 10^4$ the ratio $E_{\rm diss}/E_0$ is very close to unity. This indicates that the initial perturbation energy E0 is almost completely dissipated inside the magnetic equilibrium structure. On the contrary, the fraction of energy lost in waves (both Alfvénic and magnetosonic) transmitted through the discontinuity at the base is very low (less than 5%), while the fraction of energy of reflected magnetosonic waves is negligible ($\sim $ $2{-}3 \times 10^{-3}$). With increasing S the dissipated energy decreases, while the energy lost in transmitted waves increases. This behaviour is the result of the balance between the dissipation time $t_{\rm d}$ of a single packet and the crossing time $t_{\rm c}$, i.e., the time the packet takes to propagate all along a magnetic line connecting two points of the base z=0. The dissipation time $t_{\rm d}$ can be defined as the time the wavevector k of a packet takes to reach a typical value $k_{\rm d}=(S/2)^{1/2}$. For $k\sim k_{\rm d}$the time scale of energy dissipation is $\sim $1 (see Eq. (7)); then, energy dissipation for the packet takes place at $t\sim t_{\rm d}$. The dissipation time $t_{\rm d}$ increases roughly proportional to $\log S$ (see Fig. 4; and Malara et al. 2000). On the contrary, the crossing time $t_{\rm c}$ for a given packet is independent of S. For sufficiently small values of S, $t_{\rm d} < t_{\rm c}$. In this case, most of the packet energy is dissipated during the first trip along the magnetic line. Then, only a small fraction of the initial energy is lost in transmitted waves during the reflection process. With increasing S, the dissipation time becomes increasingly larger with respect to $t_{\rm c}$; then, the packet undergoes multiple reflections before being dissipated. In this case, the fraction of initial energy lost in transmitted waves becomes larger. This explains why $E_{\rm diss}$ decreases while $ E_{\rm t}^{\rm A} $ and $ E_{\rm t}^{\rm M} $ increase with increasing S. Clearly, both $t_{\rm d}$ and $t_{\rm c}$ are different for different packets and they depend on the form of the equilibrium magnetic field; moreover, the wavevector of the Alfvénic packet is modified by the reflection process. For these reasons, a quantitative evaluation of $E_{\rm diss}$, $ E_{\rm t}^{\rm A} $, $ E_{\rm t}^{\rm M} $, and $E_{\rm r}^{\rm M} $ has required the explicit solutions of both Eqs. (5)-(7) for all the packets, as well as the calculation of quantities involved in each reflection process. In particular, for $S\sim 10^{10}$ the dissipated energy is $E_{\rm diss}\sim 0.2{-}0.4$; then, a non negligible fraction of the initial wave energy is dissipated inside the equilibrium structure, even for extremely high values of the Reynolds number. This fact is due to the fast dissipation process (Malara et al. 2000), which is related to the three-dimensional geometry of the equilibrium magnetic field $\vec{B}_0$.

In Fig. 6 it is seen that the fraction $E_{\rm diss}$ of dissipated energy is larger in the case of a nonvanishing equilibrium current ( $\alpha = 7\pi $) with respect to the potential field case ($\alpha =0$). This could be due to two different effects: either the dissipation time $t_{\rm d}$ of single packets on average decreases with increasing $\alpha$, which would correspond to a larger wavevector growth rate; or the crossing time $t_{\rm c}$ increases with $\alpha$, which could be due to more twisted and then longer magnetic lines, and/or to a smaller average value of the Alfvén velocity. This point needs further investigation which is left for a future work.

In the lower panel of Fig. 6 the quantities $E_{\rm diss}/E_0$, $E_{\rm t}^{\rm A}/E_0$, $E_{\rm t}^{\rm M}/E_0$, and $E_{\rm r}^{\rm M}/E_0$ are shown for r=0.032. In this case both transmission coefficients are larger with respect to the case r=0.01 (Figs. 5); then, the energy lost in transmitted waves is larger, while the dissipated energy is smaller. However, the dependence of these quantities on the Reynolds number S and on the parameter $\alpha$ is qualitatively similar to that found for r=0.01. In particular, for large values of S, ( $S\sim 10^{10}$) the ratio $E_{\rm diss}/E_0\simeq 0.2$. The energy of reflected magnetosonic waves is negligible also for r=0.032.

4 Summary and conclusions

In this paper we have studied the evolution and dissipation of Alfvénic perturbations, which propagate inside a 3D inhomogeneous equilibrium magnetic field. This magnetic field belongs to a wide class of solutions of the force-free equilibrium equations (Nakagawa & Raadu 1972), obtained under some assumptions: namely, constant $\alpha$; periodicity in the horizontal directions; vanishing magnetic field for $z \rightarrow +\infty$. The magnetic field is a superposition of contributions (harmonics) with different amplitudes and at different spatial scales. We used this solution to represent a typical structure of a magnetic field above a quiet-Sun region, which is characterized by several regions of the two polarities randomly distributed in different locations and at various spatial scales. To model this kind of configuration we used a power-law spectrum for the harmonic amplitudes and random phases, while $\alpha$, which determines the current density, has been used as a free parameter of the model. This equilibrium magnetic field does not represent any particular structure deduced from observations: it is intended only as a model for a typical configuration of the coronal magnetic field. The underlying assumption is that our results on wave dissipation would not be modified if another equilibrium magnetic field of the same class would be used (for example, using a different choice for the phases of harmonics).

The propagation and dissipation of Alfvénic perturbations within the equilibrium structure has been described using an approach based on a WKB expansion of the cold plasma MHD equations (Malara et al. 2003). Some assumptions have been used in this approach: first, it has been assumed that the typical wavelength of fluctuations $\lambda$is smaller than the spatial scale L0 of the equilibrium structure. Since the wavevector of perturbations tends to increase in time, if this condition is verified at the initial time, then it will be satisfied for all the subsequent times. The initial condition $\lambda \ll L_0$ indicates that the wavelengths considered in our model do not represent the main contribution in the spectrum directly produced by photospheric motions. However, in consequence of the coupling with both transverse and longitudinal inhomogeneities at scale L0, waves at wavelengths smaller than L0 are formed. This has been suggested, for instance, for waves generated by the chromospheric network activity (Axford & McKenzie 1992). Moreover, there are indications that in closed magnetic structures, Alfvénic perturbations at wavelengths smaller than the loop length are generated by resonance phenomena (Nigro et al. 2004; Veltri et al. 2004); such wavelengths would represent the main contribution to the spectrum of Alfvénic fluctuations inside closed structures. The dissipation mechanism here studied can be applied to these waves.

The second assumption is that perturbations have a small amplitude with respect to the equilibrium structure: $\delta v/c_{\rm A} \ll 1$. It is commonly assumed that in corona the Alfvén velocity is $c_{\rm A} \sim 10^3$ km s-1 , while typical values of velocity fluctuation $\delta v \sim 20{-}40$ km s-1 can be deduced from measures of nonthermal broadening of coronal lines (Acton et al. 1981; Warren et al. 1997; Chae et al. 1998). Thus, the small amplitude assumption for fluctuations is verified in the corona. From a theoretical point of view, this assumption is commonly used to neglect nonlinear terms in the MHD equations which could give rise to a transfer of the energy of Alfvénic fluctuations to compressive modes. The efficiency of such transfer does not depend only on the amplitude but also on the wavelength of the Alfvén wave. Actually, Nakariakov et al. (1997) have shown that phase-mixing of Alfvén waves on a purely transverse inhomogeneity, and the consequent increase of magnetic pressure gradient, gives rise to magnetosonic waves, even for a small amplitude wave. A similar effect could be present also in the case considered here, since it is driven by the increase of the wavevector of Alfvén wave. However, provided that the Alfvén wave amplitude is not strongly nonlinear ( $\delta v/c_{\rm A} \lesssim 0.1$), the growth of the driven compressive fluctuation due to phase-mixing saturates at amplitudes much smaller than that of the Alfvén wave (Botha et al. 2000; Tsiklauri et al. 2001). This suggests that also in our case the coupling of Alfvén waves with compressive modes due to nonlinear terms can be neglected.

Concerning this point, Tsiklauri & Nakariakov (2002a) and Tsiklauri et al. (2003) have studied the time evolution of a 3D MHD pulse propagating in a transversely inhomogeneous medium, finding a coupling between the Alfvénic and the magnetosonic part of the perturbation; such coupling is present at the linear level. However, their results are not directly comparable with that of the present paper because their initial condition already contains a mix of a magnetosonic and an Alfvénic component. On the contrary, here we considered only purely Alfvénic perturbations, which are polarized perpendicularly both to $\vec{B}_0$ and to the local wavevector (Malara et al. 2003). Our configuration corresponds to the case $\alpha_y =0$ of Tsiklauri & Nakariakov (2002a) and Tsiklauri et al. (2003). In that case, these authors also find no energy transfer from the Alfvén wave to magnetosonic perturbation. More recently, Del Zanna et al. (2005) have studied the propagation of a purely Alfvénic pulse within an arcade. The equilibrium magnetic field they used can be obtained from our expression (17) for $\alpha =0$ and retaining the contribution of a single wavenumber $\kappa$. These authors find a coupling between Alfvén and magnetosonic waves. However, in that case the coupling is nonlinear, due to the large amplitude of the initial pulse, the velocity perturbation being comparable with the background Alfvén velocity.

Some general result which had been obtained for simpler background configurations (Malara et al. 2000) have been recovered also in the more complex equilibrium magnetic field considered here. In particular, the wavevector of Alfvénic packets tends to grow exponentially, giving origin to the scaling law (1) for the dissipation time $t_{\rm d}$, typical of fast dissipation. The dissipation time determines the typical distance over which a packet is dissipated; then, in open magnetic configurations the determination of $t_{\rm d}$ is crucial in order to establish whether an Alfvénic disturbance is dissipated before it leaves the corona. In a closed magnetic structure like that considered here the dissipation time plays a slightly different role, since Alfvénic perturbations are confined to propagate within the structure. However, since the corona is dynamically coupled to lower layers in the solar atmosphere, this confinement is not perfect and a certain leakage of fluctuating energy takes place at the coronal base (e.g., Ofman 2002). Then, the fraction of the initial perturbation energy which is dissipated inside the corona is determined by a balance between the dissipation time $t_{\rm d}$, which depends on the Reynolds number S, and the efficiency of the wave reflection process, which is quantified by the reflection coefficient $r_{\rm A}$. In our model we found that even for very large values of the Reynolds number $S\sim 10^{10}$ a non negligible fraction of the initial perturbation energy is dissipated within the coronal structure.

The value of the reflection coefficient $r_{\rm A}$ depends on how the transition between the corona and lower atmospheric layers is modelled. In this paper we represented such a transition by a sharp discontinuity in the background density. Other models have assumed an exponential dependence of the Alfvén velocity on the vertical coordinate z (e.g., De Pontieu et al. 2001; Ofman 2002) below the corona. Although this assumption is more realistic than a discontinuous density, in the other models both the background magnetic field and the perturbation wavevector have been assumed to be vertically directed, so that no magnetosonic perturbations are generated during the reflection process. On the contrary, in our model we tried to account for the geometry of the background magnetic field in calculating the wave reflection process; in particular, we took into account the local orientation of both $\vec{B}_0$ and of the wavevector of the incident wave. Moreover, the generation of both reflected and transmitted magnetosonic waves has been included in the model. We point out that in our model the detailed calculation of reflection and transmission coefficients of Alfvénic and magnetosonic perturbations has been possible as a consequence of the simplifying assumption of a discontinuous background density. However, DePontieu et al. (2001) found that, in the geometry they considered, the reflection coefficient obtained by using an exponential profile for the Alfvén velocity is lower than in the case of a discontinuous profile. This suggests that also in the more complex geometry considered here, a continuous profile of $c_{\rm A0}$ would give a reflection coefficient $r_{\rm A}$ lower than what we calculated, and the fraction of dissipated energy (e.g., $E_{\rm diss}/E_0 \sim 0.2{-}0.4$ at $S=2\times 10^{10}$) would be lowered, as well.

Acknowledgements
This work was partially supported by the MIUR (Ministero dell'Istruzione, dell'Universitá e della Ricerca) through a National Project Fund (Cofin 2002) and by the European Community within the Research Training Network "Turbulence in Space Plamas, Theory, Observation and Simulation''.

References

 

  
5 Online Material

Appendix A: Wave reflection at the discontinuity

We want to describe how Alfvén waves impinging over a planar density discontinuity give origin to both reflected and transmitted waves, and to calculate reflection and transmission coefficients. The discontinuity represents an extremely simplified model for the sharp density variation associated with the solar transition region. The width $\delta_{\rm t} < 10^3$ km of the transition region is much smaller than the typical variation scale $L_0 \sim 2.5 {-} 3 \times 10^4$ km of the background equilibrium structure. We consider a region $\Omega$ located across the density variation layer; the typical linear size of such a region is $l_\Omega$, with $\delta_{\rm t} \ll l_\Omega \ll l_0$. As a consequence of these inequalities, within the region $\Omega$ we will represent the background density variation as a real discontinuity, while on both sides of such a discontinuity the background equilibrium structure will be considered as uniform. The discontinuity is assumed to be planar and located at z=0. Moreover, we assume that the background magnetic field $\vec{B}_0$ is continuous through the discontinuity.

A.1 Linear modes

In the region $\Omega$, on both sides of the discontinuity we decompose physical quantities as the sum of a uniform background value plus a small amplitude perturbation:

 
                               $\displaystyle \rho (\vec{x},t) = \rho_0 + \varepsilon \rho_1 (\vec{x},t)$  
    $\displaystyle \vec{v} (\vec{x},t) = \varepsilon \vec{v}_1 (\vec{x},t)$ (A.1)
    $\displaystyle \vec{B} (\vec{x},t) = \vec{B}_0 + \varepsilon \vec{B}_1 (\vec{x},t)$  

where $\varepsilon \ll 1$ gives the amplitude of perturbations, and the indices "0'' and "1'' indicate quantities relative to the background structure and to the perturbation, respectively. The properties of perturbations can be derived following a standard procedure that we briefly outline. The expressions (A.1) are inserted into the Eqs. (2)-(4), which are expanded in powers of the parameter $\varepsilon$, retaining only terms at the lowest order $O(\varepsilon^1)$. A Fourier transform is carried out, according to

 \begin{displaymath}
f_1(\vec{x},t) = \int {{\rm d}\vec{k}{\rm d}\omega \over (2\pi)^4}
f_1 (\vec{k},\omega)
\exp (\vec{k}\cdot \vec{x}-\omega t)
\end{displaymath} (A.2)

where $ f(\vec{x},t)$ indicates $\rho_1$ or a component of $\vec{v}_1$ or of $\vec{B}_1$, and $f_1 (\vec{k},\omega)$ is the corresponding Fourier transform. For any given value $\vec{k}$ of the wavevector, nonvanishing solutions are found only when the frequency $\omega$ satisfies one of the following four dispersion relations:

 \begin{displaymath}
\omega = \omega^{\pm\rm A}(\vec{k}) = \sigma_{\rm A} \vec{c}_A \cdot \vec{k}
\end{displaymath} (A.3)
 \begin{displaymath}
\omega=\omega^{\pm \rm M}(\vec{k}) = \sigma_{\rm M} c_{\rm A} k
\end{displaymath} (A.4)

which give the frequency of Alfvén waves and of fast magnetosonic waves, respectively. $\sigma_{\rm A}=\pm 1$ and $\sigma_{\rm M}=\pm 1$ determine the propagation direction of Alfvén waves (parallel or antiparallel to $\vec{B}_0$) and of magnetosonic waves (parallel or antiparallel to $\vec{k}$), respectively.

The solution of the linearized equations is a sum of contributions, each due to a single mode:

 \begin{displaymath}
f_1(\vec{k},\omega ) = \sum_\alpha f_1^\alpha (\vec{k}) 2\pi
\delta (\omega - \omega^\alpha (\vec{k}))
\end{displaymath} (A.5)

where the upper Greek index identifies the contribution of a particular mode to the total perturbed field. Inserting the form (A.5) into Eq. (A.2) we obtain the equation

 \begin{displaymath}
f_1(\vec{x},t) = \sum_\alpha \int {{\rm d}\vec{k}\over (2 \p...
...c{k} \cdot \vec{x} - \omega^\alpha (\vec{k})t \right] \right\}
\end{displaymath} (A.6)

expressing the field $f_1(\vec{x},t)$ as a superposition of waves, each with different wavevector $\vec{k}$ and belonging to a different mode.


  \begin{figure}
\par\includegraphics[width=6.1cm,clip]{2272fi12.eps}
\end{figure} Figure A.1: The $(\xi , \eta , \zeta )$ reference frame is represented, along with the wavevector $\vec{k}$ and the magnetic field $\vec{B}_0$, whose direction is given by the unit vectors $\vec{e}_k$, $\vec{e}_B$, respectively. The unit vector $\vec{e}_{\rm A}$ gives the polarizarion of Alfvén waves, while the unit vectors $\vec{e}_1$ and $\vec{e}_2$ are in the $\xi \eta $ plane, perpendicular to $\vec{k}$ and to $\vec{B}_0$, respectively. The dashed line is parallel to $\vec{e}_2$.

In order to express the solution of the linearized system, for any given value of $\vec{k}$ it is useful to define a Cartesian reference frame $(\xi , \eta , \zeta )$ with the $\xi$-axis along the wavevector $\vec{k}$, and the $\xi \eta $-plane containing the direction of the background magnetic field $\vec{B}_0$ (see Fig. A.1). With respect to such reference frame, for the two magnetosonic modes we find:

 
                                          $\displaystyle \rho^{\pm\rm M}_1(\vec{k} )=\rho_0 \sin \theta a_{\pm \rm M}(\vec...
...M}_{1\xi}(\vec{k})= \sigma_{\rm M} c_{\rm A} \sin \theta
a_{\pm \rm M}(\vec{k})$  
    $\displaystyle v^{\pm \rm M}_{1\eta}(\vec{k})= -\sigma_{\rm M} c_{\rm A} \cos \t...
...M}(\vec{k}) , \;\;\;
B^{\pm \rm M}_{1\eta}(\vec{k})= B_0 a_{\pm \rm M}(\vec{k})$ (A.7)
    $\displaystyle v^{\pm \rm M}_{1\zeta}(\vec{k})=
B^{\pm \rm M}_{1\xi}(\vec{k})=
B^{\pm \rm M}_{1\zeta}(\vec{k})=0$  

where $\theta$ is the angle between $\vec{k}$ and $\vec{B}_0$, and $a_{\pm \rm M}(\vec{k})$ is the arbitrary (dimensionless) amplitude of the magnetosonic mode. Note that the velocity perturbation of magnetosonic waves is directed in the $\xi \eta $-plane, perpendicular to $\vec{B}_0$. In a similar way, for the two Alfvénic modes we find:
 
    $\displaystyle v^{\pm\rm A}_{1\zeta}(\vec{k})=- \sigma_{\rm A} c_{\rm A} a_{\pm\rm A}(\vec{k})
, \;\;\;
B^{\pm\rm A}_{1\zeta}(\vec{k})= B_0 a_{\pm\rm A}(\vec{k})$  
    $\displaystyle \rho^{\pm\rm A}_1(\vec{k}) =
v^{\pm\rm A}_{1\xi}(\vec{k})=
v^{\pm...
...{1\eta}(\vec{k})=
B^{\pm\rm A}_{1\xi}(\vec{k})=
B^{\pm\rm A}_{1\eta}(\vec{k})=0$ (A.8)

where $a_{\pm\rm A}(\vec{k})$ is the amplitude of the Alfvénic mode. Eqs. (A.7)-(A.8) indicate that the total velocity perturbation due to both Alfvénic and to magnetosonic perturbations lies in the plane perpendicular to $\vec{B}_0$. Actually, at the order $O(\varepsilon^1)$ the magnetic force $\vec{f}_{\rm M}=(1/4\pi)(\nabla \times \vec{B}_1)\times \vec{B}_0$ is perpendicular to $\vec{B}_0$, while the pressure gradient vanishes in a cold plasma. Thus, no velocity perturbations can be generated parallel to $\vec{B}_0$.

The above expressions of the perturbed fields can be written in a more general form, independent of the particular reference frame we used. For this purpose, we introduce the following unit vectors, which are represented in Fig. A.1: $\vec{e}_k(\vec{k})=\vec{k}/k$ and $\vec{e}_B=\vec{B}_0/B_0$, directed along $\vec{k}$ and $\vec{B}_0$, respectively; $\vec{e}_{\rm A}(\vec{k})=(\vec{e}_k(\vec{k}) \times \vec{e}_B)/
\vert\vec{e}_k(\vec{k}) \times \vec{e}_B\vert$ directed parallel to the Alfvénic perturbations; $\vec{e}_1(\vec{k})=\vec{e}_{\rm A}(\vec{k}) \times
\vec{e}_k(\vec{k})$ and $\vec{e}_2(\vec{k})=\vec{e}_B \times \vec{e}_{\rm A}(\vec{k})$ in the direction of magnetosonic magnetic field and velocity perturbation, respectively. Using such unit vectors, the perturbations relative to magnetosonic and Alfvénic modes can be written in the form:

 \begin{displaymath}
\rho^{\pm \rm M}_1(\vec{k})=\rho_0
\left [ \vec{e}_2(\vec{k...
...m M}_1(\vec{k})=B_0
a_{\pm \rm M}(\vec{k}) \vec{e}_1(\vec{k})
\end{displaymath} (A.9)

and

 \begin{displaymath}
\rho^{\pm\rm A}_1(\vec{k})=0 , \;\;\;
\vec{v}^{\pm\rm A}_1(...
...ec{k},)= B_0
a_{\pm\rm A}(\vec{k}) \vec{e}_{\rm A}(\vec{k}).
\end{displaymath} (A.10)

The unit vectors $\vec{e}_{\rm A}(\vec{k})$, $\vec{e}_1(\vec{k})$ and $\vec{e}_2(\vec{k})$ are not defined when $\vec{k}$ is parallel to $\vec{B}_0$. In this particular case the magnetosonic modes have the same properties as the Alfvén modes: (i) $\omega^{\pm \rm M}=\omega^{\pm\rm A}$; (ii) the velocity and magnetic field perturbations of magnetosonic mode are parallel and they are both directed perpendicular to $\vec{B}_0$ and to Alfvénic perturbations (Eqs. (A.4) and (A.7) for $\theta =0$). Thus, for $\vec{k}$ parallel to $\vec{B}_0$ we can still use the expressions (A.9) and (A.10), where $\vec{e}_1(\vec{k})=\vec{e}_2(\vec{k})$ and $\vec{e}_{\rm A}(\vec{k})$ are two arbitrary unit vectors perpendicular to each other and to $\vec{B}_0$.


  \begin{figure}
\par\includegraphics[width=6.5cm,clip]{2272fi13.eps}
\end{figure} Figure A.2: The (X, Y, Z ) reference frame is represented. The XY-plane corresponds to the base of the corona; the background magnetic field is in the YZ plane, at an angle $\phi $ with the horizontal direction.

Within the region $\Omega$ we introduce another Cartesian reference frame (X,Y,Z), represented in Fig. A.2. This reference frame is obtained by rotating the reference frame xyz around the z-direction, so that the background magnetic field $\vec{B}_0$ lies in the YZ-plane, forming an angle $\phi $ with the Y-axis. Thus, the Z-axis is directed parallel to z, the density discontinuity lies in the XY-plane and the corona lies in the half-space Z>0. In this reference frame the wavevectors $\vec{k}$ of the different modes can have an arbitrary orientation. The XYZ frame is defined within the small region $\Omega$, where the background magnetic field can be considered uniform; however, in different points of the coronal base z=0 the orientation of X and Y axes is different, according to local the direction of the background magnetic field $\vec{B}_0$. In the (X,Y,Z) reference frame the above-defined unit vectors have the following expressions:

 \begin{displaymath}
\vec{e}_{\rm A} = {\left( k_Y \sin \phi - k_Z \cos \phi \rig...
..._Y \sin \phi - k_Z \cos \phi \right)^2 + k_X^2
\right]^{1/2}}
\end{displaymath} (A.11)


 \begin{displaymath}
\vec{e}_1 = { - k_X \left[ k_Y \cos \phi + k_Z \sin \phi \ri...
...Y \sin \phi - k_Z \cos \phi \right)^2 + k_X^2
\right]^{1/2} }
\end{displaymath} (A.12)


 \begin{displaymath}
\vec{e}_2 = { k_X \vec{e}_X +
\left( k_Y \sin^2 \phi - k_Z ...
...k_Y \sin \phi - k_Z \cos \phi \right)^2 + k_X^2 \right]^{1/2}}
\end{displaymath} (A.13)

where $\vec{e}_X$, $\vec{e}_Y$ and $\vec{e}_Z$ are the unit vectors of the X, Y and Z axes, respectively. Using these expressions and the Eqs. (A.7), the perturbations of the magnetosonic mode in the XYZ reference frame are written as:

 \begin{displaymath}
\rho^{\pm \rm M}_1(\vec{k})={\rho_0 \over k^2}
\left[ \left(...
...Z \cos \phi \right)^2 + k_X^2
\right] a'_{\pm \rm M}(\vec{k})
\end{displaymath} (A.14)


 \begin{displaymath}
\vec{v}^{\pm \rm M}_1(\vec{k})= {\sigma_{\rm M} c_{\rm A} \o...
...s \phi k_Y \right) \vec{e}_Z \right ]
a'_{\pm \rm M}(\vec{k})
\end{displaymath} (A.15)


 \begin{displaymath}
\vec{B}^{\pm \rm M}_1(\vec{k})= {B_0 \over k^2}
\left \{ - k...
...cos \phi \right ]
\vec{e}_Z \right \}
a'_{\pm \rm M}(\vec{k})
\end{displaymath} (A.16)

and the perturbations of the Alfvénic mode are:

 \begin{displaymath}
\rho^{\pm\rm A}_1(\vec{k})=0
\end{displaymath} (A.17)


 \begin{displaymath}
\vec{v}^{\pm\rm A}_1(\vec{k})={- \sigma_{\rm A} c_{\rm A} \o...
...e}_Y + k_X \cos \phi \vec{e}_Z \right ]
a'_{\pm\rm A}(\vec{k})
\end{displaymath} (A.18)


 \begin{displaymath}
\vec{B}^{\pm\rm A}_1(\vec{k})= {B_0 \over k}
\left [ \left( ...
...}_Y + k_X \cos \phi \vec{e}_Z \right ]
a'_{\pm\rm A}(\vec{k}).
\end{displaymath} (A.19)

To simplify the expressions (A.14)-(A.19) we have re-defined the dimensionless wave amplitude as

 \begin{displaymath}
a'_{\rm\pm M,A}(\vec{k}) =
{k \over \left [ \left( k_Y \sin...
...i \right) ^2 + k_X^2
\right ] ^{1/2}} a_{\rm\pm M,A}(\vec{k}).
\end{displaymath} (A.20)

The above expressions cannot be used for $\vec{k}$ parallel to $\vec{B}_0$. In fact, in such a case we have kX=0, $k_Y=k \cos \phi$, and $k_Z=k \sin \phi$; thus, the denominator of Eq. (A.20) would vanish.

In the particular case when $\vec{k}$ is parallel to $\vec{B}_0$, magnetosonic and Alfvén waves cannot be distinguished. So, we conventionally choose:

 \begin{displaymath}
\vec{e}_1(\vec{k})=-\vec{e}_2(\vec{k})=- \sin \phi \vec{e}_Y...
...\phi \vec{e}_Z , \;\;\;\;
\vec{e}_{\rm A}(\vec{k})=\vec{e}_X.
\end{displaymath} (A.21)

Thus, for $\vec{k}$ parallel to $\vec{B}_0$ the perturbations of the magnetosonic mode are given by

 \begin{displaymath}
\rho^{\pm \rm M}_1(\vec{k})=0
\end{displaymath} (A.22)


 \begin{displaymath}
\vec{v}^{\pm \rm M}_1(\vec{k})= -\sigma_{\rm M} c_{\rm A}
\...
...ec{e}_Y + \cos \phi \vec{e}_Z \right)
a'_{\pm \rm M}(\vec{k})
\end{displaymath} (A.23)


 \begin{displaymath}
\vec{B}^{\pm \rm M}_1(\vec{k})= B_0
\left( - \sin \phi \vec{e}_Y + \cos \phi \vec{e}_Z \right)
a'_{\pm \rm M}(\vec{k})
\end{displaymath} (A.24)

while the perturbations of the Alfvénic mode are:

 \begin{displaymath}
\rho^{\pm\rm A}_1(\vec{k})=0
\end{displaymath} (A.25)


 \begin{displaymath}
\vec{v}^{\pm\rm A}_1(\vec{k})= - \sigma_{\rm M} c_{\rm A}
a'_{\pm\rm A}(\vec{k}) \vec{e}_X
\end{displaymath} (A.26)


 \begin{displaymath}
\vec{B}^{\pm\rm A}_1(\vec{k})= B_0
a'_{\pm\rm A}(\vec{k}) \vec{e}_X.
\end{displaymath} (A.27)

A.2 Continuity conditions

The discontinuity in the equilibrium density is located at the surface Z=0. The values of the background density for Z>0 and Z<0 are indicated by $\rho_0^+$ and $\rho_0^-$, respectively. Some physical quantity must be continuous across the discontinuity. In order to find these continuity conditions, it is useful to re-write the cold plasma MHD Eqs. (2)-(4) in the following conservative form:

 \begin{displaymath}
{\partial \rho \over \partial t} + {\partial \over \partial X_j}
\left( \rho v_j \right) = 0
\end{displaymath} (A.28)


 \begin{displaymath}
{\partial \over \partial t} \left( \rho v_i \right) +
{\part...
...k \over 8 \pi} \delta_{ij} -
{B_i B_j \over 4 \pi} \right) = 0
\end{displaymath} (A.29)


 \begin{displaymath}
{\partial B_i \over \partial t} + {\partial \over \partial X_j}
\left( B_i v_j - B_j v_i \right) =0.
\end{displaymath} (A.30)

Each component of Eqs. (A.28)-(A.30) can be written in the form:

 \begin{displaymath}
{\partial f \over \partial t} + {\partial \phi_j \over \partial X_j}=0
\end{displaymath} (A.31)

where f represents either the mass density or a component of the momentum density or of the magnetic field, and ${\bf\phi}$ is the corresponding flux density. We consider a cylindrical volume V located across the discontinuity surface Z=0, with height $\delta$. The bases of the cylinder, which are parallel to the discontinuity surface, have area S and perimeter p. We integrate the Eq. (A.31) over the volume V, transforming the volume integral of the second term into a surface integral, by means of the divergence theorem. Dividing by S we find:

 \begin{displaymath}
{{\rm d} \langle f \rangle_V \over {\rm d}t} \delta + \langl...
...hi_Z^- \rangle_S + \langle \phi \rangle_L {p\delta \over S}=0.
\end{displaymath} (A.32)

where $\langle f \rangle_V$ is the average of f over the volume V; $\langle \phi_Z^\pm \rangle_S$ is the the Z-component of the flux ${\bf\phi}$ averaged over the bases S+ and S-, above and below the surface Z=0, respectively; and $\langle \phi \rangle_L$ is average of the flux over the lateral surface of the cylinder. In the limit of vanishing $\delta$, assuming that all the above averaged quantities remain finite, we obtain:

\begin{displaymath}\langle \phi_Z^+ \rangle_S - \langle \phi_Z^- \rangle_S =
\int_S \left( \phi_Z^+ - \phi_Z^- \right) {\rm d}s = 0.
\end{displaymath}

Since the basis S of the cylinder is arbitrary, the above equation implies

 \begin{displaymath}
\phi_Z^+ = \phi_Z^-
\end{displaymath} (A.33)

indicating that the flux density component $\phi_Z$ perpendicular to the discontinuity is continuous across the discontinuity.

We apply the above argument both to the momentum Eq. (A.29) and to the magnetic field Eq. (A.30), and we will examine the mass density equation later. The continuity of the momentum flux gives the equations:

 \begin{displaymath}
\left [ \rho v_i v_Z - {B_i B_Z \over 4 \pi} \right ]^+ =
\l...
...[ \rho v_i v_Z - {B_i B_Z \over 4 \pi} \right ]^- , \;\; i=X,Y
\end{displaymath} (A.34)


 \begin{displaymath}
\left [ \rho v_Z^2 + {B_X^2 + B_Y^2 - B_Z^2 \over 8\pi}\righ...
...\rho v_Z^2 + {B_X^2 + B_Y^2 - B_Z^2 \over 8\pi}\right ]^-\cdot
\end{displaymath} (A.35)

The continuity of the magnetic field flux gives the equations:

 \begin{displaymath}
\left [ B_i v_Z - B_Z v_i \right ]^+ = \left [ B_i v_Z - B_Z v_i \right ]^-
, \;\; i=X,Y
\end{displaymath} (A.36)

while the third equation, corresponding to i=Z, is identically satisfied.

We decompose the quantities as a sum of a background equilibrium plus a small amplitude perturbation, as in Eqs. (A.1). Taking into account that the background magnetic field has been assumed continuous across the discontinuity ( $\vec{B}_0^+ = \vec{B}_0^-$), at the first order in the perturbation amplitude $\varepsilon$ Eqs. (A.34) and (A.36) give:

 
B1X+ = B1X- (A.37)


 \begin{displaymath}
B_{0Z} \left( B_{1Y}^+ - B_{1Y}^- \right) +
B_{0Y} \left( B_{1Z}^+ - B_{1Z}^- \right) = 0
\end{displaymath} (A.38)


 \begin{displaymath}
B_{0Y} \left( B_{1Y}^+ - B_{1Y}^- \right) -
B_{0Z} \left( B_{1Z}^+ - B_{1Z}^- \right) = 0
\end{displaymath} (A.39)


 
v1X+ = v1X- (A.40)


 \begin{displaymath}
B_{0Y} \left( v_{1Z}^+ - v_{1Z}^- \right) -
B_{0Z} \left( v_{1Y}^+ - v_{1Y}^- \right) = 0.
\end{displaymath} (A.41)

Being $B_0 \ne 0$, Eqs. (A.38) and (A.39) imply

 
B1Y+ = B1Y- (A.42)

and

 
B1Z+ = B1Z-. (A.43)

This result, along with Eq. (A.37), indicates that the perturbation magnetic field $\vec{B}_1$ is continuous across the discontinuity. The Eq. (A.41) can be written in the form

 \begin{displaymath}
v_{1Z}^+ \cos \phi - v_{1Y}^+ \sin \phi =
v_{1Z}^- \cos \phi - v_{1Y}^- \sin \phi
\end{displaymath} (A.44)

which is equivalent to

\begin{displaymath}\vec{v}_1^+ \cdot \vec{e}_\perp = \vec{v}_1^- \cdot \vec{e}_\perp
\end{displaymath}

where $\vec{e}_\perp = -\sin \phi \vec{e}_Y + \cos \phi \vec{e}_Z$ is the unit vector lying in the YZ-plane and perpendicular to $\vec{B}_0$(Fig. A.2). Since the velocity perturbation $\vec{v}_1$ is in the plane perpendicular to $\vec{B}_0$, the conditions (A.40) and (A.44) indicate that $\vec{v}_1$ is also continuous across the discontinuity.

Concerning the condition (A.43), it expresses the fact that $\vec{B}_1$ is divergenceless across the discontinuity. However, the Eq. (A.43) does not need to be explicitly used, since it can be obtained as a direct consequence of the other continuity conditions. To show this fact, we consider the Z-component of the induction Eq. (4), which is linearized on both sides of the discontinuity surface Z=0. Subtracting the two resulting equations we find

 \begin{displaymath}
{\partial B_{1Z} \over \partial t}^+ - {\partial B_{1Z} \ove...
...tial Y}^+ -
{\partial v_{1Y} \over \partial Y}^- \right)\cdot
\end{displaymath} (A.45)

Since all the components of the velocity perturbation $\vec{v}_1$ are continuous across the surface Z=0 for any value of X and Y, this implies that all the space derivatives in the RHS of Eq. (A.45) are continuous, as well. Then, the RHS of Eq. (A.45) vanishes, giving

\begin{displaymath}{\partial B_{1Z} \over \partial t}^+ = {\partial B_{1Z} \over \partial t}^-
\end{displaymath}

which implies

B1Z+(X,Y,t) = B1Z-(X,Y,t) + c(X,Y)

where c(X,Y) depends on the initial condition. If at the initial time B1Z+ = B1Z-, then c(X,Y)=0, and the Eq. (A.43) follows. Thus, it is not necessary to impose the B1Z continuity condition (A.43) , because it will be automatically satisfied once the continuity conditions on the perturbation velocity will be imposed.

Finally, we consider the mass flux $\rho_0 v_{1Z}$ across the discontinuity surface Z=0. Since the equilibrium density $\rho_0$ has different values on the two sides, while v1Z is continuous, it is clear that the mass flux is necessarily discontinuous: $(\rho_0 v_{1Z})^+ \ne (\rho_0 v_{1Z})^-$. Differences in the mass flux on the two sides of the discontinuity generate density fluctuations localized within the width $\delta_{\rm t}$ of the discontinuity. These density fluctuations are simply due to the undulations of the discontinuity surface, generated by the component v1Z of the velocity perturbation. The amplitude of $\delta \rho$ of such fluctuations can be estimated by Eq. (A.32), finding $\delta \rho \sim \vert \rho_0^+ - \rho_0^-\vert v_{1Z} \tau /\delta_{\rm t}$, where $\tau$ is the wave period of the velocity perturbation. However, since in a cold plasma the gradient pressure term has been dropped, these driven density variations have no effects on the dynamics of the plasma.

In summary, the continuity conditions which will be imposed at the discontinuity surface are the following:

 
B1X+(X,Y,Z=0,t) = B1X-(X,Y,Z=0,t) (A.46)


 
B1Y+(X,Y,Z=0,t) = B1Y-(X,Y,Z=0,t) (A.47)


 
v1X+(X,Y,Z=0,t) = v1X-(X,Y,Z=0,t) (A.48)


 \begin{displaymath}
v_{1Z}^+(X,Y,Z=0,t) \cos \phi - v_{1Y}^+(X,Y,Z=0,t) \sin \ph...
...v_{1Z}^-(X,Y,Z=0,t) \cos \phi - v_{1Y}^-(X,Y,Z=0,t) \sin \phi.
\end{displaymath} (A.49)

A.3 Reflected and transmitted waves

Within volumes smaller than the variation scale l0 of the background structure, the perturbations considered by the WKB theory can be approximated as a superposition of planar waves, each one with a given wavevector (e.g., Malara et al. 2003). Thus, within the region $\Omega$ we consider an Alfvén wave with a given arbitrary wavevector, incident on the discontinuity from the coronal side Z>0. At the discontinuity, four waves are generated: an Alfvén wave and a fast magnetosonic wave, both reflected back into the corona (Z>0); an Alfvén wave and a fast magnetosonic wave, both transmitted to Z<0. We will indicate quantities relative to these waves by the upper indices: "IA'' for the incident Alfvén wave; "RA'' and "TA'' for the reflected and transmitted Alfvén waves; "RM'' and "TM'' for the reflected and transmitted magnetosonic waves. Each of the five waves involved in this process have a given wavevector $\vec{k}^\beta$ and a corresponding frequency $\omega^\beta$ determined by the dispersion relations (A.3) and (A.4), where the index $\beta = \rm IA, RA,TA,RM,TM$. The velocity and magnetic field perturbations associated with such waves have the form given by Eq. (A.6):

 \begin{displaymath}
\vec{v}_1^\beta (\vec{X},t) = \vec{v}_1^\beta
\exp \left [ ...
...t(\vec{k}^\beta \cdot \vec{X} - \omega^\beta t\right) \right ]
\end{displaymath} (A.50)

The total perturbation above and below the discontinuity is a superposition of the perturbed fields associated with the single waves. Then, the continuity condition (A.46) takes the form
 
$\displaystyle B_{1X}^{\rm IA} \exp \left [ {\rm i} \left(k_X^{\rm IA}X + k_Y^{\...
...rm i} \left(k_X^{\rm RM}X + k_Y^{\rm RM}Y - \omega^{\rm RM} t\right)
\right ] =$      
$\displaystyle B_{1X}^{\rm TA} \exp \left [ {\rm i} \left(k_X^{\rm TA}X + k_Y^{\...
...rm i} \left(k_X^{\rm TM}X + k_Y^{\rm TM}Y - \omega^{\rm TM} t\right) \right ] .$     (A.51)

This equation must hold for any value of X, Y and t, which implies that the X and Y components of all wavevectors, as well as the frequencies, must be equal:

 \begin{displaymath}
k_X^{\rm RA} = k_X^{\rm RM} = k_X^{\rm TA} = k_X^{\rm TM} = k_X^{\rm IA} \equiv k_X
\end{displaymath} (A.52)


 \begin{displaymath}
k_Y^{\rm RA} = k_Y^{\rm RM} = k_Y^{\rm TA} = k_Y^{\rm TM} = k_Y^{\rm IA} \equiv k_Y
\end{displaymath} (A.53)


 \begin{displaymath}
\omega^{\rm RA} = \omega^{\rm RM} = \omega^{\rm TA} = \omega^{\rm TM} =
\omega^{\rm IA}
\end{displaymath} (A.54)

while the components $k_Z^\alpha$ of the five wavevectors are not necessarily equal. The same results (A.52)-(A.54) are obtained considering the other continuity conditions (A.47)-(A.49). The components $k_Z^\alpha$ of the wavevectors of reflected and transmitted waves can be calculated in terms of the wavevector $\vec{k}^{\rm IA}$ of the incident wave using the dispersion relations (A.3) and (A.4) and the Eqs. (A.52)-(A.54). Let us consider the single waves:

A.3.1 Reflected Alfvén wave
The equation $\omega^{\rm RA}=\omega^{\rm IA}$, along with the dispersion relation (A.3) of Alfvén waves, gives the equation

 \begin{displaymath}
\sigma_{\rm IA}\left(k_Y B_{0Y} + k_Z^{\rm IA} B_{0Z}\right) =
\sigma_{\rm RA}\left(k_Y B_{0Y} + k_Z^{\rm RA} B_{0Z}\right)
\end{displaymath} (A.55)

where $\sigma_{\rm IA}=\pm 1$ and $\sigma_{\rm RA}=\pm 1$ determine the propagation direction of the incident and of the reflected Alfvén wave, with respect to the direction of the equilibrium magnetic field. Since the incident wave propagates downward, while the reflected wave propagates upward, there are two possible situations: (i) if B0Z>0 the incident and the reflected Alfvén waves propagate antiparallel and parallel to $\vec{B}_0$, respectively; thus $\sigma_{\rm IA}=-1$ and $\sigma_{\rm RA}=1$; (ii) if B0Z<0, then $\sigma_{\rm IA}=1$ and $\sigma_{\rm RA}=-1$. The condition B0Z=0 can be verified only along lines or isolated points in the Z=0 plane, and it will not be considered. In summary, we have

 \begin{displaymath}
\sigma_{\rm IA}=-\sigma_B ; \;\;\; \sigma_{\rm RA}=\sigma_B
\end{displaymath} (A.56)

where $\sigma_B= {\rm sgn} (B_{0Z})=B_{0Z}/\vert B_{0Z}\vert$ is the sign of B0Z. Inserting this relations into the Eq. (A.55) we find the third component of the wavevector of the reflected Alfvén wave

 \begin{displaymath}
k_Z^{\rm RA}= - k_Z^{\rm IA}- 2{k_Y \over \tan \phi}\cdot
\end{displaymath} (A.57)

A.3.2 Transmitted Alfvén wave
In a similar way, the equation $\omega^{\rm TA}=\omega^{\rm IA}$and the dispersion relation (A.3) give the equation

 \begin{displaymath}
{\sigma_{\rm IA}\left(k_Y B_{0Y} + k_Z^{\rm IA} B_{0Z}\right...
...Y} + k_Z^{\rm TA} B_{0Z}\right) \over ({\rho_0^-})^{1/2}}\cdot
\end{displaymath} (A.58)

In this case the incident and the transmitted Alfvén waves both propagate downward. Then, the relation between the sign of B0Z and the propagation directions gives

 \begin{displaymath}
\sigma_{\rm IA}=\sigma_{\rm TA}=-\sigma_B .
\end{displaymath} (A.59)

Inserting this relation into the Eq. (A.58) we find

 \begin{displaymath}
k_Z^{\rm TA}={k_Z^{\rm IA} \over r} - \left( 1-{1 \over r} \right)
{k_Y \over \tan \phi}
\end{displaymath} (A.60)

where

 \begin{displaymath}
r={c_{\rm A}^- \over c_{\rm A}^+} = \left( {\rho_0^+ \over \rho_0^-} \right)^{1/2}
\end{displaymath} (A.61)

is the ratio between the Alfvén velocities on the two sides of the discontinuity. In the considered configuration, we have $r \ll 1$.

A.3.3 Reflected magnetosonic wave
The equation $\omega^{\rm IA}=\omega^{\rm RM}$ and the dispersion relations (A.3) and (A.4) give the equation

 \begin{displaymath}
\sigma_{\rm IA} \left( k_Y B_{0Y} + k_Z B_{0Z} \right) =
\s...
...RM} B_0 \left( k_X^2 + k_Y^2 + {k_Z^{\rm RM}}^2 \right)^{1/2}.
\end{displaymath} (A.62)

Since both B0 and $k^{\rm RM}$ are positive quantities, from this equation it follows that

 \begin{displaymath}
\sigma_{\rm RM}= \sigma_{\rm IA} {\rm sgn} \left(\vec{k}^{\rm IA} \cdot \vec{B}_0 \right).
\end{displaymath} (A.63)

Moreover, from Eq. (A.62) we derive:

 \begin{displaymath}
{k_Z^{\rm RM}}^2 = \left( {k_Z^{\rm IA}}^2 - k_Y^2 \right) \sin^2 \phi +
2 k_Y k_Z^{\rm IA} \sin \phi \cos \phi - k_X^2.
\end{displaymath} (A.64)

This equation indicates that $k_Z^{\rm RM}$ can be either real or imaginary, according to the sign of the RHS. In the former case the reflected magnetosonic wave is propagating, while in the latter case it is evanescent along the Z direction.

Let us suppose first that the reflected magnetosonic wave is propagating (real $k_Z^{\rm RM}$). In this case, the sign of $k_Z^{\rm RM}$ can be determined in the following way. The propagation direction of the reflected magnetosonic wave (parallel or antiparallel to $\vec{k}^{\rm RM}$) is determined by $\sigma_{\rm RM}=\pm 1$. Since the reflected magnetosonic wave propagates upward, $\sigma_{\rm RM}$ depends on the orientation of $\vec{k}^{\rm RM}$. In particular,

 \begin{displaymath}
\sigma_{\rm RM}= {\rm sgn} \left( k_Z^{\rm RM} \right).
\end{displaymath} (A.65)

This condition and the Eq. (A.63) determine the sign of $k_Z^{\rm RM}$ as a function of the incident wave:

 \begin{displaymath}
{\rm sgn} \left( k_Z^{\rm RM} \right) =
\sigma_{\rm IA} {\rm sgn} \left(\vec{k}^{\rm IA} \cdot \vec{B}_0 \right).
\end{displaymath} (A.66)

Then, from Eq. (A.64) the following expression for the Z-component of the wavevector of propagating reflected magnetosonic wave is obtained:

 \begin{displaymath}
k_Z^{\rm RM}=
\sigma_{\rm IA} {\rm sgn} \left(\vec{k}^{\rm ...
...2 k_Y k_Z^{\rm IA} \sin \phi \cos \phi - k_X^2 \right ]^{1/2}.
\end{displaymath} (A.67)

If the sign of the RHS of Eq. (A.64) is negative, then the reflected magnetosonic wave is evanescent (imaginary $k_Z^{\rm RM}$). In this case we write $ k_Z^{\rm RM}= {\rm i} \gamma^{\rm RM}$, and from Eq. (A.64) we find

 \begin{displaymath}
{\gamma^{\rm RM}}^2 =
\left( {k_Y^2 - k_Z^{\rm IA}}^2 \right) \sin^2 \phi -
2 k_Y k_Z^{\rm IA} \sin \phi \cos \phi + k_X^2 .
\end{displaymath} (A.68)

Since the wave amplitude must not diverge for $Z \rightarrow + \infty$, then $\gamma^{\rm RM}$ must be positive:

 \begin{displaymath}
\gamma^{\rm RM}=
\left [ \left( {k_Y^2 - k_Z^{\rm IA}}^2 \ri...
...2 k_Y k_Z^{\rm IA} \sin \phi \cos \phi + k_X^2 \right ]^{1/2}.
\end{displaymath} (A.69)

It is also useful to obtain a simple expression for the modulus of the wavevector $k^{\rm RM}= k_X^2 + k_Y^2 + k_Z^{\rm RM}$ of the reflected magnetosonic wave. From Eq. (A.62), we get

\begin{displaymath}k^{\rm RM} = \sigma_{\rm RM} \sigma_{\rm IA} {\vec{k}^{\rm IA} \cdot \vec{B}_0 \over B_0}\cdot
\end{displaymath}

Taking the absolute value of this equation, and using the fact that $ k^{\rm RM} \ge 0$, we obtain the expression for $k^{\rm RM}$:

 \begin{displaymath}
k^{\rm RM}={\vert\vec{k}^{\rm IA} \cdot \vec{B}_0\vert \over...
...}=
\big\vert k_Y \cos \phi + k_Z^{\rm IA} \sin \phi \big\vert.
\end{displaymath} (A.70)

Note that the Eq. (A.70) is valid both for non-evanescent and for evanescent magnetosonic wave, since in both cases $k^{\rm RM}$ is a real positive quantity.

A.3.4 Transmitted magnetosonic wave
The condition $\omega^{\rm IA}=\omega^{\rm TM}$ and the dispersion relations (A.3) and (A.4) give the equation

 \begin{displaymath}
\sigma_{\rm IA} {k_Y B_{0Y} + k_Z B_{0Z} \over {\rho_0^+}^{1...
... + {k_Z^{\rm TM}}^2 \right)^{1/2} \over {\rho_0^-}^{1/2}}\cdot
\end{displaymath} (A.71)

From this equation, it follows that

 \begin{displaymath}
\sigma_{\rm TM}= \sigma_{\rm IA} {\rm sgn} \left(\vec{k}^{\rm IA} \cdot \vec{B}_0 \right).
\end{displaymath} (A.72)

Moreover, from Eq. (A.71) we derive

 \begin{displaymath}
{k_Z^{\rm TM}}^2 = k_Y^2 \left( {\cos^2 \phi \over r^2} -1 \...
...+
{2 k_Y k_Z^{\rm IA} \sin \phi \cos \phi \over r^2} - k_X^2.
\end{displaymath} (A.73)

Also in this case, the transmitted magnetosonic wave can be either propagating or evanescent along Z, according to the sign of the RHS of Eq. (A.73).

Considering first a propagating transmitted wave, the propagation direction (parallel or antiparallel to $\vec{k}^{\rm TM}$) is determined by  $\sigma_{\rm TM}$. Since the transmitted magnetosonic wave propagates downward, the sign $\sigma_{\rm TM}$ is related to the orientation of $\vec{k}^{\rm TM}$ by

\begin{displaymath}\sigma_{\rm TM} = - {\rm sgn} \left( k_Z^{\rm TM} \right).
\end{displaymath}

This condition and the Eq. (A.72) determine the sign of $k_Z^{\rm TM}$:

 \begin{displaymath}
{\rm sgn} \left( k_Z^{\rm TM} \right) =
- \sigma_{\rm IA} {\rm sgn} \left(\vec{k}^{\rm IA} \cdot \vec{B}_0 \right).
\end{displaymath} (A.74)

Then, the expression for the z component of the wavevector of the transmitted magnetosonic wave, as a function of the incident wave, follows from Eq. (A.73):

 \begin{displaymath}
k_Z^{\rm TM} = - \sigma_{\rm IA} {\rm sgn} \left(\vec{k}^{\r...
...\rm IA} \sin \phi \cos \phi \over r^2} - k_X^2 \right ]^{1/2}.
\end{displaymath} (A.75)

In the case of evanescent transmitted magnetosonic wave (negative RHS of Eq. (A.73)), we write $k_Z^{\rm TM}= {\rm i} \gamma^{\rm TM}$, and from Eq. (A.73) we find

 \begin{displaymath}
{\gamma^{\rm TM}}^2 = k_Y^2 \left( 1 - {\cos^2 \phi \over r^...
...-
{2 k_Y k_Z^{\rm IA} \sin \phi \cos \phi \over r^2} + k_X^2.
\end{displaymath} (A.76)

Since the wave amplitude must not diverge for $Z \rightarrow - \infty$, then $\gamma^{\rm TM}$ must be negative. Thus

 \begin{displaymath}
\gamma^{\rm TM}=
- \left [ k_Y^2 \left( 1 - {\cos^2 \phi \ov...
...\rm IA} \sin \phi \cos \phi \over r^2} + k_X^2 \right ]^{1/2}.
\end{displaymath} (A.77)

In a similar way as for the reflected wave, a simple expression for the modulus of the wavevector $k^{\rm TM}= k_X^2 + k_Y^2 + k_Z^{\rm TM}$ of the transmitted magnetosonic wave can be obtained. From Eq. (A.71) we get

\begin{displaymath}k^{\rm TM}= \sigma_{\rm TM} \sigma_{\rm IA} {\vec{k}^{\rm IA}...
...over B_0}
\left( {\rho_0^- \over \rho_0^+} \right) ^{1/2}\cdot
\end{displaymath}

Taking the absolute value of the above equation and using the fact that $k^{\rm TM} \ge 0$ , we obtain

 \begin{displaymath}
k^{\rm TM}={1 \over r} {\vert\vec{k}^{\rm IA} \cdot \vec{B}_...
...r} \big\vert k_Y \cos \phi + k_Z^{\rm IA} \sin \phi \big\vert.
\end{displaymath} (A.78)

This equation is valid both for non-evanescent and for evanescent magnetosonic wave, since in both cases $k^{\rm TM}$ is a real positive quantity.

A.4 Equations for the wave amplitudes

The amplitudes of the waves involved in the reflection process can be determined by the continuity conditions at the surface Z=0. The continuity of B1X is expressed by Eq. (A.51); using the conditions (A.52)-(A.54) we obtain the equation

\begin{displaymath}B_{1X}^{\rm RA} + B_{1X}^{\rm RM} - B_{1X}^{\rm TA} - B_{1X}^{\rm TM}=-
B_{1X}^{\rm IA}.
\end{displaymath}

Using the expressions (A.16) and (A.19) for the perturbation magnetic field, we obtain:
 
             $\displaystyle \left( {a^{\rm RA}\over a^{\rm IA}} \right)
{k_Y \sin \phi - k_Z^{\rm RA} \cos \phi \over k^{\rm RA}} +$   $\displaystyle \left( {a^{\rm RM}\over a^{\rm IA}}\right)
{-k_X k_Y \cos \phi -
...
... a^{\rm IA}}\right)
{- k_Y \sin \phi + k_Z^{\rm TA} \cos \phi \over k^{\rm TA}}$  
    $\displaystyle + \left( {a^{\rm TM}\over a^{\rm IA}}\right) {k_X k_Y \cos \phi +...
...er {k^{\rm TM}}^2} =
{-k_Y \sin \phi + k_Z^{\rm IA} \cos \phi \over k^{\rm IA}}$ (A.79)

where the apostrophe and the dependence on the wavevector have been dropped from the wave amplitudes $a^{\alpha}$. In a similar way, the continuity of B1Y at the surface Z=0 (Eq. (A.47)), along with the expressions (A.16) and (A.19) for the perturbation magnetic field, give the equation
 
    $\displaystyle \left( {a^{\rm RA}\over a^{\rm IA}} \right)
{- k_X \sin \phi \ove...
...+
\left( {a^{\rm TA}\over a^{\rm IA}} \right)
{ k_X \sin \phi \over k^{\rm TA}}$  
    $\displaystyle \qquad\qquad+ \left( {a^{\rm TM}\over a^{\rm IA}} \right)
{ - \le...
...rm TM} \sin \phi \over {k^{\rm TM}}^2}
= {k_X \sin \phi \over k^{\rm IA}} \cdot$ (A.80)

Finally, using the continuity of the velocity perturbation (Eqs. (A.48) and (A.49)), the expressions (A.15) and (A.18) of the perturbation velocity, and the relations (A.56), (A.59), (A.63), and (A.72), we find the equations
 
    $\displaystyle \left( {a^{\rm RA}\over a^{\rm IA}} \right)
{- k_Y \sin \phi + k_...
...{\rm IA}} \right)
r {- k_Y \sin \phi + k_Z^{\rm TA} \cos \phi \over k^{\rm TA}}$  
    $\displaystyle \qquad\qquad+\left( {a^{\rm TM}\over a^{\rm IA}} \right)
{\rm sgn...
... \over k^{\rm TM}} =
{-k_Y \sin \phi + k_Z^{\rm IA} \cos \phi \over k^{\rm IA}}$ (A.81)

and
 
$\displaystyle \left( {a^{\rm RA}\over a^{\rm IA}}\right) {- k_X \over k^{\rm RA...
...m RM}} +
\left( {a^{\rm TA}\over a^{\rm IA}} \right)
r {- k_X \over k^{\rm TA}}$      
$\displaystyle +\left( {a^{\rm TM}\over a^{\rm IA}} \right)
{\rm sgn}\left(\vec{...
...TM} \cos \phi - k_Y \sin \phi \over k^{\rm TM}}
= {-k_X \over k^{\rm IA}} \cdot$     (A.82)

The four Eqs. (A.79)-(A.82) allow to determine the ratios $(a^{\rm RA}/a^{\rm IA})$, $(a^{\rm TA}/a^{\rm IA})$, $(a^{\rm RM}/a^{\rm IA})$, and $(a^{\rm TM}/a^{\rm IA})$ of the reflected and transmitted wave amplitudes to the incident wave amplitude. These ratios determine the reflection and transmission coefficients of the various waves. Knowing the wavevector of the incident wave, these amplitude ratios are calculated by numerically solving the linear system (A.79)-(A.82). In particular, the Z-components and the modulus of wavectors of transmitted and reflected waves which appear in Eqs. (A.79)-(A.82) are determined by the relations (A.57), (A.60), (A.67), (A.75) (or (A.69), (A.77) for evanescent magnetosonic waves), (A.70), and (A.78) in terms of the wavevector of the incident wave.

The system (A.79)-(A.82) is valid both for propagating and for evanescent magnetosonic waves. In the latter case, the component $k_Z^{\rm RM}$ and/or $k_Z^{\rm TM}$ are imaginary. Thus, in the general case the coefficients of Eqs. (A.79)-(A.82) are complex, and the ratios $(a^{\rm RA}/a^{\rm IA})$, $(a^{\rm TA}/a^{\rm IA})$, $(a^{\rm RM}/a^{\rm IA})$ , and $(a^{\rm TM}/a^{\rm IA})$ are complex quantities, as well.

A.5 Energy of waves

In a cold plasma, where the internal energy is neglected, the energy density is given by

 \begin{displaymath}
\epsilon={1 \over 2}\rho v^2 + {B^2 \over 8\pi}\cdot
\end{displaymath} (A.83)

We decompose physical quantities as a sum of the background and perturbation value, according to Eq. (A.1). Thus, the energy density is divided as a term due to the background structure plus the energy density $\epsilon_{\rm f}$ due to fluctuations. The latter has the following form:

 \begin{displaymath}
\epsilon_{\rm f}=\varepsilon {\vec{B}_0 \cdot \vec{B}_1 \ove...
...}\varepsilon^2 \rho_0 v_1^2 +
\varepsilon^2 {B_1^2 \over 8\pi}
\end{displaymath} (A.84)

where terms of order $O(\varepsilon^3)$ have been neglected. According to Eq. (A.6), the velocity and magnetic field perturbations are written as a sum of contributions due to the different propagating modes and expanded in Fourier integrals:
 
    $\displaystyle \vec{v}_1(\vec{x},t)=\sum_\alpha \int {{\rm d}\vec{k} \over (2\pi...
... i} \left [ \vec{k} \cdot \vec{x} -
\omega^\alpha (\vec{k})t \right ] \right \}$  
    $\displaystyle \vec{B}_1(\vec{x},t)=\sum_\alpha \int {{\rm d}\vec{k} \over (2\pi...
...i} \left [ \vec{k} \cdot \vec{x} - \omega^\alpha (\vec{k})t \right ] \right \}.$ (A.85)

Considering, for instance, the velocity perturbation, we require that it is a real quantity: $\vec{v}_1(\vec{x},t)= \vec{v}_1^*(\vec{x},t)$. Using the expression (A.85), this condition implies
 
    $\displaystyle \int {{\rm d}\vec{k} \over (2\pi)^3} \left [
\vec{v}_1^{\rm +A}(\...
...rm i} \omega^{-\rm M}(\vec{k})t)
\right ] \exp ({\rm i}\vec{k} \cdot \vec{x}) =$  
    $\displaystyle \qquad\int {{\rm d}\vec{k} \over (2\pi)^3} \left [
{\vec{v}_1^{\r...
...c{v}_1^{-\rm A}}^* (-\vec{k}) \exp ({\rm i} \omega^{-\rm A}(-\vec{k})t)
\right.$  
    $\displaystyle \left.\qquad+
{\vec{v}_1^{+\rm M}}^* (-\vec{k}) \exp ({\rm i} \om...
...m i} \omega^{-\rm M}(-\vec{k})t)
\right ] \exp ({\rm i} \vec{k} \cdot \vec{x}).$ (A.86)

The dispersion relations of Alfvénic and magnetosonic modes give the following symmetry properties:

 \begin{displaymath}
\omega^{\pm\rm A}(-\vec{k})=-\omega^{\pm\rm A}(\vec{k}) , \;...
...\vec{k})=\omega^{\pm \rm M}(\vec{k})=-\omega^{\mp M}(\vec{k}).
\end{displaymath} (A.87)

Using these relations, from Eq. (A.86) we derive the following reality conditions for the Fourier amplitudes of velocity perturbation of the different modes:

 \begin{displaymath}
\vec{v}_1^{\rm +A} (-\vec{k}) = {\vec{v}_1^{\rm +A}}^* (\vec...
...{v}_1^{\rm +M} (-\vec{k}) = {\vec{v}_1^{\rm -M}}^* (\vec{k}) .
\end{displaymath} (A.88)

The same reality conditions hold also for the Fourier amplitudes of magnetic field perturbation.

Inserting the Fourier expansions (A.85) into Eq. (A.84), the fluctuation energy density is written in the form

 
                                              $\displaystyle \epsilon_{\rm f}(\vec{x},t)={\varepsilon \over 4\pi} \vec{B}_0 \c...
...t [ {\rm i} (\vec{k}^\alpha \cdot \vec{x} - \omega^\alpha (\vec{k}) t) \right ]$  
    $\displaystyle \qquad+ \varepsilon^2 \sum_{\alpha ,\beta}
\int {{\rm d}\vec{k} \...
...t(\omega^\alpha (\vec{k}) + \omega^\beta (\vec{k}')\right)t \right ] \right \}.$ (A.89)

Assuming that the perturbation has a duration time T, we calculate the time average of the total fluctuation energy

 \begin{displaymath}
E_{\rm f} = {1 \over T} \int_{-T/2} ^{T/2} {\rm d}t
\int {\rm d}^3 \vec{x} \epsilon_{\rm f}(\vec{x},t).
\end{displaymath} (A.90)

Using the relations

 \begin{displaymath}
\int {\rm d}^3 \vec{x} \exp ( {\rm i} \vec{k} \cdot \vec{x})...
...} ^{T/2}
\exp (-{\rm i} \omega t) {\rm d}t = \delta_{\omega,0}
\end{displaymath} (A.91)

and the condition that fluctuating fields have a vanishing space average, we see that the first term of Eq. (A.89) is averaged out. Then, the average fluctuation energy can be written in the form

 \begin{displaymath}
E_{\rm f} = \varepsilon^2 \sum_{\alpha ,\beta} \int {{\rm d...
... ]
\delta_{\omega^\alpha (\vec{k}), -\omega^\beta (-\vec{k})}.
\end{displaymath} (A.92)

At any given $\vec{k}$, the frequencies $\omega^\alpha (\vec{k})$ and $\omega^\beta (\vec{k})$ of two modes are equal only when $\alpha =\beta$. Using this fact and the relations (A.87) we can express the average fluctuating energy in the form
 
    $\displaystyle E_{\rm f} = \varepsilon^2 \int {{\rm d} \vec{k}\over (2\pi)^3 }
\...
...
\vec{v}_1^{\rm +M}(\vec{k}) \cdot \vec{v}_1^{\rm -M}(-\vec{k}) \right) \right.$  
    $\displaystyle \left.\qquad +{1 \over 8\pi} \left(\vec{B}_1^{\rm +A}(\vec{k}) \c...
...vec{B}_1^{\rm +M}(\vec{k}) \cdot \vec{B}_1^{\rm -M}(-\vec{k})
\right) \right ].$ (A.93)

Finally, using the reality conditions (A.88) for the velocity perturbations and the analogue for the magnetic field perturbations in Eq. (A.93) we find the time-averaged fluctuation energy

 \begin{displaymath}
E_{\rm f} = \sum_\alpha \int {{\rm d} \vec{k}\over (2\pi)^3 ...
...{{\vert\vec{B}_1^\alpha (\vec{k})\vert}^2 \over 8 \pi} \right)
\end{displaymath} (A.94)

where

 \begin{displaymath}
E^\alpha (\vec{k}) = \varepsilon^2
\left( {1\over 2} \rho_0...
...{{\vert\vec{B}_1^\alpha (\vec{k})\vert}^2 \over 8 \pi} \right)
\end{displaymath} (A.95)

is the perturbation energy density in the $\vec{k}$-space of the $\alpha$-th mode at wavevector $\vec{k}$.

The energy of the modes involved in the reflection process can be calculated as a function of the amplitude $a_\alpha$ and of the wavevector. Considering first the incident Alfvén waves, we insert the expressions (A.18) and (A.19) of the velocity and magnetic field perturbation into Eq. (A.95), obtaining the energy density of the incident Alfvén wave:

 \begin{displaymath}
E^{\rm IA}(\vec{k}^{\rm IA})=\varepsilon^2 {B_0^2 \over 4\pi...
...Z^{\rm IA} \cos \phi\vert^2 + k_X^2 \over {k^{\rm IA}}^2}\cdot
\end{displaymath} (A.96)

In a similar way, the energy density of the reflected and transmitted Alfvén waves is

 \begin{displaymath}
E^{\rm RA,TA}\left(\vec{k}^{\rm RA,TA}\right)=\varepsilon^2 ...
... RA,TA} \cos \phi\vert^2 + k_X^2 \over {k^{\rm RA,TA}}^2}\cdot
\end{displaymath} (A.97)

The energy density of the reflected and transmitted magnetosonic waves are obtained inserting the expressions (A.15) and (A.16) of the velocity and magnetic field perturbation into Eq. (A.95):
 
    $\displaystyle E^{\rm RM,TM}\left(\vec{k}^{\rm RM,TM}\right) = \varepsilon^2 {B_...
...rt k_Z^{\rm RM,TM} \cos^2 \phi - k_Y \sin \phi \cos \phi\vert^2 \right)
\right.$  
    $\displaystyle \left.\qquad\qquad+ \vert k_X k_Y \cos \phi + k_X k_Z^{\rm RM,TM}...
...TM}}^2 - k_Y^2\right) \cos \phi - k_Y k_Z^{\rm RM,TM} \sin \phi \vert^2 \right.$  
    $\displaystyle \left.\qquad\qquad+
\vert\left({k^{\rm RM,TM}}^2 - k_Z^2\right) \sin \phi - k_Y k_Z^{\rm RM,TM} \cos \phi \vert^2
{k^{\rm RM,TM}}^4 \right ].$ (A.98)

Thus, the ratios between the reflected and transmitted wave energy density and the incident wave energy density are

 \begin{displaymath}
{E^{\rm RA,TA} \over E^{\rm IA}} = {\vert a^{\rm RA,TA}\vert...
...r
\vert k_Y \sin \phi - k_Z^{\rm IA} \cos \phi\vert^2 + k_X^2}
\end{displaymath} (A.99)

and
 
                                     $\displaystyle {E^{\rm RM,TM} \over E^{\rm IA}} = {\vert a^{\rm RM,TM}\vert^2 \o...
...^{\rm IA}}^2 \over \vert k_Y \sin \phi - k_Z^{\rm IA} \cos \phi\vert^2 + k_X^2}$  
    $\displaystyle \times\left [ {k^{\rm RM,TM}}^2 \left( k_X^2 + \vert k_Y \sin^2 \...
...rt k_Z^{\rm RM,TM} \cos^2 \phi - k_Y \sin \phi \cos \phi\vert^2 \right) \right.$  
    $\displaystyle \left.
+\vert k_X k_Y \cos \phi + k_X k_Z^{\rm RM,TM} \sin \phi \...
...}^2 - k_Z^2\right) \sin \phi - k_Y k_Z^{\rm RM,TM} \cos \phi \vert^2 \right ] .$ (A.100)

Such ratios depend on the amplitude ratios and on the wavevectors of the involved waves.

A.6 Reflection and transmission coefficients

The energy flux associated with a wave at a given wavevector $\vec{k}$, belonging to the $\alpha$-th mode can be defined as

 \begin{displaymath}
\vec{F}^\alpha (\vec{k}) = E^\alpha (\vec{k}) \vec{v}_{\rm g}^\alpha (\vec{k})
\end{displaymath} (A.101)

where $\vec{v}_{\rm g}^\alpha $ is the group velocity of the $\alpha$-th mode. In particular:

 \begin{displaymath}
\vec{v}_{\rm g}^{\pm\rm A} = \sigma_{\rm A} {\vec{B}_0 \over...
... \rho_0)^{1/2}}
{\vec{k}^{\pm \rm M} \over k^{\pm \rm M}}\cdot
\end{displaymath} (A.102)

The reflection coefficient for the Alfvén waves is given by the ratio

 \begin{displaymath}
r_{\rm A} = {\vert F_z^{\rm RA}(\vec{k}^{\rm RA})\vert \over...
... RA})\vert \over \vert E^{\rm IA}(\vec{k}^{\rm IA})\vert}\cdot
\end{displaymath} (A.103)

The transmission coefficient for the Alfvén waves is given by:

 \begin{displaymath}
t_{\rm A} = {\vert F_z^{\rm TA}(\vec{k}^{\rm TA})\vert \over...
... TA})\vert \over \vert E^{\rm IA}(\vec{k}^{\rm IA})\vert}\cdot
\end{displaymath} (A.104)

The reflection and transmission coefficients for the magnetosonic waves are

 \begin{displaymath}
r_{\rm M} = {\vert F_z^{\rm RM}(\vec{k}^{\rm RM})\vert \over...
...vert k_z^{\rm RM}\vert \over \vert k^{\rm RM} \sin \phi \vert}
\end{displaymath} (A.105)

and

 \begin{displaymath}
t_{\rm M} = {\vert F_z^{\rm TM}(\vec{k}^{\rm TM})\vert \over...
...vert k_z^{\rm TM}\vert \over \vert k^{\rm RM} \sin \phi \vert}
\end{displaymath} (A.106)

respectively. The energy density ratios in the above equations are given by the relations (A.99) and (A.100). The Eq. (A.101) holds only in the case of a non evanescent wave. If a wave is evanescent along a given direction, then the corresponding component of the energy flux $\vec{F}^\alpha$ is vanishing. Thus, when the reflected and/or transmitted magnetosonic wave is evanescent (along z), then the reflection coefficient $r_{\rm M}$ and/or the transmission coefficient $t_{\rm M}$ will be set equal to zero.



Copyright ESO 2005