© ESO, 2011

1. Introduction

Phase mixing, a mechanism for dissipating shear Alfvén waves, was first proposed by Heyvaerts & Priest (1983). The basic concept is straightforward: consider shear Alfvén waves propagating in an inhomogeneous plasma, such that on each magnetic fieldline each wave propagates with its own local Alfvén speed. Thus, after propagating a certain distance, these neighbouring perturbations will be out of phase, which will lead to the generation of smaller and smaller transverse spatial scales, and thus to the growth of strong currents. This ultimately results in a strong enhancement of the dissipation of Alfvén-wave energy via viscosity and/or resistivity. Alfvén wave phase mixing has also been studied extensively as a possible mechanism for heating the corona (e.g. Heyvaerts & Priest 1983; Browning 1991; Ireland 1996; Malara et al. 1996; Narain & Ulmschneider 1990, 1996).

Nakariakov et al. (1997) extended the model of Heyvaerts & Priest (1983) to include compressibility and nonlinear effects. They found that fast magnetoacoustic waves are generated continuously by Alfvén-wave phase mixing at a frequency twice that of the driven Alfvén wave, and that these generated waves propagate across magnetic fieldlines and away from the phase mixing layer. Since there is a permanent leakage of energy away from the phase mixing layer, these fast waves can cause indirect heating of the plasma as they propagate away and dissipate far from the layer itself, thus spreading the heating across the domain. Nakariakov et al. (1997) found that the amplitude of these fast waves grows linearly in time, according to weakly-nonlinear, β = 0, analytical theory. Nakariakov et al. (1998) further extended this model to include a background steady flow and found that the findings of their 1997 model persist.

Hood et al. (2002) investigated the phase mixing of single Alfvén-wave pulses and found that this results in a slower power-law damping (as opposed to the standard  ~exp(−t³) for harmonic Alfvén waves, i.e. Heyvaerts & Priest 1983). Hence, Alfvénic-pulse perturbations will be able to transport energy to a greater coronal height than that of a harmonic Alfvén wavetrain.

Botha et al. (2000) considered a developed stage of Alfvén-wave phase mixing and found that the growth of the generated fast waves saturates at amplitudes much lower than that of the driven Alfvén wave. They concluded that the nonlinear generation of fast waves (Nakariakov et al. 1997) saturates due to destructive wave interference, and has little effect on the standard phase mixing model of Heyvaerts & Priest (1983). The numerical simulations of Botha et al. assumed wave propagation in an ideal plasma, and thus heating was absent from their model. Tsiklauri and co-authors repeated these numerical experiments for both a weakly and strongly-nonlinear Alfvénic pulse (Tsiklauri et al. 2001, 2002).

De Moortel et al. (1999, 2000) investigated how gravitational density stratification and magnetic field divergence changes the efficiency of phase mixing. They report that the resultant dissipation can be either enhanced or diminished depending on the specific choice of equilibrium. Ruderman et al. (1998) investigated phase mixing in open magnetic equilibria, and found similar results using the WKB method. These 1998 analytical results were repeated and corrected by numerical and semi-analytical work by Smith et al. (2007). All these authors found the following: a diverging magnetic field enhances the efficiency of phase mixing, whereas gravitational stratification diminishes the mechanism. Hood et al. (1997a,b) derive analytical, self-similar solutions of Alfvén-wave phase mixing in both open and closed magnetic topologies.

Ofman & Davila (1995) found that in an inhomogeneous coronal hole with an enhanced dissipation parameter (S = 1000−10   000), Alfvén waves can dissipate within several solar radii, which can provide significant energy for the heating and acceleration of the solar wind. This model was later extended to include nonlinear effects (Ofman & Davila 1997).

Parker (1991) pointed out that phase mixing requires an ignorable coordinate, an assumption which is expected to be unphysical in the corona. Parker found that including all three coordinates results in the driven Alfvén waves coupling with a fast magnetoacoustic mode, and that this elimates the growth in transverse spatial scales, and thus phase mixing is absent from the system. Instead, such a system exhibits resonant absorption (e.g. Lee & Roberts 1986; Hollweg & Yang 1988) on the surfaces where the phase velocity equals the Alfvén velocity. However, Parker’s conclusions are in disagreement with the work of Tsiklauri and co-authors who considered the interaction of an impulsively-generated, weakly-nonlinear MHD pulse with a one-dimensional density inhomogeniety, considered in the three-dimensional regime (i.e. without an ignorable coordinate) in both an ideal (Tsiklauri & Nakariakov 2002) and resistive (Tsiklauri et al. 2003) plasma. Tsiklauri and co-authors found that phase mixing remains a relevant paradigm and that the dynamics can still be qualitatively understood in terms of the classic 2.5D models. Mocanu et al. (2008) have revisited the Heyvaerts & Priest (1983) model using anisotopic viscosity (i.e. incorporating the Braginskii 1965, stress tensor) and report that this significantly increases the damping lengths, i.e. compared to those obtained for isotropic dissipation. More recently, Threlfall et al. (2010) have investigated the effect of the Hall term on phase mixing in the ion-cyclotron range of frequencies.

Another key concept that we shall invoke in this paper is the ponderomotive force: a nonlinear force proportional to spatial gradients in magnetic pressure, also referred to as the Alfvén wave-pressure force. The ponderomotive effect has been considered in a solar context initially by Hollweg (1971) and later by Verwichte and co-authors (Verwichte 1999; Verwichte et al. 1999). Hollweg (1971) considered linearly-polarised Alfvén waves propagating in a direction parallel to the magnetic field, and found that the transverse behaviour of the Alfvén wave was identical when comparing the linear and nonlinear (to second order) solutions, but that longitudinal wave velocity and density fluctuations appear in the nonlinear solutions driven by gradients in the wave magnetic-field pressure, i.e. the Alfvén wave is no longer purely transverse and is compressive (through nonlinear coupling to magnetoacoustic waves). Thus, the ponderomotive force can be used as an extra acceleration mechanism and as an explanation for density fluctuations in the solar wind. Verwichte (1999) presented a mechanical analogy for the ponderomotive effect by considering the resulting motion of discrete particles (beads) on an oscillating string. Verwichte et al. (1999) considered the temporal evolution of a weakly-nonlinear, Alfvén wave in a β = 0 homogeneous plasma. These authors showed that the an initially-excited, gaussian-pulse perturbation in transverse velocity splits into two Alfvén wave pulses, each propagating in opposite directions (as naturally expected). Furthermore, Verwichte et al. found that the ponderomotive force produces a shock in longitudinal velocity at the starting position. Note that in a cold plasma, there is no force to counteract this ponderomotive acceleration.

In this paper, we investigate the nonlinear, nonideal behaviour of Alfvén-wave propagation and phase mixing over long timescales. Thus, this paper can be seen as an extension of the model of Botha et al. (2000) to include visco-resistive effects. We also seek to address a fundamental question of phase mixing: by considering nonlinear, nonideal phase mixing over long timescales, is it possible to observe a drifting of the heating layer? In other words, phase mixing will occur due to the density inhomogenity in our system, and this process will generate strong, localised heating due to enhanced dissipation. This localised heat deposition is expected to modify the equilibrium density profile, and thus may change the location of maximum heating. However, it is unclear what the actual result will be: we may observe a change in the location of the heating layer, or the phase mixing mechanism may break down, or the heating may bifurcate spatially (since our density profile will no longer be monotonic). Is is also unclear how this may affect the indirect heating of the plasma, due to the coupling to the fast magnetoacoustic mode (Nakariakov et al. 1997). Such an investigation requires nonlinear and nonideal effects to be considered and, as we shall see, observed over long timescales.

The work of Ofman et al. (1998) is also relevant here. These authors investigated a model of resonant absorption that incorporated the dependence of loop density on the heating rate, and studied the spatial and temporal dependence of the heating layer. Ofman et al. find that the heating occurs in multiple resonance layers, rather than the single layer of the classic resonant absorption models (e.g. Ionson 1978; Ulmschneider et al. 1991; Ruderman & Roberts 2002) and that these layers drift throughout the loop to heat the entire volume. Poedts & Boynton (1996) also investigated resonant absorption using nonlinear, resistive MHD simulations and found a spreading of the heat deposition, i.e. a broadening of the resonant layer due to changes in the background inhomogeneity.

This paper has the following outline: Sect. 2 describes the governing equations, assumptions and analytical and numerical details of our investigation, Sect. 3 investigates the 1D nonlinear, nonideal system (with no density inhomogeneity) and focuses on the underlying physical processes. Section 4 details the bulk-flow phenomenon found to be present in our system, and investigates its density dependence. Section 5 considers a 2D model with an inhomogeneous density profile, and details the long-term evolution and coupled nature of the three MHD waves present in our system. The conclusions are presented in Sect. 6 and there are two appendicies.

2. Basic equations

