Issue 
A&A
Volume 527, March 2011



Article Number  A149  
Number of page(s)  16  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201015552  
Published online  15 February 2011 
Phase mixing of nonlinear viscoresistive Alfvén waves
School of Mathematics and Statistics, University of St Andrews, KY16 9SS, UK
email: james.a.mclaughlin@northumbria.ac.uk
Received:
9
August
2010
Accepted:
11
January
2011
Aims. We investigate the behaviour of nonlinear, nonideal Alfvén wave propagation within an inhomogeneous magnetic environment.
Methods. The governing MHD equations are solved in 1D and 2D using both analytical techniques and numerical simulations.
Results. We find clear evidence for the ponderomotive effect and viscoresistive heating. The ponderomotive effect generates a longitudinal component to the transverse Alfvén wave, with a frequency twice that of the driving frequency. Analytical work shows the addition of resistive heating. This leads to a substantial increase in the local temperature and thus gas pressure of the plasma, resulting in material being pushed along the magnetic field. In 2D, our system exhibits phase mixing and we observe an evolution in the location of the maximum heating, i.e. we find a drifting of the heating layer.
Conclusions. Considering Alfvén wave propagation in 2D with an inhomogeneous density gradient, we find that the equilibrium density profile is significantly modified by both the flow of density due to viscoresistive heating and the nonlinear response to the localised heating through phase mixing.
Key words: magnetohydrodynamics (MHD) / magnetic fields / waves / shock waves / Sun: corona / Sun: oscillations
Appendices and movies are only available in electronic form at http://www.aanda.org.
© 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énwave 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énwave 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 weaklynonlinear, β = 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énwave pulses and found that this results in a slower powerlaw damping (as opposed to the standard ~exp(−t^{3}) for harmonic Alfvén waves, i.e. Heyvaerts & Priest 1983). Hence, Alfvénicpulse 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énwave 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 coauthors repeated these numerical experiments for both a weakly and stronglynonlinear 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 semianalytical 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, selfsimilar solutions of Alfvénwave 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 coauthors who considered the interaction of an impulsivelygenerated, weaklynonlinear MHD pulse with a onedimensional density inhomogeniety, considered in the threedimensional regime (i.e. without an ignorable coordinate) in both an ideal (Tsiklauri & Nakariakov 2002) and resistive (Tsiklauri et al. 2003) plasma. Tsiklauri and coauthors 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 ioncyclotron 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 wavepressure force. The ponderomotive effect has been considered in a solar context initially by Hollweg (1971) and later by Verwichte and coauthors (Verwichte 1999; Verwichte et al. 1999). Hollweg (1971) considered linearlypolarised 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 magneticfield 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 weaklynonlinear, Alfvén wave in a β = 0 homogeneous plasma. These authors showed that the an initiallyexcited, gaussianpulse 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énwave 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 viscoresistive 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 bulkflow 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 longterm 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: (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^{7}Hm^{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 where ϵ is the rateofstrain tensor, given by: The viscous heating term is given by 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 nonisotropic tensor. The mathematical details of this nonisotropic 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 nonisotropic 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 nondimensionalise all variables. Let v = v_{0}v^{ ∗ }, B = BB^{ ∗ }, x = Lx^{ ∗ }, y = Ly^{ ∗ }, ρ = ρ_{0}ρ^{ ∗ }, p = p_{0}p^{ ∗ }, , t = τ_{A}t^{ ∗ }, η = η_{0} and ν = ν_{0}, where ∗ denotes a dimensionless quantity and v_{0}, B, L, ρ_{0}, p_{0}, τ_{A}, η_{0} and ν_{0} are constants with the dimensions of the variable they are scaling. We then set and v_{0} = L / τ_{A} (this sets v_{0} as a constant equilibrium Alfvén speed). We also set , where R_{m} is the magnetic Reynolds number. This process nondimensionalises Eqs. (1) and under these scalings, t^{ ∗ } = 1 (for example) refers to t = τ_{A} = L / v_{0}; 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 nondimensionalised is understood.
2.1. Basic equilibrium
We consider Eqs. (1) in Cartesian coordinates and assume there are no variations in the zdirection (∂ / ∂z = 0). We consider an inhomogeneous plasma in the xdirection embedded in a uniform magnetic field in the ydirection, i.e.: (2)with finite amplitude perturbations of the form: (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): where the lefthandsides involve only linear terms and N_{1} − N_{8} denote the nonlinear terms. The nonlinear terms are: the resistive terms (where R_{1} − R_{3} are linear, R_{4} is nonlinear) are: and the viscous terms (V_{1} − V_{3} are linear, V_{4} is nonlinear) are 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: where is the unperturbed Alfvén speed and is the unperturbed sound speed. In Eqs. (28) − (30) the nonlinear and nonideal terms can be found on the righthandside. 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 (righthandside), 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_{2} (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_{3} (i.e. ~B_{z}∂B_{z} / ∂y) corresponds to the Alfvénwave 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 v_{y} for β = 0.1 plasma with η = ν = 0. Dotted and dashed lines represent y = v_{A}t and y = c_{s}t respectively. The blue line denotes the maximum amplitude: . b) Blowup of region 100 ≤ y ≤ 200 for v_{y}. In both subfigures, t = 1000τ_{A}. 
2.3. Numerical Setup
We solve MHD Eqs. (1) using the LARE2D Lagrangianremap code (Arber et al. 2001).
We consider a numerical domain with − 100 ≤ x ≤ 100, 0 ≤ y ≤ 10 000 with uniformgrid spacing in the xdirection and a stretched grid in the ydirection. The stretched grid places the majority of the ydirection 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 linearlypolarised, sinusoidal Alfvén waves. This harmonic wave train is driven at the y = 0 boundary such that: (31)where A is the amplitude and ω is the frequency. All other quantities have zero gradient boundary conditions at the yboundaries, and we utilise periodic boundary conditions at the xboundaries. 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^{2} (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. Onedimensional system
Let us start by considering the nonlinear and nonideal aspects of the Alfvén wave. In this section, we set ρ_{0} = constant, i.e. we remove the density inhomogeneity. This effectively reduces Eqs. (4) − (30) to a onedimensional (1D) system where variations in the xdirection 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: (32)where the arbitrary function ℱ can include a form for the rampup of the driver as well as the periodic driver. To illustrate the terms, we ignore the rampup and consider the sinusodial driver; the solutions are: and 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_{s}t ≈ 280 and y = v_{A}t = 1000. These longitudinal motions have a maximum amplitude of (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_{s}t. 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 boundarydriven slow wave. Note that the longitudinal motions for 0 ≤ y ≤ c_{s}t 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 blowup 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_{s}t at t = 1000τ_{A} and so, after y = c_{s}t ≈ 280, only one wave is present.
The analytical solution for these nonlinearlygenerated longitudinal motions is found by substituting Eq. (32) into the weaklynonlinear form of Eq. (29), namely: Thus, the analytical solution for v_{y} is: (33)which is valid only for early times (i.e. until nonlinear response becomes nonnegligible). Note that these longitudinal motions were first reported in Botha et al. (2000).
Note that for 0 < y < c_{s}t, 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: Thus, for the typical numerical resolutions considered in this paper (δy ~ 0.1), both these perturbations are well resolved.
Fig. 2 a) Transverse perturbations v_{z} for η = 0.01, β = 0.1 plasma (ν = 0), where dashed lines represent the envelope v_{z} = ± Aexp(k_{I}y). b) Longitudinal perturbations v_{y} for same plasma, where dashed line represents y = c_{s}t and blue horizontal line denotes v_{y} = 0. In both subfigures, t = 1000τ_{A} and dotted line represents y = v_{A}t. 
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: (34)where ω = 0.25, , η = 0.01 and k = k_{R} + ik_{I}. In Fig. 2a, the dashed lines represent the envelope v_{z} = ± Aexp(k_{I}y), where k_{I} < 0. Assuming (or equivalently that the magnetic Reynolds number is large) we can estimate that: (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 ydirection with v_{y} > 0 and that has a maximum around y = c_{s}t. 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: (36)The righthandside of this equation has two contributions: the first term depends upon the rateofchange of the magnetic pressure gradient ) of the driven wave and the second depends upon pressure gradients created by the ohmic heating term (i.e. ) 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: (37)Thus, from Eqs. (7) and (10), B_{z} has the form: (38)The derivation of Eq. (38) is given in Appendix B.
Fig. 3 a) Longitudinal perturbations v_{y} 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) Blowup of v_{y} over region 0 ≤ y ≤ 100 for ν = 0.01, η = 0 plasma only, where crosses indicate grid point resolution. Black line indicates v_{y} = 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 = c_{s}t (y = v_{A}t) and the blue line denotes ρ = ρ_{0}. 
Hence, the righthandside of Eq. (36) can be simplified such that: where the terms representing the Ponderomotive force, , the resistive heating force, , and the nonoscillatory resistive bulk flow force, , are given by: which is valid for c_{s}t < y < v_{A}t. 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, nonoscillatory and it is this term that will create the bulk flow in the longitudinal direction.
Fig. 5 Longterm evolution of a) density, b) temperature, c) pressure gradient (∂p / ∂y), and d) k_{I}, 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: (42)Note the similarity to Eq. (34).
Figure 3a shows a comparison of v_{y} for the purelyresistive (blues, Sect. 3.2, Fig. 2b) and purelyviscous (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 ydirection, and the acoustic perturbations (from 0 ≤ y ≤ 280). However, note that the acoustic perturbations now take a different form to those in the purelyresistive 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 blowup of v_{y} in the region 0 ≤ y ≤ 100 in the viscous system.
In the viscous plasma, the bulkflow 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. Bulkflow 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 viscoresistive (η = ν = 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 ydirection, leading to a decrease in the density profile around y = 0.
The physical interpretation of this bulk flow is as follows: viscoresistive 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 v_{z} for η = ν = 0.01, β = 0.1 plasma with ρ_{0} = 5, where red dashed lines represent the envelope v_{z} = ± Aexp(k_{I}y). b) Longitudinal perturbations v_{y} for same plasma, where black dashed line represents y = c_{s}t and blue horizontal line denotes v_{y} = 0. In both subfigures, t = 1000τ_{A} and the vertical black dotted line represents y = v_{A}t. 
Let us now look look at the longterm 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 viscoresistive heating has substantially changed the background density profile, and that there is a substantial flow of plasma in the positive ydirection, causing a reduction in the density since the boundary conditions do not allow for this plasma to be replaced by an inflow. As can be seen from Fig. 5a, it appears that the rateofchange (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 rateofchange: either the bulk flow (in the positive ydirection) has moved so much density that an opposing pressure gradient (in the negative ydirection) has been set up, or it could be that at later times there is less viscoresistive heating in our system.
To examine the buildup 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 buildup of an opposing pressure gradient in the negative ydirection.
Let us now investigate the alternative hypothesis; at later times there is less viscoresistive 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 rateofchange of the density profile: as the simulation proceeds, the viscoresistive 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 viscoresistive heating decreases. This in turn leads to a decrease in the rateofchange 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 Longterm 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 ρ_{0} = 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 ρ_{0} = 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 ρ_{0} = 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 (boundarydriven) component and a monotonicallyincreasing quadratic component, i.e. the bulk flow in the positive ydirection, as described previously. Note that the bulk flow is significantly stronger than that examined in Sect. 4. This is because in the ρ_{0} = 5 system we have a larger value of  k_{I}  and thus more viscoresistive damping. Hence, there is more viscoresistive heating, which increases the local gas pressure by a greater degree, and it is this thermalpressure gradient that drives the bulk flow.
Figures 7a and 7b show the longterm evolution of the density and temperature profile, respectively, in the ρ_{0} = 5 system. As in the ρ_{0} = 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 ydirection. 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 ρ_{0}, the greater the (initial) value of  k_{I}  and hence the greater the magnitude of viscoresistive 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 viscoresistive damping, and so on (strong feedback effect). Thus, we can understand the nature of the bulkflow 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 rateofchange of density (∂ρ / ∂t) at y = 1 for the ρ_{0} = 5 and ρ_{0} = 1 systems, respectively. Let us first consider Fig. 8a. it is clear that rateofchange 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 rateofchange of density is decreasing in the ρ_{0} = 1 system but at a much slower rate than that seen in Fig. 8a. Again, this is in agreement with our explanation.
5. Twodimensional 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: where, in this paper we set a = 35. This equilibrium density profile (ρ_{0}) can be seen in Fig. 9a (black line); the red line represents the corresponding equilibrium Alfvénspeed profile. The equilibrium temperature profile (T_{0}) 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^{4}] with a uniform/stretched grid in the x / ydirection. We drive our system with linearlypolarised, sinusoidal Alfvén waves at y = 0 (boundary condition 31). All other quantities have zero gradient boundary conditions at the yboundaries and we utilise periodic boundary conditions at the xboundaries.
The numerical results can be considered over two times scales: shortterm evolution, over which the classical phase mixing solution dominates, and the longterm evolution over which nonlinear effects can become important.
Fig. 10 Contours of v_{z} 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 online 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 zdirection. Three different behaviours are apparent: the oscillations along x = 0 correspond to the propagation of the Alfvén wave in a ρ_{0} = 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 ρ_{0} = 1 plasma and a cut along x = 100 is identical to that in Fig. 2a. As expected, the propagation along x = 0 (where ρ_{0} [0,y] = 5) is damped much more rapidly than that along x = 100 (where ρ_{0} [100,y] = 1). Thus, the propagation along x = 0 and x = 100 and the (viscoresistive) damping mechanism can be fully understood from the corresponding onedimensional 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 viscoresistive damping).
The classical phase mixing solution of Heyvaerts & Priest (1983) is given by: (43)where k_{ ∥ } = ω / v_{A}(x) and . 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 v_{z}(35,y). Heyvaerts and Priest solution overplotted. 
Fig. 12 Shaded surface of v_{y}(x,y) at t = 1000τ_{A}. The temporal evolution is shown in the movie available in the online edition. 
5.2. Evolution of v_{y}: slow magnetoacoustic waves
Figure 12 shows a shadedsurface 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 onedimensional plasma with β = 0.1, η = ν = 0.01 at t = 1000τ_{A}, i.e. similar to Fig. 3a. Thus, as in our onedimensional 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 ydirection. As expected, all amplitudes are of the order A^{2}.
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 v_{x}(x,y) at t = 100τ_{A}. b) Shaded surface of v_{x}(x,y) at t = 1000τ_{A}. The temporal evolution of the right panel is shown in the movie available in the online 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énspeed 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 shadedsurface 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 setup 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^{2} 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 shadedsurfaces 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 onedimensional 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 bulkflow 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_{0} 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 viscoresistive dissipation associated with phase mixing, producing a distinctive teardrop shape. Independently, the heating in the bottom leftcorner is the viscoresistive heating associated Alfvénwave damping along ρ_{0} = 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_{0} (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_{0} 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_{0} 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_{0} changes during our simulation, and that this drifting is in the direction of decreasing Alfvénspeed 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_{0} sometimes occurs due to the bulkflow 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_{0} 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 − T_{0} 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 − T_{0} 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, onedimensional 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: boundarydriven 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 viscoresistive damping (through the dispersion relation) and thus to viscoresistive heating, which in turn leads to the introduction of a new phenomenon: a bulk flow in the positive ydirection. Physically, this bulk flow is a direct response to the viscoresistive 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 ydirection because the boundary conditions are held fixed (v_{y} = 0). As a result of this outflow, a significant reduction in density occurs since the boundary conditions do not allow for the plasma to be replaced by an inflow. 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 viscoresistive heating is strongly dependent upon k_{I} and thus on our choice of ρ_{0}; the greater the value of ρ_{0}, the more viscoresistive heating that occurs. More viscoresistive 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 bulkflow phenomenon (which is itself directly due to the viscoresistive heating). For a ρ_{0} = 1 plasma, we find a decrease in density of about 4.6% at y = 0 but for a ρ_{0} = 5 plasma we find a decrease of about 46.4%. We also find that the rateofchange 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 viscoresistive damping, which slows the rateofchange 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 bulkflow 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 twodimensional 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 (ρ_{0} = 5) and x = 100 (ρ_{0} = 1) were identical to those of the corresponding viscoresistive damped Alfvén waves seen in our onedimensional analysis.
Our onedimensional analysis also explained the longitudinal motions (v_{y}) in our system and we again observed the boundarydriven acoustic waves (interpreted as slow magnetoacoustic waves in 2D, with speed c_{s}), our ponderomotive waves (propagating at v_{A}) and the bulkflow 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énwave 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énwave 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^{2}. 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 bulkflow phenomenon is naturally present in our system, with the density changes as predicted by our onedimensional 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^{2}) 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
Access here
Movie 2
Access here
Movie 3
Access hereAppendix A: Alfvén wave propagation in a onedimensional, ideal, β = 0 plasma
Fig. A.1 a) Transverse perturbations v_{z} for ideal (η = ν = 0) β = 0 plasma. b) Longitudinal perturbations v_{y} for ideal β = 0 plasma. In both subfigures, t = 1000τ_{A} and dotted line represents y = v_{A}t. 
Consider a cold (β = 0), ideal (η = ν = 0) plasma in our onedimensional 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 yaxis using the relation t = y / v_{A}. The rampup 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 secondorder 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 wavepressure gradients). At early times, these longitudinal motions are governed by a simplified version of Eq. (29): (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: (A.2)which is valid for 0 < y < v_{A}t. Note that these perturbations have an amplitude of and are always positive for an Alfvén wave propagating in the positive ydirection (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 rampup 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: (B.1)
Let us assume that our linear Alfvén wave can be represented as:
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_{3} = 0) we have:
The righthandside of Eq. (B.1) has two contributions, which can be calculated using our expression for b_{z}. The first term has the form:
and the second has the form:
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.:
where as before k = k_{R} + ik_{I}. Thus (from Eq. (7)) b_{z} has the form: (B.5)
Now we can explicitly determine the two contributions:
and the second has the form:
Equations (B.5) − (B.7) give equivalent solutions to Eqs. (B.2) − (B.4), and thus also lead to Eqs. (38) − (41).
Note that in both derivations, we have assumed an infinite harmonic solution for v_{z}, whereas our numerical simulations are for driven harmonic wavetrains. Thus, our solutions are only valid for c_{s}t < y < v_{A}t.
Acknowledgments
J.A.M. and I.D.M. acknowledge financial assistance from the Leverhulme Trust and Royal Society, respectively. J.A.M. acknowledges IDL support provided by STFC. J.A.M. also wishes to thank Valery Nakariakov,
Tony Arber and Hugh Hudson for helpful and insightful discussions and suggestions. The computational work for this paper was carried out on the joint STFC and SFC (SRIF) funded cluster at the University of St Andrews (Scotland, UK).
References
 Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, J. Comp. Phys., 171, 151 [Google Scholar]
 Botha, G. J. J., Arber, T. D., Nakariakov, V. M., & Keenan, F. P. 2000, A&A, 363, 1186 [NASA ADS] [Google Scholar]
 Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205 [NASA ADS] [Google Scholar]
 Browning, P. K. 1991, Plasma Phys. Cont. Nucl. Fusion, 33, 539 [Google Scholar]
 De Moortel, I., Hood, A. W., Ireland, J., & Arber, T. D. 1999, A&A, 346, 641 [NASA ADS] [Google Scholar]
 De Moortel, I., Hood, A. W., & Arber, T. D. 2000, A&A, 354, 334 [NASA ADS] [Google Scholar]
 Heyvaerts, J., & Priest, E. R. 1983, A&A, 117, 220 [NASA ADS] [Google Scholar]
 Hollweg, J. V. 1971, J. Geophys. Res., 76, 5155 [NASA ADS] [CrossRef] [Google Scholar]
 Hollweg, J. V., & Yang, G. 1988, J. Geophys. Res., 93, 5423 [NASA ADS] [CrossRef] [Google Scholar]
 Hood, A. W., Brooks, S. J., & Wright, A. N. 2002, Proc. Roy. Soc. A, 458, 2307 [Google Scholar]
 Hood, A. W., Ireland, J., & Priest, E. R. 1997a, A&A, 318, 957 [NASA ADS] [Google Scholar]
 Hood, A. W., GonzalésDelgado, D., & Ireland, J. 1997b, A&A, 324, 11 [NASA ADS] [Google Scholar]
 Ionson, J. A. 1978, ApJ, 226, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Ireland, J. 1996, Ann. Geophys., 14, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, M. A., & Roberts, B. 1986, ApJ, 301, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Malara, F., Primavera, L., & Veltri, P. 1996, ApJ, 459, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Mocanu, G., Marcu, A., Ballai, I., & Orza, B. 2008, Astron. Nachr., 329, 780 [Google Scholar]
 Narain, U., & Ulmschneider, P. 1990, Space Sci. Rev., 54, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Narain, U., & Ulmschneider, P. 1996, Space Sci. Rev., 75, 453 [NASA ADS] [CrossRef] [Google Scholar]
 Nakariakov, V. M., & Oraevsky V. N. 1995, Sol. Phys., 160, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Nakariakov, V. M., Roberts B., & Murawski, K. 1997, Sol. Phys., 175, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Nakariakov, V. M., Roberts B., & Murawski, K. 1998, A&A, 332, 795 [NASA ADS] [Google Scholar]
 Ofman, L., & Davila, J. M. 1995, J. Geophys. Res., 100, 23413 [NASA ADS] [CrossRef] [Google Scholar]
 Ofman, L., & Davila, J. M. 1997, ApJ, 476, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Ofman, L., Klimchuk, J. A., & Davila, J. M. 1998, ApJ, 493, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1991, ApJ, 376, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Poedts, S., & Boynton, G. C. 1996, A&A, 306, 610 [NASA ADS] [Google Scholar]
 Priest, E. R. 1982, Solar magnetohydrodynamics (D. Reidel Publishing. Company) [Google Scholar]
 Ruderman, M. S., & Roberts B. 2002, ApJ, 577, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S., Nakariakov, V. M., & Roberts B. 1998, A&A, 338, 1118 [NASA ADS] [Google Scholar]
 Smith, P. D., Tsiklauri, D., & Ruderman, M. S. 2007, A&A, 475, 1111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Threlfall, J., McClements, K. G., & De Moortel, I. 2010, A&A, 525, A155 [Google Scholar]
 Tsiklauri, D., & Nakariakov, V. M. 2002, A&A, 393, 321 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tsiklauri, D., Arber, T. D., & Nakariakov, V. M. 2001, A&A, 379, 1098 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tsiklauri, D., Nakariakov, V. M., & Arber, T. D. 2002, A&A, 395, 285 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tsiklauri, D., Nakariakov, V. M., & Rowlands, G. 2003, A&A, 400, 1051 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ulmschneider, P., Priest, E. R., & Rosner, R. 1991, Mechanisms of Chromospheric and Coronal Heating (Berlin: Springer) [Google Scholar]
 Verwichte, E. 1999, PhD Thesis: Aspects of Nonlinearity and Dissipation in Magnetohydrodynamics, The Open University [Google Scholar]
 Verwichte, E., Nakariakov, V. M., & Longbottom, A. W. 1999, J. Plasma Phys., 62, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Wilkins, M. L. 1980, J. Comp. Phys., 36, 281 [Google Scholar]
All Figures
Fig. 1 a) Longitudinal perturbations v_{y} for β = 0.1 plasma with η = ν = 0. Dotted and dashed lines represent y = v_{A}t and y = c_{s}t respectively. The blue line denotes the maximum amplitude: . b) Blowup of region 100 ≤ y ≤ 200 for v_{y}. In both subfigures, t = 1000τ_{A}. 