We consider the nonlinear, compressible, viscous and resistive MHD equations appropriate to the solar corona: $\begin{array}{ccc}\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left(\mathit{\rho }{v}\right)& \mathrm{=}& \mathrm{0}\mathit{,}\\ \mathit{\rho }\left[\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\left({v}\mathrm{·}\mathrm{\nabla }\right){v}\right]& \mathrm{=}& \mathrm{-}\mathrm{\nabla }\mathit{p}\mathrm{+}\left(\frac{\mathrm{\nabla }\mathrm{\times }{B}}{\mathit{\mu }}\right)\mathrm{\times }{B}\mathrm{+}\mathit{\nu }\mathrm{\nabla }\mathrm{·}{S}\mathit{,}\\ \frac{\mathit{\partial }{B}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{\nabla }\mathrm{\times }\mathrm{(}{v}\mathrm{\times }{{B}}^{\mathrm{)}}\mathrm{+}\mathit{\eta }{\mathrm{\nabla }}^{\mathrm{2}}{B}\mathit{,}\\ \frac{\mathit{\partial p}}{\mathit{\partial t}}\mathrm{+}\left({v}\mathrm{·}\mathrm{\nabla }\right)\mathit{p}& \mathrm{=}& \end{array}$(1)where ρ is the mass density, v is the plasma velocity, B the magnetic induction (usually called the magnetic field), p is the plasma pressure, μ = 4π × 10^-7Hm^-1 is the magnetic permeability, σ is the electrical conductivity, η = 1 / μσ is the magnetic diffusivity, γ = 5 / 3 is the ratio of specific heats and j = ∇ × B / μ is the electric current density. In this paper, η and ν are assumed to be constants.

The viscous stress tensor, S, is given by $\begin{array}{ccc}{\mathit{S}}_{\mathit{ij}}\mathrm{=}{\mathit{ϵ}}_{\mathit{ij}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{\nabla }\mathrm{·}{v}\mathit{,}& & \end{array}$where ϵ is the rate-of-strain tensor, given by: $\begin{array}{ccc}{\mathit{ϵ}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{i}}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\mathrm{+}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{j}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\right)\mathrm{·}& & \end{array}$The viscous heating term is given by $\begin{array}{ccc}{\mathit{Q}}_{\mathrm{visc}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathrm{3}}{\mathit{ϵ}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}\mathit{.}& & \end{array}$Note that in the presence of strong magnetic fields, the classical viscous term used in Eqs. (1) is not the most appropriate for the solar corona, since the viscosity takes the form of a non-isotropic tensor. The mathematical details of this non-isotropic tensor can be found in Braginskii (1965). In this paper, we have chosen to implement the simplest, isotropic version of the stress tensor as a first step in investigating nonlinear phase mixing, and non-isotropic viscosity will be considered in future investigations. Our choice of viscosity term also allows a direct comparison with the results of Heyvaerts & Priest (1983) and allows us to derive analytical solutions to our governing equations.

We now consider a change of scale to non-dimensionalise all variables. Let v = v₀v ^∗ , B = BB ^∗ , x = Lx ^∗ , y = Ly ^∗ , ρ = ρ₀ρ ^∗ , p = p₀p ^∗ , $\mathrm{\nabla }\mathrm{=}\frac{\mathrm{1}}{\mathit{L}}{\mathrm{\nabla }}^{\mathrm{\ast }}$, t = τ_At ^∗ , η = η₀ and ν = ν₀, where  ∗  denotes a dimensionless quantity and v₀, B, L, ρ₀, p₀, τ_A, η₀ and ν₀ are constants with the dimensions of the variable they are scaling. We then set $\mathit{B}\mathit{/}\sqrt{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\mathrm{=}{\mathit{v}}_{\mathrm{0}}$ and v₀ = L / τ_A (this sets v₀ as a constant equilibrium Alfvén speed). We also set ${\mathit{\eta }}_{\mathrm{0}}{\mathit{\tau }}_{\mathrm{A}}\mathit{/}{\mathit{L}}^{\mathrm{2}}\mathrm{=}{\mathit{R}}_{\mathrm{m}}^{-1}$, where R_m is the magnetic Reynolds number. This process non-dimensionalises Eqs. (1) and under these scalings, t ^∗  = 1 (for example) refers to t = τ_A = L / v₀; i.e. the time taken to travel a distance L at the background Alfvén speed. For the rest of this paper, we drop the star indices; the fact that all variables are now non-dimensionalised is understood.

2.1. Basic equilibrium

We consider Eqs. (1) in Cartesian coordinates and assume there are no variations in the z-direction (∂ / ∂z = 0). We consider an inhomogeneous plasma in the x-direction embedded in a uniform magnetic field in the y-direction, i.e.: $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{0}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathrm{=}\mathrm{constant}\mathit{,}{{B}}_{\mathrm{0}}\mathrm{=}\mathrm{\left(}\mathrm{0}\mathit{,}{\mathit{B}}_{\mathrm{0}}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{,}{{v}}_{\mathrm{0}}\mathrm{=}{0}& & \end{array}$(2)with finite amplitude perturbations of the form: $\begin{array}{ccc}\mathit{\rho }& \mathrm{=}& \mathit{\rho }\mathrm{\left(}\mathit{x,y,t}\mathrm{\right)}\mathit{,}{v}\mathrm{=}\mathrm{\left(}{\mathit{v}}_{\mathit{x}}\mathit{,}{\mathit{v}}_{\mathit{y}}\mathit{,}{\mathit{v}}_{\mathit{z}}\mathrm{\right)}\mathrm{=}{v}\mathrm{\left(}\mathit{x,y,t}\mathrm{\right)}\\ \mathit{p}& \mathrm{=}& \end{array}$(3)

2.2. Analytical description

Following the analysis of Nakariakov et al. (1995, 1997) and Botha et al. (2000), Eqs. (1) are written in component form using equilibrium (2) and perturbations (3): $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial x}}\left({\mathit{\rho }}_{\mathrm{0}}{\mathit{v}}_{\mathit{x}}\right)\mathrm{+}{\mathit{\rho }}_{\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial y}}\mathrm{=}{\mathit{N}}_{\mathrm{1}}\mathit{,}\\ & & {\mathit{\rho }}_{\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial p}}{\mathit{\partial x}}\mathrm{+}\frac{{\mathit{B}}_{\mathrm{0}}}{\mathit{\mu }}\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{y}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{x}}}{\mathit{\partial y}}\right)\mathrm{+}{\mathit{V}}_{\mathrm{1}}\mathrm{=}{\mathit{N}}_{\mathrm{2}}\mathit{,}\\ & & {\mathit{\rho }}_{\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial p}}{\mathit{\partial y}}\mathrm{+}{\mathit{V}}_{\mathrm{2}}\mathrm{=}{\mathit{N}}_{\mathrm{3}}\mathit{,}\\ & & {\mathit{\rho }}_{\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial t}}\mathrm{-}\frac{{\mathit{B}}_{\mathrm{0}}}{\mathit{\mu }}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\mathrm{+}{\mathit{V}}_{\mathrm{3}}\mathrm{=}{\mathit{N}}_{\mathrm{4}}\mathit{,}\\ & & \frac{\mathit{\partial }{\mathit{B}}_{\mathit{x}}}{\mathit{\partial t}}\mathrm{-}{\mathit{B}}_{\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial y}}\mathrm{-}{\mathit{R}}_{\mathrm{1}}\mathrm{=}{\mathit{N}}_{\mathrm{5}}\mathit{,}\\ & & \frac{\mathit{\partial }{\mathit{B}}_{\mathit{y}}}{\mathit{\partial t}}\mathrm{+}{\mathit{B}}_{\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x}}\mathrm{-}{\mathit{R}}_{\mathrm{2}}\mathrm{=}{\mathit{N}}_{\mathrm{6}}\mathit{,}\\ & & \frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial t}}\mathrm{-}{\mathit{B}}_{\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial y}}\mathrm{-}{\mathit{R}}_{\mathrm{3}}\mathrm{=}{\mathit{N}}_{\mathrm{7}}\mathit{,}\\ & & \frac{\mathit{\partial p}}{\mathit{\partial t}}\mathrm{+}\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial y}}\right)\mathrm{=}{\mathit{N}}_{\mathrm{8}}\mathit{,}\end{array}$where the left-hand-sides involve only linear terms and N₁ − N₈ denote the nonlinear terms. The nonlinear terms are: $\begin{array}{ccc}{\mathit{N}}_{\mathrm{1}}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial x}}\left(\mathit{\rho }{\mathit{v}}_{\mathit{x}}\right)\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial y}}\mathrm{(}\mathit{\rho }{{\mathit{v}}_{\mathit{y}}}^{\mathrm{)}}\mathit{,}\\ {\mathit{N}}_{\mathrm{2}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{B}}_{\mathit{z}}}{\mathit{\mu }}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial x}}\mathrm{-}\mathit{\rho }\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial t}}\mathrm{-}\left({\mathit{\rho }}_{\mathrm{0}}\mathrm{+}\mathit{\rho }\right)\left({\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{+}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }}{\mathit{\partial y}}\right){\mathit{v}}_{\mathit{x}}\\ & & \mathrm{-}\frac{{\mathit{B}}_{\mathit{y}}}{\mathit{\mu }}\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{y}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{x}}}{\mathit{\partial y}}\right)\mathit{,}\\ {\mathit{N}}_{\mathrm{3}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{B}}_{\mathit{z}}}{\mathit{\mu }}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\mathrm{-}\mathit{\rho }\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial t}}\mathrm{-}\left({\mathit{\rho }}_{\mathrm{0}}\mathrm{+}\mathit{\rho }\right)\left({\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{+}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }}{\mathit{\partial y}}\right){\mathit{v}}_{\mathit{y}}\\ & & \mathrm{+}\frac{{\mathit{B}}_{\mathit{x}}}{\mathit{\mu }}\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{y}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{x}}}{\mathit{\partial y}}\right)\mathit{,}\\ {\mathit{N}}_{\mathrm{4}}& \mathrm{=}& \mathrm{-}\mathit{\rho }\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial t}}\mathrm{-}\left({\mathit{\rho }}_{\mathrm{0}}\mathrm{+}\mathit{\rho }\right)\left({\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{+}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }}{\mathit{\partial y}}\right){\mathit{v}}_{\mathit{z}}\\ & & \mathrm{+}\frac{{\mathit{B}}_{\mathit{y}}}{\mathit{\mu }}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\mathrm{+}\frac{{\mathit{B}}_{\mathit{x}}}{\mathit{\mu }}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial x}}\mathit{,}\\ {\mathit{N}}_{\mathrm{5}}& \mathrm{=}& \frac{\mathit{\partial }}{\mathit{\partial y}}\mathrm{(}{\mathit{v}}_{\mathit{x}}{\mathit{B}}_{\mathit{y}}\mathrm{-}{\mathit{v}}_{\mathit{y}}{\mathit{B}}_{\mathit{x}}\mathrm{)}\mathit{,}\\ {\mathit{N}}_{\mathrm{6}}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{(}{\mathit{v}}_{\mathit{x}}{\mathit{B}}_{\mathit{y}}\mathrm{-}{\mathit{v}}_{\mathit{y}}{\mathit{B}}_{\mathit{x}}\mathrm{)}\mathit{,}\\ {\mathit{N}}_{\mathrm{7}}& \mathrm{=}& \frac{\mathit{\partial }}{\mathit{\partial x}}\left({\mathit{v}}_{\mathit{z}}{\mathit{B}}_{\mathit{x}}\mathrm{-}{\mathit{v}}_{\mathit{x}}{\mathit{B}}_{\mathit{z}}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial y}}\mathrm{(}{\mathit{v}}_{\mathit{z}}{\mathit{B}}_{\mathit{y}}\mathrm{-}{\mathit{v}}_{\mathit{y}}{\mathit{B}}_{\mathit{z}}\mathrm{)}\mathit{,}\\ {\mathit{N}}_{\mathrm{8}}& \mathrm{=}& \mathrm{-}\left({\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{+}{\mathit{v}}_{\mathit{y}}\frac{\mathit{\partial }}{\mathit{\partial y}}\right)\mathit{p}\mathrm{-}\mathit{\gamma p}\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial y}}\right)\mathrm{+}{\mathit{R}}_{\mathrm{4}}\mathrm{+}{\mathit{V}}_{\mathrm{4}}\mathit{,}\end{array}$the resistive terms (where R₁ − R₃ are linear, R₄ is nonlinear) are: $\begin{array}{ccc}{\mathit{R}}_{\mathrm{1}}& \mathrm{=}& \mathit{\eta }\left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{B}}_{\mathit{x}}\mathit{,}\\ {\mathit{R}}_{\mathrm{2}}& \mathrm{=}& \mathit{\eta }\left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{B}}_{\mathit{y}}\mathit{,}\\ {\mathit{R}}_{\mathrm{3}}& \mathrm{=}& \mathit{\eta }\left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{B}}_{\mathit{z}}\mathit{,}\\ {\mathit{R}}_{\mathrm{4}}& \mathrm{=}& \left(\mathit{\gamma }\mathrm{-}\mathrm{1}\right)\frac{\mathit{\eta }}{\mathit{\mu }}\left[{\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{y}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{x}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial x}}\right)}^{\mathrm{2}}\right]\mathit{,}\end{array}$and the viscous terms (V₁ − V₃ are linear, V₄ is nonlinear) are $\begin{array}{ccc}{\mathit{V}}_{\mathrm{1}}& \mathrm{=}& \\ {\mathit{V}}_{\mathrm{2}}& \mathrm{=}& \\ {\mathit{V}}_{\mathrm{3}}& \mathrm{=}& \\ {\mathit{V}}_{\mathrm{4}}& \mathrm{=}& \mathit{\nu }\left(\mathit{\gamma }\mathrm{-}\mathrm{1}\right)\{{\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x}}\right)}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial y}}\mathrm{+}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial x}}\right)}^{\mathrm{2}}\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial x}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\left[{\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\right]\left(\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial y}}\right)\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathrm{·}\end{array}$Equations (4) − (11) describe all the waves occurring in the plasma. These equations can be combined to obtain the evolution expressions for the velocity components: $\begin{array}{ccc}& & \left[\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}\mathrm{(}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}{{\mathit{c}}_{\mathit{S}}^{\mathrm{2}}}^{\mathrm{)}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{-}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right]{\mathit{v}}_{\mathit{x}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{v}}_{\mathit{y}}}{\mathit{\partial x\partial y}}\\ & & \mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}[\mathrm{-}\frac{\mathit{\partial }{\mathit{V}}_{\mathrm{1}}}{\mathit{\partial t}}\mathrm{-}\frac{\mathit{\partial }{\mathit{N}}_{\mathrm{8}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }{\mathit{N}}_{\mathrm{2}}}{\mathit{\partial t}}\\ & & \mathrm{-}\frac{{\mathit{B}}_{\mathrm{0}}}{\mathit{\mu }}\left(\frac{\mathit{\partial }{\mathit{N}}_{\mathrm{6}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{2}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{N}}_{\mathrm{5}}}{\mathit{\partial y}}\mathrm{-}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{1}}}{\mathit{\partial y}}\right)]\mathit{,}\\ & & \left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{v}}_{\mathit{y}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x\partial y}}\\ & & \mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\left(\mathrm{-}\frac{\mathit{\partial }{\mathit{V}}_{\mathrm{2}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }{\mathit{N}}_{\mathrm{3}}}{\mathit{\partial t}}\mathrm{-}\frac{\mathit{\partial }{\mathit{N}}_{\mathrm{8}}}{\mathit{\partial y}}\right)\mathit{,}\\ & & \left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{v}}_{\mathit{z}}\\ & & \end{array}$where ${\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}$ is the unperturbed Alfvén speed and ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}$ is the unperturbed sound speed. In Eqs. (28) − (30) the nonlinear and nonideal terms can be found on the right-hand-side. Note that, Eqs. (4) − (30) reduce to the corresponding equations in Botha et al. (2000) in an ideal plasma.

Considering Eq. (30) and neglecting the nonlinear and nonideal terms (right-hand-side), we see that v_z corresponds to the linear Alfvén wave. Thus, if v_A = v_A(x), then both v_z and B_z will depend on x. Hence, there will be coupling to both magnetoacoustic modes as follows: Firstly, we can see that the first term in N₂ (i.e.  ~B_z∂B_z / ∂x) will generate a fast wave (i.e. drives v_x term) via phase mixing. Secondly, we see that the first term in N₃ (i.e.  ~B_z∂B_z / ∂y) corresponds to the Alfvén-wave pressure gradient, or the ponderomotive force, and that this (nonlinear) force generates a slow wave (i.e. drives v_y term) along the background magnetic field. Also note that these two terms we have described are quadratic in terms of B_z and so will have an amplitude that is related to the square of the amplitude of the linear Alfvén wave.

Fig. 1 a) Longitudinal perturbations vy for β = 0.1 plasma with η = ν = 0. Dotted and dashed lines represent y = vAt and y = cst respectively. The blue line denotes the maximum amplitude: $\mathit{y}\mathrm{=}{\mathit{A}}^{\mathrm{2}}\mathit{/}\mathrm{4}{}^{\mathrm{(}}\mathit{v}_{\mathrm{A}}^{\mathrm{2}}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}$. b) Blow-up of region 100 ≤ y ≤ 200 for vy. In both subfigures, t = 1000τA.

2.3. Numerical Set-up

We solve MHD Eqs. (1) using the LARE2D Lagrangian-remap code (Arber et al. 2001).

We consider a numerical domain with  − 100 ≤ x ≤ 100, 0 ≤ y ≤ 10   000 with uniform-grid spacing in the x-direction and a stretched grid in the y-direction. The stretched grid places the majority of the y-direction gridpoints near low values of y. Results presented in this paper have a typical numerical resolution of 2000 × 20   000 which means that δx ≈ δy ≈ 0.1 (due to the stretched grid). We chose such a large domain to ensure that no wave energy is reflected back into our system during our numerical simulations.