In the text 
Fig. 2 a) Transverse perturbations v_{z} for η = 0.01, β = 0.1 plasma (ν = 0), where dashed lines represent the envelope v_{z} = ± Aexp(k_{I}y). b) Longitudinal perturbations v_{y} for same plasma, where dashed line represents y = c_{s}t and blue horizontal line denotes v_{y} = 0. In both subfigures, t = 1000τ_{A} and dotted line represents y = v_{A}t. 

In the text 
Fig. 3 a) Longitudinal perturbations v_{y} 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) Blowup of v_{y} over region 0 ≤ y ≤ 100 for ν = 0.01, η = 0 plasma only, where crosses indicate grid point resolution. Black line indicates v_{y} = 0. In both subfigures, t = 1000τ_{A}. 

In the text 
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 = c_{s}t (y = v_{A}t) and the blue line denotes ρ = ρ_{0}. 

In the text 
Fig. 5 Longterm evolution of a) density, b) temperature, c) pressure gradient (∂p / ∂y), and d) k_{I}, over 0 ≤ t ≤ 100 000τ_{A}, where the different colours represent different times. The colour bar denotes the different times. 

In the text 
Fig. 6 a) Transverse perturbations v_{z} for η = ν = 0.01, β = 0.1 plasma with ρ_{0} = 5, where red dashed lines represent the envelope v_{z} = ± Aexp(k_{I}y). b) Longitudinal perturbations v_{y} for same plasma, where black dashed line represents y = c_{s}t and blue horizontal line denotes v_{y} = 0. In both subfigures, t = 1000τ_{A} and the vertical black dotted line represents y = v_{A}t. 

In the text 
Fig. 7 Longterm 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. 

In the text 
Fig. 8 Evolution of ∂ρ / ∂t at y = 1 for a) ρ_{0} = 5 and b) ρ_{0} = 1 plasmas. 

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

In the text 
Fig. 10 Contours of v_{z} 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 online edition. 

In the text 
Fig. 11 Plot of v_{z}(35,y). Heyvaerts and Priest solution overplotted. 

In the text 
Fig. 12 Shaded surface of v_{y}(x,y) at t = 1000τ_{A}. The temporal evolution is shown in the movie available in the online edition. 

In the text 
Fig. 13 a) Contour of v_{x}(x,y) at t = 100τ_{A}. b) Shaded surface of v_{x}(x,y) at t = 1000τ_{A}. The temporal evolution of the right panel is shown in the movie available in the online edition. 

In the text 
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}. 

In the text 
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}. 

In the text 
Fig. 16 Contour of T − T_{0} 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 − T_{0} 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. 

In the text 
Fig. A.1 a) Transverse perturbations v_{z} for ideal (η = ν = 0) β = 0 plasma. b) Longitudinal perturbations v_{y} for ideal β = 0 plasma. In both subfigures, t = 1000τ_{A} and dotted line represents y = v_{A}t. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.