We drive our system with linearly-polarised, sinusoidal Alfvén waves. This harmonic wave train is driven at the y = 0 boundary such that: $\begin{array}{ccc}{\mathit{v}}_{\mathit{z}}\mathrm{\left(}\mathit{x,}\mathrm{0}\mathit{,t}\mathrm{\right)}& \mathrm{=}& \mathit{A}\sin \mathrm{(}\mathit{\omega }{\mathit{t}}^{\mathrm{)}}\mathit{,}{\mathit{v}}_{\mathit{x}}\mathrm{\left(}\mathit{x,}\mathrm{0}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathit{y}}\mathrm{\left(}\mathit{x,}\mathrm{0}\mathit{,t}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ {\mathit{B}}_{\mathit{z}}\mathrm{\left(}\mathit{x,}\mathrm{0}\mathit{,t}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\mathit{A}\sqrt{\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}}\sin \mathrm{(}\mathit{\omega }{\mathit{t}}^{\mathrm{)}}\mathit{,}{\frac{\mathit{\partial }{\mathit{B}}_{\mathit{x}}}{\mathit{\partial y}}|}_{\mathit{y}\mathrm{=}\mathrm{0}}\mathrm{=}{\frac{\mathit{\partial }{\mathit{B}}_{\mathit{y}}}{\mathit{\partial y}}|}_{\mathit{y}\mathrm{=}\mathrm{0}}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(31)where A is the amplitude and ω is the frequency. All other quantities have zero gradient boundary conditions at the y-boundaries, and we utilise periodic boundary conditions at the x-boundaries. In this paper, we set A = 0.01 and ω = 0.25. At the start of the numerical run, the driven Alfvén wave is ramped up to A over the first four wavelengths. It should be noted that boundary conditions (31) do not drive a pure Alfvén wave into the numerical domain; slow magnetoacoustic modes are nonlinearly excited as well with amplitude of order A² (see Sect. 3.1).

Note that our choice of equilibrium (Eq. (2)) is formally only an equilibrium in ideal MHD. However, numerical tests with no driving term show that this equilibrium is still valid over long timescales, i.e. the results that follow are due to the driven Alfvén waves, not due to our choice of equilibrium.

3. One-dimensional system

Let us start by considering the nonlinear and nonideal aspects of the Alfvén wave. In this section, we set ρ₀ = constant, i.e. we remove the density inhomogeneity. This effectively reduces Eqs. (4) − (30) to a one-dimensional (1D) system where variations in the x-direction are ignorable (∂ / ∂x = 0). This 1D system will clearly illustrate the nonlinear driver in the equations as well as the effect of viscosity and resistivity. These effects must be identified before trying to interpret the results with a density inhomogeneity.

The driven, linear Alfvén wave has the form: ${\mathit{v}}_{\mathit{z}}\mathrm{=}\mathrm{ℱ}\mathrm{\left(}\mathit{\omega }\mathrm{\right[}\mathit{t}\mathrm{-}\mathit{y}\mathit{/}{\mathit{v}}_{\mathrm{A}}\mathrm{\left]}\mathrm{\right)}\mathit{,} {\mathit{B}}_{\mathit{z}}\mathrm{=}\mathrm{-}\sqrt{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}{\mathit{v}}_{\mathit{z}}\mathit{,}$(32)where the arbitrary function ℱ can include a form for the ramp-up of the driver as well as the periodic driver. To illustrate the terms, we ignore the ramp-up and consider the sinusodial driver; the solutions are: $\begin{array}{ccc}{\mathit{v}}_{\mathit{z}}& \mathrm{=}& \mathit{A}\sin \left(\mathit{\omega }\mathrm{[}\mathit{t}\mathrm{-}\mathit{y}\mathit{/}{\mathit{v}}_{\mathrm{A}}\mathrm{]}\right) \mathrm{for} \mathrm{0}\mathit{<}\mathit{y}\mathit{<}{\mathit{v}}_{\mathrm{A}}\mathit{t,}\\ {\mathit{B}}_{\mathit{z}}& \mathrm{=}& \mathrm{-}\sqrt{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\mathit{A}\sin \left(\mathit{\omega }\mathrm{[}\mathit{t}\mathrm{-}\mathit{y}\mathit{/}{\mathit{v}}_{\mathrm{A}}\mathrm{]}\right) \mathrm{for} \mathrm{0}\mathit{<}\mathit{y}\mathit{<}{\mathit{v}}_{\mathrm{A}}\mathit{t,}\end{array}$and $\begin{array}{ccc}{\mathit{v}}_{\mathit{z}}\mathrm{=}{\mathit{B}}_{\mathit{z}}\mathrm{=}\mathrm{0} \mathrm{for} \mathit{y}\mathit{>}{\mathit{v}}_{\mathrm{A}}\mathit{t}\mathit{.}& & \end{array}$The nonlinear evolution of the Alfvén wave in an ideal, β = 0 plasma can be found in Appendix A. However, this result can be easily derived from the β ≠ 0 case (Sect. 3.1) and so Appendix A is only included for completeness.

3.1. Ideal β > 0 plasma

Consider a β ≠ 0 plasma, where we set β = 0.1. The addition of finite-β has no effect on the v_z component of the driven Alfvén wave (i.e. Eq. (32) is still valid) but does influence the longitudinal component of velocity generated by the ponderomotive force. This can be seen in Fig. 1a. Here, we see that there are two kinds of longitudinal motions in the system: the first can be seen between y = c_st ≈ 280 and y = v_At = 1000. These longitudinal motions have a maximum amplitude of ${\mathit{A}}^{\mathrm{2}}\mathit{/}\mathrm{4}{}^{\mathrm{(}}\mathit{v}_{\mathrm{A}}^{\mathrm{2}}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}$ (this maximum is denoted by the horizontal blue line) and they are always positive. Thus, there is a net outflow.

The second type of longitudinal motion can be seen between y = 0 and y = c_st. These longitudinal perturbations are acoustic waves (in 1D) with both positive and negative motions of the plasma along the field, and can be identified as a boundary-driven slow wave. Note that the longitudinal motions for 0 ≤ y ≤ c_st are a combination of these two motions: the acoustic/slow wave (speed c_s) and the ponderomotive perturbations (speed v_A).

In Fig. 1b, we can see a blow-up of the region 100 ≤ x ≤ 200 for v_y (black line). The red line represents an analytical solution describing both types of longitudinal motions (see Eqs. (33) below). The agreement is excellent. Note that the slow wave has only propagated a distance of y = c_st at t = 1000τ_A and so, after y = c_st ≈ 280, only one wave is present.

The analytical solution for these nonlinearly-generated longitudinal motions is found by substituting Eq. (32) into the weakly-nonlinear form of Eq. (29), namely: $\begin{array}{ccc}\left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{v}}_{\mathit{y}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{\partial }{\mathit{N}}_{\mathrm{3}}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial t\partial y}}\left(\frac{{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{·}& & \end{array}$Thus, the analytical solution for v_y is: $\begin{array}{ccc}{\mathit{v}}_{\mathit{y}}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \end{array}& & \end{array}$(33)which is valid only for early times (i.e. until nonlinear response becomes non-negligible). Note that these longitudinal motions were first reported in Botha et al. (2000).

Note that for 0 < y < c_st, both types of wave in Eq. (33) must have the same frequency (i.e. twice the driving frequency) but will propagate at different speeds. Thus, they must have different wavenumbers and hence different wavelengths. For the parameters chosen in this paper, these wavelengths are: $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{alfv}\begin{array}{}\mathrm{´}\\ \mathrm{e}\end{array}\mathrm{n}}\mathrm{=}\mathrm{12.6} \mathrm{and} {\mathit{\lambda }}_{\mathrm{slow}}\mathrm{=}\mathrm{3.63.}& & \end{array}$Thus, for the typical numerical resolutions considered in this paper (δy ~ 0.1), both these perturbations are well resolved.

Fig. 2 a) Transverse perturbations vz for η = 0.01, β = 0.1 plasma (ν = 0), where dashed lines represent the envelope vz =  ± Aexp(kIy). b) Longitudinal perturbations vy for same plasma, where dashed line represents y = cst and blue horizontal line denotes vy = 0. In both subfigures, t = 1000τA and dotted line represents y = vAt.

3.2. Resistive plasma

Let us now consider a resistive plasma, where we set η = 0.01 (and ν = 0, β = 0.1). As before, an Alfvén wave is driven into the numerical domain and the resultant v_z behaviour can be seen in Fig. 2a. It clear that the wave now experiences resistive damping. Such damping can be estimated from Fourier analysing the linear form of Eq. (30) to give the dispersion relation: $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}\mathrm{i}\mathit{\eta }{\mathit{\omega }}^{\mathrm{)}}{\mathit{k}}^{\mathrm{2}}& & \end{array}$(34)where ω = 0.25, ${\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{=}{\mathrm{\left|}{{B}}_{\mathrm{0}}\mathrm{\right|}}^{\mathrm{2}}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}$, η = 0.01 and k = k_R + ik_I. In Fig. 2a, the dashed lines represent the envelope v_z =  ± Aexp(k_Iy), where k_I < 0. Assuming ${\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{\gg }\mathit{\eta \omega }$ (or equivalently that the magnetic Reynolds number is large) we can estimate that: $\begin{array}{ccc}{\mathit{k}}_{\mathrm{I}}\mathrm{=}\mathrm{-}\frac{\mathit{\eta }{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{2}{\mathit{v}}_{\mathrm{A}}^{\mathrm{3}}}& & \end{array}$(35)to a high degree of accuracy.

Figure 2b shows the corresponding longitudinal motions (v_y) in our resistive system. As in Fig. 1b, we see that the acoustic and ponderomotive wave components are again present, but that there is now a third phenomenon: a bulk flow in the positive y-direction with v_y > 0 and that has a maximum around y = c_st. Physically, this bulk flow (i.e. movement of material) is due to the ohmic heating in our system, and is a natural consequence of driving an Alfvén wave in a resistive plasma.

The longitudinal behaviour is governed by Eq. (29), which in 1D (i.e. ∂ / ∂x = 0) and at early times (i.e. no nonlinear feedback) reduces to: $\begin{array}{ccc}& & \left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{v}}_{\mathit{y}}\mathrm{=}\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{B}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\right)\mathrm{-}\frac{\left(\mathit{\gamma }\mathrm{-}\mathrm{1}\right)\mathit{\eta }}{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial y}}\left[{\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\right]\mathrm{·}\end{array}$(36)The right-hand-side of this equation has two contributions: the first term depends upon the rate-of-change of the magnetic pressure gradient $\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial t\partial y}}\mathrm{(}\frac{{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}}$) of the driven wave and the second depends upon pressure gradients created by the ohmic heating term (i.e. $\frac{\mathit{\partial }}{\mathit{\partial y}}{\mathit{j}}_{\mathit{x}}^{\mathrm{2}}$) increasing the gas pressure. Note the second contribution vanishes under the ideal approximation.

Let us now evaluate these terms under the assumption that our linear Alfvén wave can be represented as: $\begin{array}{ccc}{\mathit{v}}_{\mathit{z}}\mathrm{=}\mathit{A}\sin \left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right){\mathrm{e}}^{{\mathit{k}}_{\mathrm{I}}\mathit{y}}\mathit{.}& & \end{array}$(37)Thus, from Eqs. (7) and (10), B_z has the form: $\begin{array}{ccc}{\mathit{B}}_{\mathit{z}}& \mathrm{=}& \mathrm{-}\frac{\mathit{\omega }{\mathit{B}}_{\mathrm{0}}\mathit{A}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{k}}_{\mathrm{R}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{I}}^{\mathrm{2}}\mathrm{)}}[{\mathit{k}}_{\mathrm{R}}\sin \left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)\\ & & \mathrm{-}{\mathit{k}}_{\mathrm{I}}\cos \left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)]{\mathrm{e}}^{{\mathit{k}}_{\mathrm{I}}\mathit{y}}\mathit{.}\end{array}$(38)The derivation of Eq. (38) is given in Appendix B.

Fig. 3 a) Longitudinal perturbations vy for η = 0.01, β = 0.1 plasma comparing η = 0.01, ν = 0 (blue) and ν = 0.01, η = 0 (red). Note only 0 ≤ y ≤ 350 is shown and that the agreement after y ≈ 280 is excellent (the two curves lie on each other). b) Blow-up of vy over region 0 ≤ y ≤ 100 for ν = 0.01, η = 0 plasma only, where crosses indicate grid point resolution. Black line indicates vy = 0. In both subfigures, t = 1000τA.

Fig. 4 a) Density evolution in ideal plasma. b) Density evolution in nonideal plasma (η = ν = 0.01). In both subfigures, β = 0.1, t = 1000τA, dashed (dotted) lines represent y = cst (y = vAt) and the blue line denotes ρ = ρ0.

Hence, the right-hand-side of Eq. (36) can be simplified such that: $\begin{array}{ccc}\left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{v}}_{\mathit{y}}\mathrm{=}\mathit{P̂}\mathrm{+}\mathit{R̂}\mathrm{+}\mathit{BF̂}& & \end{array}$where the terms representing the Ponderomotive force, , the resistive heating force, , and the non-oscillatory resistive bulk flow force, , are given by: $\begin{array}{ccc}\mathit{P̂}& \mathrm{=}& \frac{{\mathit{\omega }}^{\mathrm{3}}{\mathit{A}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{k}}_{\mathrm{R}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{I}}^{\mathrm{2}}\mathrm{)}}[{\mathit{k}}_{\mathrm{R}}\cos \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)\\ & & \mathrm{+}{\mathit{k}}_{\mathrm{I}}\sin \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)]{\mathrm{e}}^{\mathrm{2}{\mathit{k}}_{\mathrm{I}}\mathit{y}}\\ \mathit{R̂}& \mathrm{=}& \mathrm{-}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{A}}^{\mathrm{2}}\left(\mathit{\gamma }\mathrm{-}\mathrm{1}\right)\mathit{\eta }}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}}[{\mathit{k}}_{\mathrm{R}}\sin \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)\\ & & \mathrm{+}{\mathit{k}}_{\mathrm{I}}\cos \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)]{\mathrm{e}}^{\mathrm{2}{\mathit{k}}_{\mathrm{I}}\mathit{y}}\\ \mathit{BF̂}& \mathrm{=}& \mathrm{-}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{A}}^{\mathrm{2}}\left(\mathit{\gamma }\mathrm{-}\mathrm{1}\right)\mathit{\eta }{\mathit{k}}_{\mathrm{I}}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{2}{\mathit{k}}_{\mathrm{I}}\mathit{y}}\end{array}$which is valid for c_st < y < v_At. The derivation of these equations is given in Appendix B.

Note that the amplitude of each term is proportional to the amplitude squared of the driven Alfvén wave and that they decay with height since k_I < 0. In addition, we note that and are trigonometric terms with twice the driven frequency and twice the wavenumber. In the absence of dissipation (η = 0 and so k_I = 0) there is no decay and and vanish.

We also note that is, critically, non-oscillatory and it is this term that will create the bulk flow in the longitudinal direction.

Fig. 5 Long-term evolution of a) density, b) temperature, c) pressure gradient (∂p / ∂y), and d) kI, over 0 ≤ t ≤ 100   000τA, where the different colours represent different times. The colour bar denotes the different times.

3.3. Viscous plasma

Let us now consider a purely viscous plasma, where we set ν = 0.01, η = 0 and β = 0.1. Again, we drive an Alfvén wave into our numerical domain and the resultant v_z is identical to that seen Fig. 2a. However, the damping mechanism is now due to viscosity rather than resistivity, as in Sect. 3.2. Such damping can be estimated from Fourier analysing the linear form of Eq. (30) to give the dispersion relation: $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}\mathrm{i}\mathit{\nu }{\mathit{\omega }}^{\mathrm{)}}{\mathit{k}}^{\mathrm{2}}\mathit{.}& & \end{array}$(42)Note the similarity to Eq. (34).

Figure  3a shows a comparison of v_y for the purely-resistive (blues, Sect. 3.2, Fig. 2b) and purely-viscous (red line) systems. Focusing on the viscous propagation, we see that three components are present: the ponderomotive wave (between 280 ≤ y ≤ 1000), the bulk flow in the positive y-direction, and the acoustic perturbations (from 0 ≤ y ≤ 280). However, note that the acoustic perturbations now take a different form to those in the purely-resistive system (black line) and we note that in the viscous system the smaller wavelengths are rapidly damped out by viscosity. Above y ≈ 280, the agreement is excellent (the two curves lie ontop of one another). Note that the oscillation is fully resolved, as can be seen in Fig. 3b, which shows a blow-up of v_y in the region 0 ≤ y ≤ 100 in the viscous system.

In the viscous plasma, the bulk-flow phenomenon now comes from viscous heating, as opposed to ohmic heating as in Sect. 3.2. Comparing the two systems (Fig. 3a) we see that the agreement for the ponderomotive component and bulk flow is excellent (unsurprising, since Eqs. (34) and (42) are identical for ν = η). However, the viscous damping of the acoustic component is more pronounced (with an estimated damping length of 13.5).

4. Bulk-flow phenomenon

Let us now investigate the effect that our longitudinal bulk flow has on the equilibrium density profile, starting with the ideal case. Figure 4a shows the evolution of the density perturbations at t = 1000τ_A in our ideal plasma. We can clearly see that the contributions from the ponderomotive component and acoustic wave components. We also see that the density profile is perturbed about the equilibrium but that these perturbations do not grow in time.

Figure 4b shows the density profile at t = 1000τ_A in our visco-resistive (η = ν = 0.01) system. Here, in addition to the contributions from the ponderomotive component and acoustic wave, we see there is also a third component; i.e. the bulk flow has shifted the density profile in the positive y-direction, leading to a decrease in the density profile around y = 0.

The physical interpretation of this bulk flow is as follows: visco-resistive heating increases the local temperature (with maximum around y = 0). This increase in temperature increases the thermal pressure, and the resultant pressure gradient drives the bulk flow. Note that the bulk flow is not due to the ponderomotive wave pressure force, otherwise it would be apparent in our ideal plasma (i.e. in Fig. 4a). Thus, the bulk flow is a natural consequence of driving an Alfvén wave in a nonideal plasma.

Fig. 6 a) Transverse perturbations vz for η = ν = 0.01, β = 0.1 plasma with ρ0 = 5, where red dashed lines represent the envelope vz =  ± Aexp(kIy). b) Longitudinal perturbations vy for same plasma, where black dashed line represents y = cst and blue horizontal line denotes vy = 0. In both subfigures, t = 1000τA and the vertical black dotted line represents y = vAt.

Let us now look look at the long-term evolution of the density profile. Figure 5a shows the evolution of the density profile over 0 ≤ t ≤ 100   000τ_A, where the different colours represent different times (as denoted by the colour bar). It is clear that the density profile is modified over time, and at t = 100   000τ_A the density at y = 0 is approximately 95.5% of its equilibrium value. Figure 5b shows the evolution of the temperature profile, and we see that there is a steady increase in temperature over the whole evolution, with maximum increase at y = 0.

It is clear that the visco-resistive heating has substantially changed the background density profile, and that there is a substantial flow of plasma in the positive y-direction, causing a reduction in the density since the boundary conditions do not allow for this plasma to be replaced by an in-flow. As can be seen from Fig. 5a, it appears that the rate-of-change (i.e. the decrease) in the density is decreasing (this is confirmed below in Sect. 4.1). There are two possible explanations for this decrease in the rate-of-change: either the bulk flow (in the positive y-direction) has moved so much density that an opposing pressure gradient (in the negative y-direction) has been set up, or it could be that at later times there is less visco-resistive heating in our system.

To examine the build-up of opposing pressure gradients, we can look at the evolution of ∂p / ∂y in our system (where p is total pressure), and this can be seen in Fig. 5c. We can see that there are both positive and negative perturbations, related to the oscillatory motions, but crucially we find no evidence for the build-up of an opposing pressure gradient in the negative y-direction.

Let us now investigate the alternative hypothesis; at later times there is less visco-resistive heating in our system. Recall that Eq. (41) highlighted the role of k_I in our system; k_I is relevant both to the damping of the driven Alfvén wave and to the enhanced bulk flow. Figure 5d shows the evolution of k_I in our system, and we see that the magnitude of k_I is decreasing over time. The behaviour of k_I and the density profile are directly related by Eq. (35).

This provides the explanation for the decrease in the rate-of-change of the density profile: as the simulation proceeds, the visco-resistive damping heats the plasma, leading to a decrease in density. Thus (through Eq. (35)), the magnitude of k_I decreases, thus there is less damping in the system, and thus the rate of visco-resistive heating decreases. This in turn leads to a decrease in the rate-of-change of density. There is a strong feedback effect: a decrease in density leads to a decrease in  | k_I |  and thus a decrease in the nonideal heating, and thus a decrease in the bulk flow.

Fig. 7 Long-term evolution of a) density and b) temperature, over 0 ≤ t ≤ 100   000τA, where the different colours represent different times. The colour bar denotes the different times.

Fig. 8 Evolution of ∂ρ / ∂t at y = 1 for a) ρ0 = 5 and b) ρ0 = 1 plasmas.

4.1. Dependence of equilibrium density profile

From Eq. (35) and the discussion in Sect. 4, it is clear that the larger the value of  | k_I |  in our system, the more nonideal heating occurs at that time. Thus, there should be a strong dependence upon our choice of equilibrium density profile.

Let us consider a 1D system similar to that of Sect. 4 (i.e. β = 0.1, η = ν = 0.01, boundary conditions given by Eq. (31)) where we set ρ₀ = 5. Figure 6a shows the behaviour of v_z in this system at t = 1000τ_A. As in Fig. 2a, we see that the driven Alfvén wave is damped. The magnitude of this resistive damping can be reproduced using dispersion relation (34), where we replace η by η + ν. Note that the wave is more rapidly damped than that in Fig. 2a, since in the ρ₀ = 5 system there is a larger value of  |k_I|  (|k_I|  = 6.99 × 10^-3 as opposed to  |k_I|  = 6.25 × 10^-4 in the ρ₀ = 1 system). This can also be seen from Eq. (35).

Figure 6b shows the behaviour of v_y at t = 1000τ_A. As in Fig. 3a, we observe three components to the propagation: a ponderomotive component, an acoustic (boundary-driven) component and a monotonically-increasing quadratic component, i.e. the bulk flow in the positive y-direction, as described previously. Note that the bulk flow is significantly stronger than that examined in Sect. 4. This is because in the ρ₀ = 5 system we have a larger value of  | k_I |  and thus more visco-resistive damping. Hence, there is more visco-resistive heating, which increases the local gas pressure by a greater degree, and it is this thermal-pressure gradient that drives the bulk flow.

Figures 7a and 7b show the long-term evolution of the density and temperature profile, respectively, in the ρ₀ = 5 system. As in the ρ₀ = 1 system, it is clear from Fig. 7a that the density profile has decreased substantially during the simulation, and the density at y = 0 is approximately 53.6% its initial value after t = 100   000τ_A. Figure 7a also clearly shows the bulk flow of density in the positive y-direction. Figure 7b shows that there is a monotonic increase in temperature over the whole simulation with maximum increase at y = 0.

Thus, it appears that a key variable in our system is the choice of equilibrium density profile: the greater the value of ρ₀, the greater the (initial) value of  | k_I |  and hence the greater the magnitude of visco-resistive damping (which heats the plasma, which increases the thermal pressure, which drives the bulk flow due to pressure gradients). The bulk flow decreases the value of ρ, which decreases the value of  | k_I | , which reduces the amount of visco-resistive damping, and so on (strong feedback effect). Thus, we can understand the nature of the bulk-flow phenomenon over a range of density values.

Fig. 9 a) Equilibrium density profile (black) and equilbrium Alfvén speed profile (red). b) Equilibrium temperature profile.

Figures 8a and 8b show the rate-of-change of density (∂ρ / ∂t) at y = 1 for the ρ₀ = 5 and ρ₀ = 1 systems, respectively. Let us first consider Fig. 8a. it is clear that rate-of-change is decreasing as time evolves, i.e.  | ∂ρ / ∂t |  is decreasing. Thus, as the simulation proceeds, density decreases at a slower rate, as expected from our explanation that as density decreases, there is less heating, leading to less thermal expansion. Figure 8b shows that, again, the rate-of-change of density is decreasing in the ρ₀ = 1 system but at a much slower rate than that seen in Fig. 8a. Again, this is in agreement with our explanation.

5. Two-dimensional inhomogeneous plasma

We now extend the model of Sect. 3 to include an inhomogeneous density profile. This will naturally introduce the phase mixing phenomenon into our system, in addition to the others already discussed.

We consider the following equilibrium density profile: $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{4}{\mathrm{sech}}^{\mathrm{2}}\left[{\left(\frac{\mathit{x}}{\mathit{a}}\right)}^{\mathrm{2}}\right]\mathit{,}& & \end{array}$where, in this paper we set a = 35. This equilibrium density profile (ρ₀) can be seen in Fig. 9a (black line); the red line represents the corresponding equilibrium Alfvén-speed profile. The equilibrium temperature profile (T₀) is shown in Fig. 9b, which has been chosen such that the equilibrium gas pressure is constant.

As detailed in Sect. 2.3, we solve MHD Eqs. (1) in a numerical domain of (x,y)  ∈  [−100,100]  ×  [0,10⁴]  with a uniform/stretched grid in the x- / y-direction. We drive our system with linearly-polarised, sinusoidal Alfvén waves at y = 0 (boundary condition 31). All other quantities have zero gradient boundary conditions at the y-boundaries and we utilise periodic boundary conditions at the x-boundaries.

The numerical results can be considered over two times scales: short-term evolution, over which the classical phase mixing solution dominates, and the long-term evolution over which nonlinear effects can become important.

Fig. 10 Contours of vz for 0 ≤ x ≤ 100, 0 ≤ y ≤ 500, at t = 1000τA. White dashed line denotes x = 35. The temporal evolution is shown in the movie available in the on-line edition.

5.1. Evolution of v_z: Alfvén waves

Figure 10 shows a contour plot of v_z for 0 ≤ x ≤ 100, 0 ≤ y ≤ 500 at time t = 1000τ_A. This contour plot shows the propagation of the Alfvén waves in our system, where the light and dark bands denote the peaks and troughs, respectively, of the wave motions in the z-direction. Three different behaviours are apparent: the oscillations along x = 0 correspond to the propagation of the Alfvén wave in a ρ₀ = 5 plasma and a cut along x = 0 reveals a profile identical to Fig. 6a. Similarly, the oscillations along x = 100 correspond to the propagation of an Alfvén wave in a ρ₀ = 1 plasma and a cut along x = 100 is identical to that in Fig. 2a. As expected, the propagation along x = 0 (where ρ₀ [0,y]  = 5) is damped much more rapidly than that along x = 100 (where ρ₀ [100,y]  = 1). Thus, the propagation along x = 0 and x = 100 and the (visco-resistive) damping mechanism can be fully understood from the corresponding one-dimensional results.

Let us now consider a cut along x = 35 in Fig. 10, i.e. v_z(35,y); this can be seen in Fig. 11. Here, we see that the wave is rapidly damped (more rapidly than standard visco-resistive damping).

The classical phase mixing solution of Heyvaerts & Priest (1983) is given by: $\begin{array}{ccc}{\mathit{v}}_{\mathit{z}}\mathrm{=}\mathrm{\pm }\mathit{A}\exp \left({\mathit{k}}_{\mathrm{I}}\mathit{y}\right)\exp \left[\mathrm{-}\frac{\mathrm{1}}{\mathrm{6}}{\left(\frac{{\mathit{k}}_{\mathrm{\parallel }}\mathit{y}}{{\mathit{R}}^{\mathrm{1}\mathit{/}\mathrm{3}}}\right)}^{\mathrm{3}}\right]\mathit{,}& & \end{array}$(43)where k _∥  = ω / v_A(x) and $\mathit{R}\mathrm{=}\mathit{\omega }\mathit{/}\mathrm{(}\mathit{\eta }\mathrm{+}\mathit{\nu }\mathrm{)}{{}^{\mathrm{\right[}}\frac{\mathrm{d}}{\mathrm{d}\mathit{x}}\left(\ln {\mathit{k}}_{\mathrm{\parallel }}\right)^{\mathrm{\left]}}}^{\mathrm{2}}$. This solution is overplotted in Fig. 11 (red dashed line) and the agreement is excellent. Thus, at early times, the damping of v_z is well understood as that corresponding to the classic phase mixing mechanism.

A movie of the evolution of v_z over the whole simulation (up to t = 100   000τ_A) is associated with Fig. 10. It is clear that the behaviour of v_z changes very little over the whole simulation, and that Alfvén waves are continuously driven into the numerical domain; there is no evidence of the system choking itself off. Note that the small changes in v_z are dictated by the changes in the density profile (see Sect. 5.4 below).

Fig. 11 Plot of vz(35,y). Heyvaerts and Priest solution overplotted.

Fig. 12 Shaded surface of vy(x,y) at t = 1000τA. The temporal evolution is shown in the movie available in the on-line edition.

5.2. Evolution of v_y: slow magnetoacoustic waves

Figure 12 shows a shaded-surface of v_y(x,y) at t = 1000τ_A. As in the analysis of v_z, a cut along x = 0 reveals behaviour identical to that of Fig. 6b, and a cut along x = 100 reveals behaviour identical to that of v_y in a one-dimensional plasma with β = 0.1, η = ν = 0.01 at t = 1000τ_A, i.e. similar to Fig. 3a. Thus, as in our one-dimensional analyses, we see that our simulation contains three types of longitudinal motions: slow magnetoacoustic waves (previously called acoustic waves in Sect. 3), ponderomotive wave components (propagating at the Alfvén speed) and a bulk flow in the positive y-direction. As expected, all amplitudes are of the order A².

A movie of the evolution of v_y over the whole simulation (up to t = 100   000τ_A) is associated with Fig. 12. It is clear that the growth of v_y saturates, and that the behaviour of v_y is dictated by that of v_z and the heating profile. Thus, the complete evolution of v_y is of secondary importance in understanding the evolution of the density and temperature profiles.

Fig. 13 a) Contour of vx(x,y) at t = 100τA. b) Shaded surface of vx(x,y) at t = 1000τA. The temporal evolution of the right panel is shown in the movie available in the on-line edition.

5.3. Evolution of v_x: fast magnetoacoustic waves

In addition to Alfvén waves and slow magnetoacoustic waves, there is now a third type of MHD wave present in our system: fast magnetoacoustic waves, which can be seen in Fig. 13. These fast waves are nonlinearly generated by transverse gradients in the Alfvén-speed profile, in agreement with the analytical work of Nakariakov et al. (1997). The fast waves are permanently generated by phase mixing, and are refracted towards regions of lower Alfvén speed. This can be clearly seen in Fig. 13a, which shows a contour plot of v_x behaviour at t = 100τ_A. These oscillations have a maximum around x = 35 and are refracted in towards x = 0 (i.e. towards regions of lower Alfvén speed, also see Fig. 9a). These wave motions are generated with a frequency twice that of the driving frequency, which is in agreement with previous analytical predictions (Nakariakov et al. 1997, 1998).

A shaded-surface of the behaviour of v_x at t = 1000τ_A can be seen in Fig. 13b. Here, we can see that the generated fast waves have propagated across the magnetic fieldlines and have set-up an interference pattern. The interference pattern also results from waves generated along x = 35 overlapping with those generated along x =  − 35, and there is also overlap due to the periodic boundary conditions. However, as first noted by Botha et al. (2000), these fast waves grow initially and then saturate. Botha et al. explain this saturation in terms of simple wave interference: physically, the saturation of fast waves is due to destructive interference from incoherent sources.

A movie of the evolution of v_x over the whole simulation (up to t = 100   000τ_A) is associated with Fig. 13. As with the evolution of v_y, it is clear that v_x quickly saturates at amplitudes of order A² and, hence, it is only of secondary importance in understanding the evolution of the density and temperature profiles.

Fig. 14 a) Shaded surface of ρ(x,y) at t = 0 (i.e. equilibrium density profile). b) Shaded surface of ρ(x,y) at t = 100   000τA.

5.4. Evolution of density profile

Let us now look at the evolution of the density profile in our system. Figure 14 compares shaded-surfaces of the initial density profile (Fig. 14a) at t = 0 and the final density profile (Fig. 14b) at t = 100   000τ_A at the end of our simulation. Note that we have presented  − 100 ≤ x ≤ 0 here, for clarity, and recall that the simulation is symmetric about x = 0.

It is clear that the density profile has changed substantially during the simulation. Figure 15 shows cuts along x = 0, 35 and 100 at the end of the simulation. The cuts along x = 0 and x = 100 are in excellent agreement with those predicted by the one-dimensional model (i.e. compare to Figs. 5a and 7a at t = 100   000τ_A). Thus, the nature of the density change at these locations can be understood in terms of the bulk-flow phenomenon described in Sect. 4.

The nature of the density change along x = 35 is however completely different. Here, the bulk flow has been suppressed (very little density change can be seen after y ≈ 225) and instead the localised decrease in density is due to the strong, localised heating from phase mixing (see below).

5.5. Evolution of temperature profile

Figure 16 shows contour plots of T − T₀ at times (a) t = 1000τ_A, (b) t = 10   000τ_A and (c) t = 100   000τ_A. In Fig. 16a, we see that there is strong, localised heating (maximum occurs at x = 38.6, y = 70.4) which results from the enhanced visco-resistive dissipation associated with phase mixing, producing a distinctive teardrop shape. Independently, the heating in the bottom left-corner is the visco-resistive heating associated Alfvén-wave damping along ρ₀ = 5.

At a later time (Fig. 16b, t = 10   000τ_A) the overall temperature profile appears qualitatively similar to that in Fig. 16a, i.e. phase mixing is still the dominant heating mechanism. However, the maximum of T − T₀ (located at x = 39.6, y = 80.5) is now approximately ten times greater (note that the colour bar has changed).

Figure 16c shows the temperature profile at the end of our simulation (t = 100   000τ_A). Again, qualitatively the contours appears similar to those in (a) and (b). However, the magnitude has again increased substantially (the colour bar has again changed) and the maximum of T − T₀ is now located at x = 37.7, y = 108.8. Again, phase mixing is still the dominant heating mechanism.

Figure 16d shows the location of the maximum of T − T₀ for 2000τ_A ≤ t ≤ 100   000τ_A (note the different length scales on the axes). It is clear that the location of the maximum of T − T₀ changes during our simulation, and that this drifting is in the direction of decreasing Alfvén-speed profile. This result addresses one of the fundamental question asked in Sect. 1, i.e. we find that drifting of the heating layer in nonideal, nonlinear simulations of phase mixing is possible. However, we have also shown that such drifting is very small and takes a very long time, at least for the parameters considered in this paper.

Finally, note that before t = 2000τ, the maximum of T − T₀ sometimes occurs due to the bulk-flow heating, rather than due to phase mixing (phase mixing takes a finite amount of time to become the dominant heating mechanism) and thus, for clarity, we have not included these points in Fig. 16d. However, for comparison we have included the location of the maximum of T − T₀ at t = 1000τ (i.e. x = 38.6, y = 70.4) and this is denoted by a star.

Fig. 15 Density profiles of ρ(x,y) along the lines x = 0 (blue), x = 35 (black) and x = 100 (red), at t = 100   000τA.

Fig. 16 Contour of T − T0 at a) t = 1000τA, b) t = 10   000τA and c) t = 100   000τA. Dashed white line denotes x = 35. d) Plots (diamonds) the location of the maximum of T − T0 over 2000τA ≤ t ≤ 100   000τA. Also plotted (single star) is the location at t = 1000τA. In subfigures a) − c), each colour bar represents magnitude, whereas in (d) the corresponding colour bar donotes time.

6. Conclusions

We have investigated the nonlinear, nonideal behaviour of Alfvén wave propagation and phase mixing within an inhomogeneous environment, over long timescales. The governing MHD equations have been solved in 1D and 2D environments using both analytical techniques and numerical simulations.

In an ideal, one-dimensional study (no density inhomogeneity, ∂ / ∂x = 0) we find that by driving a linear Alfvén wave (in v_z) into our numerical domain, we also generate two types of longitudinal wave: boundary-driven acoustic waves (propagating at speed c_s) and a nonlinear perturbation (propagating at v_A) which is driven by the ponderomotive force. Both these motions have a frequency twice that of the driven Alfvén wave, and an exact mathematical solution for both wave types was derived. The acoustic wave is degenerate under the β = 0 approximation (see Appendix A) but the ponderomotive wave is always present in a nonlinear system.

We find that the addition of resistive and viscous effects naturally leads to visco-resistive damping (through the dispersion relation) and thus to visco-resistive heating, which in turn leads to the introduction of a new phenomenon: a bulk flow in the positive y-direction. Physically, this bulk flow is a direct response to the visco-resistive heating in the system: the heating increases the local temperature which increases the local thermal pressure. The resulting presure gradient drives the bulk flow. The bulk flow is purely in the positive y-direction because the boundary conditions are held fixed (v_y = 0). As a result of this out-flow, a significant reduction in density occurs since the boundary conditions do not allow for the plasma to be replaced by an in-flow. Note that the bulk flow is not due to the ponderomotive force, otherwise such a flow would be apparent in the ideal, nonlinear system.

We find that the magnitude of the visco-resistive heating is strongly dependent upon k_I and thus on our choice of ρ₀; the greater the value of ρ₀, the more visco-resistive heating that occurs. More visco-resistive heating leads to a stronger bulk flow in the longitudinal direction.

We also investigated the effect of driving an Alfvén wave in a nonideal plasma over very long timescales. We find that the equilibrium density profile changes substantially over the duration of our simulation (end of the simulation at t = 100   000τ_A) specifically due to this bulk-flow phenomenon (which is itself directly due to the visco-resistive heating). For a ρ₀ = 1 plasma, we find a decrease in density of about 4.6% at y = 0 but for a ρ₀ = 5 plasma we find a decrease of about 46.4%. We also find that the rate-of-change of density at y = 0 is decreasing as the simulation proceeds. This is because as the density decreases, the value of  |k_I|  also decreases. This smaller value of  |k_I|  reduces the amount of visco-resistive damping, which slows the rate-of-change of density, and so on. Thus, there is a strong feedback effect in our system, and we can explain the nature of the bulk flow and change of density profile over a range of density values. We conclude that the bulk-flow phenomenon is a natural consequence of driving an Alfvén wave in a nonideal plasma.

We then extended our analysis to include an inhomogeneous density profile in a two-dimensional system and, as before, we drive our system with sinusodial, linear Alfvén waves. As expected from Heyvaerts & Priest (1983), our simulation displays classic phase mixing. A cut showing v_z as a function of height along x = 35 provided an excellent match with the analytical predictions of Heyvaerts & Priest (1983). Cuts along x = 0 (ρ₀ = 5) and x = 100 (ρ₀ = 1) were identical to those of the corresponding visco-resistive damped Alfvén waves seen in our one-dimensional analysis.

Our one-dimensional analysis also explained the longitudinal motions (v_y) in our system and we again observed the boundary-driven acoustic waves (interpreted as slow magnetoacoustic waves in 2D, with speed c_s), our ponderomotive waves (propagating at v_A) and the bulk-flow phenomenon.

In addition to Alfvén waves and slow magnetoacoustic waves, our 2D system also contains fast magnetoacoustic waves (seen primarily in v_x). These fast waves are continuously generated by Alfvén-wave phase mixing at a frequency twice that of the driven Alfvén wave, and propagate across magnetic fieldlines and away from the phase mixing layer (as predicted by Nakariakov et al. 1997, 1998). These waves are also refracted towards regions of lower Alfvén speed.

Since there is a permanent leakage of energy away from the phase mixing layer, it is possible that these fast waves can cause indirect heating of the plasma as they propagate away and dissipate far from the phase mixing layer itself, thus spreading the heating across the domain (Nakariakov et al. 1997). However, due to our choice of amplitude (A = 0.01) we find that only a small fraction of the Alfvén-wave energy is converted into fast waves and, thus, only a small amount of indirect heating occurs.

Over long timescales, we find that both the slow and fast waves saturate, and always remain at amplitudes of the order A². Hence, these waves only play a secondary role in the heating of the plasma.

We find that at the end of our simulation, the density profile has changed substantially. This is due to two effects: firstly, the bulk-flow phenomenon is naturally present in our system, with the density changes as predicted by our one-dimensional analysis. Secondly, there is a substantial decrease in density at the locations of phase mixing heating.

The change in the background density profile does not affect the propagation of the Alfvén wave (amplitude of the order A) to a great degree, and the behaviour of v_x (fast waves) and v_y (slow waves) do not have a strong feedback effect on the system (both amplitudes of the order A²) and only play secondary roles in the system (i.e. A = 0.01, so feedback effect is only of the order of 1%). However, these conclusions may be different for Alfvén waves driven with a much larger amplitudes.

We have also investigated how the temperature profile changes during our simulation. We find a substantial increase in localised temperature, generated by the phase mixing mechanism, over the duration of our simulation, and find that the maximum temperature increases by a factor of about a hundred during the simulation. We also find that the location of the maximum temperature drifts during our simulation, thus providing an answer to one of the key questions asked in Sect. 1, i.e. we find that drifting of the heating layer is possible in a nonlinear phase mixing scenario.

However, we also find that this drifting is very small and occurs over a very long timescale (our simulation runs over 100   000τ_A), at least for the parameters we have considered in this paper. For typical coronal parameters, τ_A = 60 seconds and so we are considering timescales of approximately 69.4 days, which is too long to be physically important in the corona (i.e. other coronal processes act over much shorter timescales). However, it may be possible to increase the speed and magnitude of this drifting by significantly increasing the magnitude of heating in our system, i.e. attempt to increase the amount of enhanced dissipation from phase mixing. Thus, future studies could consider increasing the amplitude and frequency of the driven Alfvén wave, or increasing the steepness of the gradient in our background density inhomogeneity.

Online material

Movie 1

Movie 2

Movie 3

Appendix A: Alfvén wave propagation in a one-dimensional, ideal, β = 0 plasma

Fig. A.1 a) Transverse perturbations vz for ideal (η = ν = 0) β = 0 plasma. b) Longitudinal perturbations vy for ideal β = 0 plasma. In both subfigures, t = 1000τA and dotted line represents y = vAt.

Consider a cold (β = 0), ideal (η = ν = 0) plasma in our one-dimensional system (Sect. 3). Driving this system with boundary condition (31) generates an Alfvén wave at y = 0; this can be seen in Fig. A.1a. The Alfvén wave (v_z) propagates in the direction of increasing/positive y with amplitude A = 0.01 and there is no damping. The time of the snapshot can be read directly from the y-axis using the relation t = y / v_A. The ramp-up over the first four wavelengths is clearly visible.

Figure A.1b shows the longitudinal motions (v_y) in the system. Here, v_y is driven by the nonlinear terms of Eq. (29) (~B_z∂B_z / ∂y) and again there is no damping. In this paper, we call this second-order nonlinear effect the ponderomotive effect, and thus these longitudinal motions, which propagate at the Alfvén speed, are driven by the ponderomotive force (Alfvén wave-pressure gradients). At early times, these longitudinal motions are governed by a simplified version of Eq. (29): $\begin{array}{ccc}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{v}}_{\mathit{y}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{B}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\right)\mathrm{·}& & \end{array}$(A.1)Assuming boundary conditions (31) in our ideal, cold plasma and that v_x, v_y are initially zero, Eq. (A.1) has an analytical solution of the form: $\begin{array}{ccc}{\mathit{v}}_{\mathit{y}}\mathrm{=}\frac{{\mathit{A}}^{\mathrm{2}}}{\mathrm{4}{\mathit{v}}_{\mathrm{A}}}\left\{\mathrm{1}\mathrm{-}\cos \left(\mathrm{2}\mathit{\omega }\mathrm{[}\mathit{t}\mathrm{-}\mathit{y}\mathit{/}{\mathit{v}}_{\mathrm{A}}\mathrm{]}\right)\right\}\mathit{,}& & \end{array}$(A.2)which is valid for 0 < y < v_At. Note that these perturbations have an amplitude of and are always positive for an Alfvén wave propagating in the positive y-direction (due to the constant of integration, and thus due to the boundary conditions). Finally, note that the frequency of these longitudinal motions is twice that of the driving frequency, and that the motions do not grow with time. Interestingly, the ramp-up to maximum amplitude now occurs over eight wavelengths, i.e. twice that of the driven wave.

Appendix B: Derivation of Eqs. (38)–(41)

Equation (29) governs the driven, longitudinal motions in the system, and in 1D (∂ / ∂x = 0) has the form: $\begin{array}{ccc}& & \left(\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right){\mathit{v}}_{\mathit{y}}\mathrm{=}\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{B}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\right)\mathrm{-}\frac{\left(\mathit{\gamma }\mathrm{-}\mathrm{1}\right)\mathit{\eta }}{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial y}}\left[{\left(\frac{\mathit{\partial }{\mathit{B}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\right]\mathrm{·}\end{array}$(B.1)

Let us assume that our linear Alfvén wave can be represented as:

$\begin{array}{ccc}{\mathit{v}}_{\mathit{z}}\mathrm{=}\frac{\mathit{A}}{\mathrm{2}\mathit{i}}{\mathrm{e}}^{\mathrm{i}\left(\mathit{\omega t}\mathrm{-}\mathit{ky}\right)}\mathrm{-}\frac{\mathit{A}}{\mathrm{2}\mathit{i}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}^{\mathrm{\ast }}\mathit{y}\right)}& & \end{array}$

where A is the amplitude of our wave, ω is the driving frequency, k = k_R + ik_I is our wavenumber, with complex conjugate k ^∗  = k_R − ik_I. Thus, from Eq. (7) (with N₃ = 0) we have:

$\begin{array}{ccc}{\mathit{b}}_{\mathit{z}}\mathrm{=}\mathrm{-}\frac{\mathit{\omega }{\mathit{B}}_{\mathrm{0}}\mathit{A}}{\mathrm{2}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}}\left[\frac{\mathrm{1}}{\mathrm{i}\mathit{k}}{\mathrm{e}}^{\mathrm{i}\left(\mathit{\omega t}\mathrm{-}\mathit{ky}\right)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{i}{\mathit{k}}^{\mathrm{\ast }}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}^{\mathrm{\ast }}\mathit{y}\right)}\right]\mathit{.}& & \end{array}$(B.2)

The right-hand-side of Eq. (B.1) has two contributions, which can be calculated using our expression for b_z. The first term has the form:

$\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}& \mathrm{=}& \frac{\mathit{\omega }{\mathit{B}}_{\mathrm{0}}\mathit{A}}{\mathrm{2}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}}\mathrm{[}{\mathrm{e}}^{\mathrm{i}\left(\mathit{\omega t}\mathrm{-}\mathit{ky}\right)}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}^{\mathrm{\ast }}\mathit{y}\right)}\mathrm{]}\\ \mathrm{\Rightarrow }\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{b}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}\right)& \mathrm{=}& \end{array}$(B.3)

and the second has the form:

$\begin{array}{ccc}& & {\left(\frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{A}}^{\mathrm{2}}}{\mathrm{4}{\mathit{v}}_{\mathrm{A}}^{\mathrm{4}}}\mathrm{[}{\mathrm{e}}^{\mathrm{2}\mathrm{i}\left(\mathit{\omega t}\mathrm{-}\mathit{ky}\right)}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}^{\mathrm{\ast }}\mathit{y}\right)}\mathrm{+}\mathrm{2}{\mathrm{e}}^{\mathrm{i}\left({\mathit{k}}^{\mathrm{\ast }}\mathrm{-}\mathit{k}\right)\mathit{y}}\mathrm{]}\\ & & \mathrm{\Rightarrow }\frac{\mathit{\partial }}{\mathit{\partial y}}\left[{\left(\frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\right]\mathrm{=}\mathrm{i}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{A}}^{\mathrm{2}}}{\mathrm{2}{\mathit{v}}_{\mathrm{A}}^{\mathrm{4}}}\mathrm{[}\mathrm{-}\mathit{k}{\mathrm{e}}^{\mathrm{2}\mathrm{i}\left(\mathit{\omega t}\mathrm{-}\mathit{ky}\right)}\mathrm{+}{\mathit{k}}^{\mathrm{\ast }}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}^{\mathrm{\ast }}\mathit{y}\right)}\\ & & \mathrm{+}\left({\mathit{k}}^{\mathrm{\ast }}\mathrm{-}\mathit{k}\right){\mathrm{e}}^{\mathrm{i}\left({\mathit{k}}^{\mathrm{\ast }}\mathrm{-}\mathit{k}\right)\mathit{y}}{\mathit{.}}^{\mathrm{]}}\end{array}$(B.4)

The form of Eqs. (B.2) − (B.4) lead to Eqs. (38)−(41).

Equations (B.2) − (B.4) can also be derived explicitly from our choice of v_z, i.e.:

$\begin{array}{ccc}{\mathit{v}}_{\mathit{z}}\mathrm{=}\mathit{A}\sin \left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right){\mathrm{e}}^{{\mathit{k}}_{\mathrm{I}}\mathit{y}}& & \end{array}$where as before k = k_R + ik_I. Thus (from Eq. (7)) b_z has the form: $\begin{array}{ccc}{\mathit{b}}_{\mathit{z}}& \mathrm{=}& \mathrm{-}\frac{\mathit{\omega }{\mathit{B}}_{\mathrm{0}}\mathit{A}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{k}}_{\mathrm{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{k}}_{\mathrm{I}}^{\mathrm{2}}}^{\mathrm{)}}}[{\mathit{k}}_{\mathrm{R}}\sin \left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)\\ & & \end{array}$(B.5)

Now we can explicitly determine the two contributions:

$\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}\mathrm{=}\frac{\mathit{\omega }{\mathit{B}}_{\mathrm{0}}\mathit{A}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}}\cos \left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right){\mathrm{e}}^{{\mathit{k}}_{\mathrm{I}}\mathit{y}}\\ & & \mathrm{\Rightarrow }\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{b}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}\right)\mathrm{=}\mathrm{-}\frac{{\mathit{\omega }}^{\mathrm{3}}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{A}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{4}}\mathrm{(}{\mathit{k}}_{\mathrm{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{k}}_{\mathrm{I}}^{\mathrm{2}}}^{\mathrm{)}}}[{\mathit{k}}_{\mathrm{R}}\cos \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)\\ & & \mathrm{+}{\mathit{k}}_{\mathrm{I}}\sin \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)]{\mathrm{e}}^{\mathrm{2}{\mathit{k}}_{\mathrm{I}}\mathit{y}}\end{array}$(B.6)

and the second has the form:

$\begin{array}{ccc}& & {\left(\frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{A}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{4}}}{\cos }^{\mathrm{2}}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right){\mathrm{e}}^{\mathrm{2}{\mathit{k}}_{\mathrm{I}}\mathit{y}}\\ & & \mathrm{\Rightarrow }\frac{\mathit{\partial }}{\mathit{\partial y}}\left[{\left(\frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial y}}\right)}^{\mathrm{2}}\right]\mathrm{=}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{A}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{A}}^{\mathrm{4}}}[{\mathit{k}}_{\mathrm{R}}\sin \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)\\ & & \mathrm{+}{\mathit{k}}_{\mathrm{I}}\cos \mathrm{2}\left(\mathit{\omega t}\mathrm{-}{\mathit{k}}_{\mathrm{R}}\mathit{y}\right)\mathrm{+}{\mathit{k}}_{\mathrm{I}}]{\mathrm{e}}^{\mathrm{2}{\mathit{k}}_{\mathrm{I}}\mathit{y}}\mathit{.}\end{array}$(B.7)

Equations (B.5) − (B.7) give equivalent solutions to Eqs. (B.2) − (B.4), and thus also lead to Eqs. (38) − (41).