IMPLICATIONS OF A NON-MODAL LINEAR THEORY FOR THE MARGINAL STABILITY STATE AND THE DISSIPATION OF FLUCTUATIONS IN THE SOLAR WIND

, , and

Published 2010 April 27 © 2010. The American Astronomical Society. All rights reserved.
, , Citation Enrico Camporeale et al 2010 ApJ 715 260 DOI 10.1088/0004-637X/715/1/260

0004-637X/715/1/260

ABSTRACT

A magnetized plasma with anisotropic particle distributions may be unstable to a number of different kinetic instabilities. The solar wind is often observed in a state which is close to that implying instability, i.e., in a marginal stability state. Normal-mode linear theory predicts that fluctuations in a stable plasma damp exponentially. The non-modal approach for a linearized system differs from a normal-mode analysis by following the temporal evolution of some perturbed equilibria, and therefore includes transient effects. We employ a non-modal approach for studying the stability of a bi-Maxwellian magnetized plasma using the Landau fluid model, which we briefly describe. We show that bi-Maxwellian stable equilibria can support transient growth of some physical quantities, and we study how these transients behave when an equilibrium approaches its marginally stable condition. The results obtained with a non-modal approach are relevant to a re-examination of the concept of linear marginal stability. Moreover, we highlight some aspects of the dissipation of turbulent fluctuations, which suggest that the non-modal approach should be included in future studies.

Export citation and abstract BibTeX RIS

1. INTRODUCTION

Linear theory represents a powerful tool for the interpretation and understanding of many space plasma properties observed by in situ spacecraft. For instance, the temperature anisotropy of the solar wind is thought to be bounded by values that are consistent with the stability thresholds derived within the linear theory of the Vlasov–Maxwell set of equations. The broadly accepted view is that a macroscopic property of the plasma (such as temperature, density, or mean velocity) can be constrained by the nonlinear feedback associated with a linear instability. This is because the primary consequence of any instability is to reduce the amount of free energy that drives the instability, relaxing the plasma toward a marginally stable condition. This is equivalent to saying that the plasma is unlikely to be found in an unstable state because it tends to change its macroscopic properties in a way that would lead to a linearly stable condition.

In the solar wind, it is argued that the expansion of the plasma from the Sun in a radially decreasing magnetic field should produce much higher values of temperature anisotropy (with respect to the background magnetic field) than those observed. Also, it has been shown that the highest values of anisotropy observed in the solar wind are consistent with the thresholds of linear kinetic instabilities driven by temperature anisotropies. This is true both for protons (Kasper et al. 2002; Hellinger et al. 2006; Marsch et al. 2006; Matteini et al. 2007; Bale et al. 2009) and electrons (Gary et al. 2005; Stverak et al. 2008), for a large range of plasma beta, and both for T/T>1 (whistler and mirror instabilities) and for T/T < 1 (firehose instability).

The physical behavior of an unstable anisotropic collisionless plasma subject to the electron firehose instability, that relaxes toward marginal stability, has been elucidated in Camporeale & Burgess (2008), with fully nonlinear particle-in-cell simulations. What emerges is that the plasma state is likely to be found bouncing around the marginal stability threshold, due to the competition of two mechanisms: the reduction of the anisotropy and plasma beta caused by the development of the firehose instability (above the threshold), and the increase of these quantities due to the damping of magnetic fluctuations (below the threshold) that result in the energization of the particles, predominantly in one direction.

Despite the success of the linear theory to delineate the macroscopic properties in which the solar wind plasma is more likely to be found (i.e., in stable conditions) and the good agreement between linear theory predictions and solar wind data, there are at least two issues which are problematic from the point of view of linear theory.

First, the greater part of the protons and electrons in the solar wind, at any distance from the Sun, is observed to be isotropic or only slightly anisotropic (Hellinger et al. 2006; Stverak et al. 2008), i.e., in a condition where no anisotropy instability can be excited and therefore the aforementioned argument associated with the linear constraining mechanism cannot be invoked. In other words, the fact that linear instabilities do not allow the plasma to develop anisotropies higher than certain values does not explain why most of the plasma is found to be very far from those values. Since the effect of the expansion should result in a continuous increase of the anisotropy, one would expect to find the plasma most of the time close to marginal stability, unless an isotropizing mechanism is in operation even before it reaches the unstable condition. Observational results have suggested, as a possible explanation, that the collisional age is related to the isotropization of thermal electrons (Salem et al. 2003; Stverak et al. 2008), but the relative importance of collisions and instabilities is still unknown.

Second, a relatively high occurrence of short wavelength magnetic fluctuations with small amplitude is persistently found in the plasma, even in stable conditions (Bale et al. 2009). These fluctuations could result from a balance between a nonlinear transfer and dissipation but they could not be accounted for on the basis of the modal linear theory alone, as they would be expected to damp exponentially in time.

The aim of this paper is to present a non-modal approach to the linear stability problem for a collisionless plasma in an uniform magnetic field that, by offering a new framework for the understanding of linear marginal stability, will reconcile those two apparent observational contradictions with a consistent physical interpretation.

The non-modal approach will also have implications for the physical mechanism for kinetic damping of fluctuations in the turbulent dissipation range. This is an area of increasingly active research, and it has been argued that the linear approximation (i.e., the assumption that perturbed quantities have a much smaller amplitude than the equilibrium ones) might be used in this scenario (Gary & Borovsky 2004). There is a general agreement on the fact that the nonlinear cascade of turbulent fluctuations should take into account the onset of kinetic effects at a certain spatial scale, around the ion Larmor radius (Matthaeus et al. 2008). It is indeed observed that the power density of magnetic fluctuations undergoes an abrupt steepening for wavenumbers above a certain value (Leamon et al. 1998; Hamilton et al. 2008; Sahraoui et al. 2009), and it is thought that kinetic effects (or even simply dispersion) might be responsible for the steepening. However, a complete understanding on exactly why this happens and what are the important parameters that determine at which scale one should expect the kinetic effects to produce a change of the properties of the power spectrum is still missing. Also, it is a matter of controversy whether the power spectrum of magnetic fluctuations at high frequencies follows a power law in wavenumbers (Sahraoui et al. 2009), or is better fitted by an exponential decay (Alexandrova et al. 2009). The first result is usually interpreted as the signature that a nonlinear cascade is able to proceed below electron scales, while the latter observations would suggest a damping apparently in agreement with the linear theory. However, there is no general agreement on which normal modes would be able to remain undamped at high frequencies. The current debate is focused mainly on the oblique-propagating kinetic Alfven waves (KAWs; Howes et al. 2008a) and the parallel whistler (Shaikh 2009; Gary & Smith 2009).

Although the study of kinetic instabilities, and their application to the solar wind plasma, has been traditionally kept separate from the studies addressing the turbulent dissipation range, we regard those two aspects of plasma physics as interlocked. In fact, they describe the same problem from two distinct viewpoint, because the understanding of turbulent dissipation at small scales can be properly addressed only by including in the energy balance the injection of turbulent fluctuations due to the development of kinetic instabilities.

Different approaches to this problem include the use of gyrokinetic linear theory (Howes et al. 2006, 2008a) and simulations (Howes et al. 2008b), Hall-MHD (Krishan & Mahajan 2004; Galtier & Buchlin 2007; Shaikh and Shukla 2009) and particle-in-cell simulations (Gary et al. 2008). A simple diffusive model for the nonlinear cascade in the inertial range (i.e., for spatial scales that can be precisely studied with a MHD description) was suggested by Zhou & Matthaeus (1990). This model has been successively applied by Li et al. (2001) for addressing the transition between inertial and dissipation ranges. Their conclusion was that any kinetic linear damping mechanism would not be able to reproduce the observed power spectra of magnetic fluctuations, because it would produce a steep cutoff instead of a power law, for high wavenumbers.

We will show that, by approaching the linear theory with a non-modal formalism, the possibility that a completely linear damping mechanism is responsible for the persistence of a power-law spectrum at high wavenumbers, might still be valid. Moreover, we argue that the results based on a non-modal approach clarify the relationship between different eigenmodes, and clearly establish that it is very unlikely that the problem of turbulent dissipation will be understood by studying the damping properties of one single normal mode (or alternatively of a spectrum of a single mode), as previous works attempted to do. We will argue that, at kinetic scales, the evolution of a small perturbation depends very much on the linear interactions between many different modes, including heavily damped ones. This linear interaction is only found when the non-modal properties of the system are included. In a sense, this argument makes the competition between whistler and KAWs an ill-posed problem. Also, from a purely computational point of view, we will show that our results question the possibility of exciting one single eigenmode.

1.1. Modal and Non-modal Approaches

Plasma stability theory is traditionally treated as a normal-mode ("modal") analysis. This means that the focus is on the eigenmodes of a slightly perturbed equilibrium, which are assumed to grow or damp exponentially in time. The most rapidly growing (or the least damped) eigensolutions are the object of greatest interest in the normal-mode analysis, which therefore investigates the stability problem purely in its time-asymptotic solutions, regardless of the details of the perturbation imposed on the initial equilibrium.

The non-modal approach differs from the normal-mode analysis in treating the linear stability as an initial value problem. This allows the study of transient phenomena and to follow the evolution of a perturbed system in time. This evolution depends, of course, on how the initial equilibrium is perturbed.

The non-modal approach to the plasma stability problem can in some cases be crucial to the understanding of the evolution of a plasma in the linear regime. Especially when studying a system in stable conditions, the normal-mode analysis misses the phenomenon of transient growth in time of some physical quantities. This is a well-known effect which has been long studied in hydrodynamic flows (Schmid 2007). It is related to the spectral properties of the linear operator that describes the evolution of the system in time. In particular, when an operator is non normal, i.e., it does not commute with its adjoint, the norm of an initial perturbation applied to the equilibrium can grow in time by large factors, before decaying, even if all the eigenmodes of the operator are damping (i.e., the equilibrium is stable).

Despite the fact that transient growth has been known of for a long time (Trefethen et al. 1993), it has not been emphasized in the theory of plasma stability. We have conducted the first study of transient growth for a stable kinetic plasma in a previous paper (Camporeale et al. 2009), showing that this phenomenon appears to be related to the kinetic regime of a plasma; it is indeed more accentuated for higher plasma beta and shorter spatial scales. In order to include kinetic effects, we used a Landau fluid (LF) model, which incorporates linear Landau damping, and finite Larmor radius (FLR) corrections.

1.2. Aims of the Paper

The purpose of this paper is threefold. First, we extend the work of Camporeale et al. (2009) to an anisotropic plasma, for conditions typical of the solar wind. By doing so, we will highlight the inadequacy of the normal-mode analysis for a stable kinetic plasma.

Second, we will clarify the relation between transient growth and marginal stability, and we will propose a physical picture where the two observational contradictions mentioned above will be reconciled with the linear theory.

Third, we will discuss the importance and appropriateness of a non-modal linear theory for the understanding of the dissipation of turbulent fluctuations. It should be the goal of a future complete theory of solar wind turbulence to elaborate a unified theoretical framework that links the physics of kinetic instabilities, and the damping of turbulent fluctuations in a consistent manner, and we believe that such a theory will benefit from embracing a non-modal approach to the stability problem, such as the one presented in this paper.

The paper is organized as follows. In Section 2, we briefly describe the LF model that has been used to obtain the results presented in the paper. The model has been described and commented in length elsewhere (Passot & Sulem 2006, 2007; Sulem & Passot 2008), hence we will just present the main features of the model, and we report the set of equations in the Appendix. In Section 3, we will introduce the methodology of the non-modal approach, and we will define some important quantities for our analysis. Section 4 describes the results of our study, with emphasis on the relationship between transient growth and marginal stability and on the non-modal approach for the study of turbulent dissipative fluctuations. A final discussion, with suggestions for future work, is reported in Section 5.

2. THE FLR-LANDAU FLUID MODEL

The idea of incorporating Landau damping in a set of fluid equations was introduced by Hammett & Perkins (1990) and later developed in Snyder et al. (1997), where the fluid hierarchy obtained from the drift kinetic equation is closed at the level of the third or fourth order moment. These models are limited to scales large compared with the ion gyroradius. Other models, called gyrofluid models, consider the fluid hierarchy obtained from the gyrokinetic equation, providing a set of equations for fluid moments suitable for the description of sub-Larmor radius scales (Brizard 1992; Dorland and Hammett 1993; Beer & Hammet 1996). The equations of gyrofluid models, however, are not written in the physical coordinates but in the gyrocenter variables, making their interpretation more difficult. A simpler formulation retaining hydrodynamic nonlinearities together with a linear approximation of FLR contributions was recently developed by deriving equations for the hydrodynamic moments directly from the Vlasov–Maxwell system (Goswami et al. 2005; Sulem & Passot 2008; Passot & Sulem 2007). This is the model used in this paper. The hierarchy of fluid equations is closed at the level of the fourth order moment. In its linearized version, the plasma dispersion relation is approximated by a suitable Padé approximant, and this makes it possible to cast the linear set of equations as a standard eigenvalue problem. In addition to the hierarchy closure, this approach involves the modeling of FLR effects by expressing the non-gyrotropic part of tensors such as pressures, heat fluxes, or fourth order moments in terms of lower-rank moments, in a way consistent with the linear kinetic theory in the low-frequency limit epsilon ∼ ω/Ωi ≪ 1, for both quasi-transverse fluctuations (k||/kepsilon) with no condition on krL (as in gyrokinetic and gyrofluid approaches), but also for hydrodynamic scales with k||k ≪ 1/rL. Here, Ωi denotes the ion cyclotron frequency and rL the ion Larmor radius. At large scales, the model, which then reduces to usual anisotropic MHD, also captures the fast waves, in contrast with gyrofluids. The frequency and damping rates of low-frequency waves are accurately described in a range of scales that extends to small (sub-Larmor radius) scales when the propagation direction is almost perpendicular to the ambient magnetic field (according to the gyrokinetic scaling). The complete model is quite involved and is thoroughly described in Passot & Sulem (2007). The full set of equations is reported in the Appendix, for completeness. The total number of variables is 16, hence the linear problem reduces to the formulation of a 16 × 16 complex matrix.

3. NON-MODAL APPROACH

In this section, we introduce the mathematical tools (mostly taken from linear algebra) and the methodology that will be used in the rest of the paper. For a more complete description of the non-modal stability theory in the context of hydrodynamics, we refer to the review by Schmid (2007); several issues related to the non-normality of a linear operator (that will be a central point of our discussion) have been described in great depth in the monograph by Trefethen & Embree (2005).

The LF model can be linearized as usual, by writing each physical quantity as a sum of an equilibrium and a perturbed contribution (subscript 0 and 1, respectively): ϕ = ϕ0 + εϕ1, assuming that ε ≪ 1, and neglecting terms of order higher than one in ε. Once the first-order variables are Fourier-decomposed (dropping the subscript) $\phi (\mathbf {r},t)=\tilde{\phi }\exp [i(\mathbf {k\cdot r})]$, the linear LF model can be cast as a set of ordinary differential equations for the complex amplitudes $\tilde{\phi }$:

Equation (1)

A is a time-independent operator, that takes the form of a 16 × 16 sparse complex matrix, whose entries depend on the properties of the plasma (protons and electrons temperature anisotropy and plasma beta) and on the magnitude and angle of propagation of the wavevector k.

The fact that the LF model describes the plasma through a set of only 16 variables makes the analysis of the matrix A computationally affordable without any particular method used for large matrix manipulations. Accordingly, all the results presented in this paper have been produced using in-built routines of MATLAB.

The solution of Equation (1) is given by

Equation (2)

where $\tilde{\phi }(0)$ is the state vector of the initial perturbation and the exponential of the matrix (which is defined as $e^{\mathbf {A}t}=\mathbf {I}+\mathbf {A}t+\frac{1}{2}(\mathbf {A}t)^2+\cdots$) completely determines the evolution in time of the initial state.

If A commutes with its adjoint AA* = A*A, then A is said to be normal and its eigenvectors form a complete orthogonal set. The definition of normality depends on which norm one refers to. In fact, we recall that the norm of a vector is defined as $\Vert u\Vert =\sqrt{\langle u,u\rangle }$, where the inner product 〈 · 〉 on $\mathbb {C}^n$ is 〈x, y〉 = y*Mx and M is a positive-definite matrix. It is known (see, e.g., the Appendix of Farrell & Ioannou 1993 for a proof) that for any operator A there is at least one matrix M that defines an inner product that makes the eigenvectors of A orthogonal. In other words, any non-normal operator can be made normal by choosing a different definition of inner product, from which the definitions of norm and orthogonality follow. We use in this paper the two-norm $\Vert u\Vert _2=\sqrt{\sum _{i=1}^{16}|u_i|^2}$, which is the case for M equal to the identity matrix. A different choice of M would be equivalent to a transformation of variables $v={\bf M}^\frac{1}{2}u$, so that the new norm $\Vert u\Vert _M=\sqrt{u^*{\bf M}u}$ is equal to ||v||2. Note that the new set of variables vi, which are a linear combination of ui, do not have, in general, any physical meaning. This explains the apparent contradiction that one operator can support a transient growth of the state vector, when this is measured with one norm, but such transient growth will not appear if the definition of norm is changed.

In summary, although, by choice of a norm, it is possible to recast an operator as normal, the corresponding variables in this case are not the basic variables of the underlying physical system. Of course, the mathematical definition of normality does not affect the transient growth effect in a physical sense. A simple way to check that a transient growth is not a mathematical artifact due to a specific choice of a norm is to measure the growth of single physical quantities, such as the magnetic field, instead of the entire state vector. This will be the approach mostly used in this paper.

Although the two-norm ||u||2 (the subscript 2 will be dropped) does not provide an information on single variables, it indicates whether the perturbation complies or not with the requirement of being much smaller than the equilibrium quantities. What is important is that a large two-norm of a perturbation implies a deviation from the assumption of linearity, and therefore might result in the triggering of nonlinear effects.

3.1. Spectral Abscissa and Numerical Abscissa

In order to study how the plasma responds to a small perturbation we introduce the quantity G(t), which measures the amplification (or reduction) in time of an initial perturbation:

Equation (3)

This quantity clearly depends on the details of the initial perturbation $\tilde{\phi }(0)$, but it is always bounded from above by the quantity ||eAt||, which defines the supremum of G(t) for any possible $\tilde{\phi }(0)$. The behavior of ||eAt|| depends on the spectral properties of A. If one indicates with σ(A) the set of eigenvalues of A, then the quantity

Equation (4)

is referred to as the spectral abscissa. To determine this quantity is the main objective of the normal-mode analysis, because it gives the growth rate of the most unstable mode (for an instability) or the damping rate of the least damped mode (for a stable plasma). Hence, the time-asymptotic evolution of any initial perturbation is given by eα(A)t. Moreover, for a normal operator, the spectral abscissa α bounds the quantity G(t), for any time t ⩾ 0:

Equation (5)

In general, however, for a non-normal operator, Equation (5) holds only in the limit t, and there is no exact formula for the quantity ||eAt||, that can only be approximately estimated for t>0 (different approximations can be found in Trefethen & Embree 2005).

But a rigorous formula exists for the growth of ||eAt|| at t = 0. This is referred to as the numerical abscissa:

Equation (6)

Note that since (A + A*) is Hermitian, η(A) is real. The numerical abscissa provides the information about the highest possible growth rate of any initial perturbation, at time t = 0. For a normal operator, it follows from Equation (6) that η(A) = α(A).

What is surprising is that, for a non-normal operator, even if all the normal modes are damping (i.e., α(A) < 0) the numerical abscissa can be positive, hence allowing the quantity G(t) to grow for some particular initial perturbations. It has been shown in Camporeale et al. (2009) that G(t) can indeed reach values of about 103–104 over short periods of time, for a Maxwellian plasma described by the LF model.

The fact that a perturbation could grow while the eigenmodes of the operator decay in time may appear counter-intuitive. It is purely due to the superposition of non-orthogonal eigenvectors, and we refer to Figure 2 in Schmid (2007) for a graphical paradigmatic explanation of this effect.

3.2. Pseudospectra

A key property of systems with a non-normal operator is that the eigenvalues may be highly sensitive to small perturbations. In order to quantify the effect of small perturbations on a linear operator, the concept of pseudospectra has been introduced. Although pseudospectra will not be used for deriving the results presented in the next section, we introduce this concept here because we think that it is the easiest way to understand why normal and non-normal operators behave differently.

There are at least four definitions of pseudospectra that have been shown to be mathematically equivalent (Trefethen & Embree 2005). We report here the two that have a most immediate link with a physical interpretation. If one employs the convention that the spectrum σ(A) (the set of eigenvalues) is formed by complex values z for which ||(zIA)−1|| = , where I is the identity operator, the ε-pseudospectrum σε(A) is defined as the set of $z\in \mathbb {C}$ such that

for any ε>0.

This is equivalent to saying that z is an eigenvalue of the perturbed operator (A + E) for some operator E with ||E|| < ε, i.e.,

Let us now comment on these definitions and their meaning in a more physical sense. Our linear operator studies the evolution of a small perturbation applied to an equilibrium, namely a bi-Maxwellian plasma in an uniform magnetic field. However, one could argue that small discontinuities or inhomogeneities of the equilibrium quantities might be modeled as small disturbances to the linear operator. The question that pseudospectra quantitatively answer is: how does the set of eigenvalues change under the effect of small perturbations of the operator? From definitions 1 and 2, one see that the ε-pseudospectrum σε tends precisely to the standard spectrum σ, when ε → 0. In the complex plane, the ε-pseudospectra are the open subset that contains all the eigenvalues of the perturbed operator (A + E), for any possible perturbation E, such that ||E|| < ε. Therefore, they give a measure of the distortion of the spectrum due to the perturbation applied to the operator. It is straightforward to prove that the ε-pseudospectra are a nested set. That is, $\sigma _{\varepsilon _1}(\mathbf {A})\subseteq \sigma _{\varepsilon _2}(\mathbf {A})$ for ε1 ⩽ ε2. An example of pseudospectra of our LF operator is given in Figure 1, for an isotropic plasma with β = 10 and k = 0.5 (quasi-perpendicular propagation). Each contour corresponds to values of ε = 10−6.8, 10−6.6, ..., and10−5.6. Only a part of the complex plane is shown, with 9 of the 16 eigenvalues visible. The important point here is to understand that perturbations of a certain size make the ε-pseudospectrum so large that it would contain many eigenvalues of the original unperturbed operator. In a sense, this means that such perturbations distort the problem in such a way that the information about the modes of the unperturbed operator becomes irrelevant. In this respect, the different behavior between normal and non-normal operator is most evident. The ε-pseudospectrum of a normal operator is defined as the union of open balls about the points of the spectrum:

where "dist" indicates the distance of a point to a set in the complex plane. For a normal operator, the curves of the pseudospectra can be computed straightforwardly once the eigenvalues are known. We show in Figure 1, with dotted lines, how the contours of σε would be for the same LF operator, if it were normal, for ε = 0.00075. This contour is qualitatively similar to the (true) contour for ε = 10−5.8, embracing all the 9 eigenvalues. This is the crucial point about the different behavior between normal and non-normal operators. A perturbation of the order 10−5.8 is sufficient to achieve a distortion of the spectrum that, if the operator would have been normal, would have required a perturbation of the order 7.5 × 10−4, i.e., about 470 times higher.

Figure 1.

Figure 1. Contours of the ε-pseudospectrum for a plasma with β = 10, k = 0.5, and T/T = 1. Contours are plotted for log10ε = −6.8, − 6.6, ..., and − 5.6. The dotted line is how the ε-contour would appear if the operator were normal, for ε = 0.00075.

Standard image High-resolution image

The displacement of the spectrum of an operator subject to small perturbations implies two facts. On one hand some, if not all, of the modes become somehow coupled, i.e., their relative distance in the complex plane can be modified, and their properties (like phase speed and damping rate) changed. On the other hand, from a computational point of view, the non-normality of the operator almost completely rules out the possibility of exactly exciting one single eigenvector. Small errors, due for instance to digits truncation or approximation, can result in the excitation of a "pseudomode," that could lie in the complex space, very far from the mode that was intended to be excited, and that is in reality a superposition of different non-orthogonal modes. This might lead to an initial transient behavior, which is not unusual in numerical simulations, even though it is seldom commented.

4. RESULTS

In this section, we apply the linear LF model to a stable proton–electron (p, e) bi-Maxwellian plasma in an uniform magnetic field. We are interested in how the distance from marginal stability affects the evolution of linear perturbations in the plasma and wish to address transient behavior that is missed by the standard modal analysis. For all the results shown in this section, the protons are considered isotropic with Tp = Tp = Te. All quantities are therefore referred to electrons, and the subscript e is dropped.

We focus on quasi-perpendicular waves, since the LF model has a domain of validity that extends to small scales only for oblique wavevectors. The angle of propagation for the wavevector k with respect to the background magnetic field is θ = tan−1(1000).

It is known that the marginal stability thresholds (i.e., the curves for which the growth rate is exactly null) obey a law of the form: $\frac{T_\perp }{T_\parallel }=1+\frac{S}{\beta ^\alpha }$, where S and α are constants. We have derived the stability thresholds for the LF model, and we have computed the fitting parameters using a Levenberg–Marquardat method (Press et al. 1986). The result is

Equation (7)

and

Equation (8)

The values for T>T are not dissimilar from those obtained from the Vlasov–Maxwell equations for the electron firehose instability (S = −1.29, α = 0.98), which is known to yield highest growth rate for quasi-perpendicular propagation (Camporeale & Burgess 2008). The parameters for T < T are instead not comparable with those obtained for the mirror and whistler instabilities (Gary & Wang 1996), because in this case the whistler instability, which is parallel propagating, sets the instability threshold (Gary & Karimabadi 2006), and the LF model becomes invalid for strictly parallel wavevectors. We therefore will focus mainly on the case T>T.

4.1. Transient Growth and Marginal Stability

It has been shown in Camporeale et al. (2009) that an isotropic Maxwellian plasma is able to sustain large transient growth of an initial perturbation that will eventually decay in later times.

As we have seen in Section 3, the spectral abscissa α (which is negative, for a stable plasma) provides the information about the late-time damping rate of the perturbation. That is, any initial fluctuation will decay as eαt, for large times. An interesting point is to actually analyze the time at which the system starts to behave as predicted by the normal-mode analysis.

We show in Figure 2 the amplitude of the y-component of the magnetic field, normalized to its initial value, for two particular initial conditions. The wavevector k is chosen equal to 1 and T/T = 1. The parameter β is equal to 1 in the top panel and β = 5 in the bottom panel. What emerges is that, for both cases, the behavior is highly oscillatory, and no damping is evident before TΩi = 105. This implies that, for the initial perturbation to decay, the plasma, once perturbed, should evolve without encountering any further perturbation for a very large time, which is quite unrealistic. The damping predicted by modal analysis (i.e., eαt) is shown in dashed line.

Figure 2.

Figure 2. Evolution in time of the absolute value of the amplitude of the y-component of the magnetic field, normalized to its initial value, for one particular choice of initial condition. The parameters used are: k = 1, T/T = 1, and β = 1 (top panel) and β = 5 (bottom panel). The curve in dashed line is the evolution eαt, predicted by modal theory.

Standard image High-resolution image

Before proceeding, let us clarify how the initial perturbations used in Figure 2 were chosen. Let us recall that the singular value decomposition (SVD) of an operator A is given by

Equation (9)

where U and V are unitary matrices, and ${\rm {\bSigma }}$ is a diagonal matrix that contains the singular values of A.

If one wants to find the initial perturbation $\tilde{\phi }(0)$, that at a specific time Θ undergoes an amplification G(Θ) equal to ||eAΘ|| (i.e., the highest possible amplification at time Θ), it is sufficient to calculate the SVD of eAΘ and to identify the column vector of V associated with the highest singular value. The initial perturbations chosen to produce Figures 2 and 3 are the ones that approximately attain the highest amplification at a certain time. One could argue that these initial conditions are very special and cannot represent the totality of all the possible initial states. This is certainly true. However, it is not feasible to apply the non-modal approach to a very large set of different conditions (that will never be large enough to represent "all" cases). Besides, the methodology that we use, i.e., to focus on only certain particular cases that represent the "worst case scenario," is common also to the modal approach, whose results are supposed to be valid for any initial conditions (albeit restricted to asymptotic times), but still represent only the scenario in which the fastest growing (or the least damped) mode is actually excited. Later, we describe an attempt to approach the problem in a more statistical manner.

Figure 3.

Figure 3. Evolution in time of the absolute value of the amplitude of the y-component of the magnetic field, normalized to its initial value, for one particular choice of initial condition. The parameters used are: T/T = 0.65, β = 5, and k = 1 (top panel), k = 5 (central panel), and k = 10 (bottom panel). The curve in dashed line is the evolution eαt, predicted by modal theory.

Standard image High-resolution image

Figure 2 shows the results of two cases for an isotropic plasma, with different β, but what is perhaps more interesting is to look at what happens when the plasma is closer to a marginal stability condition, given by Equations (7) and (8). For β = 5, Equation (7) predicts that the plasma is marginally stable when T/T ∼ 0.6. We show in Figure 3 three cases of transient growth with β = 5 and T/T = 0.65, for k = 1, 5, 10. The damping rate is, respectively, α = −1.02 × 10−5, − 7.2 × 10−7, and − 3.9 × 10−7.

Since the damping rate is so low, what is expected is that any perturbation applied to the plasma would remain unchanged in amplitude for long times. This is indeed what is observed. What is rather interesting, however, is that the initial perturbations shown in Figure 3 undergo a transient growth at early times and therefore remain amplified for long times. In other words, a transient growth effect is able to amplify a perturbation even at marginal stability condition, and the fact that the damping rate is very small allows to the amplified perturbation to survive for long times. For instance, in the top panel of Figure 3, the initial value of By gets amplified by a factor between 10 and 100 and remains so at least for 105 ion gyroperiods. We have checked the results of Figures 2 and 3 against the prediction of a nonlinear code, based on the same LF model, and they agree very well. This suggests that the transient growth is not a spurious effect. Note also that the apparent large gradients in the magnetic field in Figures 2 and 3 are just an effect of the log–log plot, and the curves appear very smooth on a linear plot.

We now try to quantify the importance of transient growth effects in a statistical sense for the following range of parameters: β ∈ [1.5, 10], T/T ∈ [0.2, 1.2] and k = 1.

We have divided the (β, T/T) space in a 100 × 100 grid. For each point in the stable region, we have generated 10,000 random initial perturbations, where each of the 16 components of the state vector is independently chosen from a uniform distribution in the interval [ − 1, 1] × [ − 1, 1]i in the complex plane. In Figure 4, we show the mean value of G(t) at times t = 1, 10, 100, and1000 in log10 scale. We recall that G(t) measures the amplifications of the two-norm on the vector, hence considering all the 16 variables equally. One can note that at time t = 1, G(t) increases for higher β, almost independently of the value of the anisotropy. Transient growth is therefore not a sporadic event, even for an isotropic plasma. For later times, G(t) has still large values, and the highest peaks tend to be localized at the edge of marginal stability, for increasing time. In Figure 5, we show the average (left panels) and the maximum (right panels) value of δB/B0 calculated across the 10,000 different initial perturbations. The highest magnetic fluctuations are evident in proximity of the T>T instability threshold for any time. Once again the transient growth of magnetic fluctuations appears not to be a sporadic and unusual event that happens only for carefully chosen initial perturbations, but a feature that is persistent in time.

Figure 4.

Figure 4. Average value of G(t) computed over 10,000 randomly generated initial perturbations, shown at four different times, in the (β, 1 − T/T) space. Logarithmic scale.

Standard image High-resolution image
Figure 5.

Figure 5. Average (left panel) and maximum (right panel) values of δB/B0 computed over 10,000 randomly generated initial perturbations, shown at four different times, in the (β, 1 − T/T) space. Logarithmic scale.

Standard image High-resolution image

In the light of these results, we suggest the following scenario that might explain both the presence of fluctuations at high wavenumbers in stable conditions and the fact that most of the solar wind is observed close to temperature isotropy, despite the expansion.

The increase of anisotropy is expected, in fluid theories, due to the conservation of one or more adiabatic invariants. For instance, the CGL approximation (Chew et al. 1956) predicts that the ratio T/T increases as r2 (where r is the distance from the Sun), if the plasma is expanding radially in a magnetic field that varies as r−2. Given that a temperature anisotropy is a source of free energy for instabilities, it has been argued that the rise of anisotropy should be possible only until the free energy becomes so large that an instability might be triggered, with the consequence of not allowing a further enhancement of anisotropy, and therefore to bound also the increase of free energy.

However, we have shown that, in stable conditions, small perturbations can give rise to transient growth of magnetic fluctuations and of high order moments of the particle distribution function. We argue that the presence of these fluctuations, even far from the stability threshold, can influence the ability of the plasma to further increase the free energy via the anisotropization of the particle distribution function. In fact, the thresholds for marginal stability are derived from the linear theory as the conditions for which the time asymptotic growth rate α becomes null, assuming for the plasma to be in equilibrium in an uniform magnetic field. No consideration is made for the fact that the growth rate is effectively a function of time, and that the magnetic field might not be uniform, due to transient fluctuations. What might happen is therefore that the increase of anisotropy might be bounded even before reaching what is considered the threshold for anisotropy instabilities. In this scenario, a parcel of plasma could experience a "local" marginal stability condition due to a temporary enhanced magnetic fluctuation, and the anisotropy will be reduced in the same way as it is reduced under unstable conditions.

This scenario reconciles the observational contradictions mentioned above with a consistent interpretation of the linear theory. Of course, trying to interpret space observations purely on the base of linear theory might appear a naive approach and a complete and coherent understanding of the process of the solar wind expansion is feasible only via a nonlinear treatment and through computer simulations. However, our results show that a non-modal linear theory could already constitute a strong theoretical ground for data interpretation.

4.2. The Dissipation of Turbulent Fluctuations

In this section, we highlight the particular behavior of the numerical abscissa as a function of wavenumber and point out certain unique features which indicate that non-modal effects may play an important part in the cascade of turbulent fluctuations at short wavelengths.

An approach that is very often used in modeling the dissipation at short wavelength of turbulent fluctuations is to assume that such small fluctuations can be regarded as an ensemble of linear waves, each of them damping according to their linear damping rate. Howes et al. (2008a) have devised a gyrokinetic formalism to study turbulence in a magnetized plasma, in the framework of the Goldreich and Sridhar critical balance assumption (Sridhar & Goldreich 1994). They have shown, by performing nonlinear simulations, that, in the case of small damping rates, the linear damping does not underestimate the rate at which electromagnetic energy is dissipated (Howes et al. 2008b). However, the gyrokinetic cascade model produces an exponential roll-off of the power spectrum for high wavenumbers, instead of the observed power law (Leamon et al. 1998). The fact that a linear damping mechanism cannot reproduce a power law in the spectrum, for high wavenumbers, was previously pointed out by Li et al. (2001), using the full Vlasov–Maxwell linear theory. Our contribution here is to revisit the conclusions reached by Li et al. (2001), by showing that a non-modal approach might be more appropriate to address this problem, which is still largely unsolved.

The equation used by Li et al. (2001) to model the spectral energy W(k) in wavenumber space is of the kind

Equation (10)

where D(k) is a diffusion operator (that can be modeled in such a way to recover a Kolmogorov cascade for large wavelengths), S(k) = S(k0)δ(k −  k0) is a source term that operates at the (small) injection wavenumber k0, and γ(k) is the damping rate given by the Vlasov–Maxwell linear theory. This model is clearly oversimplified in at least three ways. First, it assumes that the nonlinear cascade is a diffusive process in wavenumber space. Second, it assumes that the cascade is isotropic in k (while it is known that anisotropy plays a very important role). Third, it describes the damping process entirely via an ensemble of independent linear waves, all propagating in the same direction, and each damping according to its respective linear damping rate.

Despite the weaknesses of the model, it still represents a good starting point for our purpose of highlighting the importance of transient, non-asymptotic dynamics. The main conclusion of Li et al. (2001) was that "a power-law spectrum in W(k), in the presence of damping, requires γ(k) to be a specific power law." However, at the light of the considerations hitherto drawn, we point out that treating the cascade process as a time-asymptotic problem rather than an initial value problem will lead to possibly different results. In fact, each wavenumber is continuously subject to injection of energy via the nonlinear cascade, and to treat the evolution of a fluctuation, as if it were injected in a remote past, and let free to damp undisturbed is, in our view, an oversimplification of the whole process.

Although it has been shown that some fluid models are able to recover the steepening of the power spectrum purely by a nonlinear mechanism (see, e.g., Galtier & Buchlin 2007), we note that if one assumes for the nonlinear terms in the cascade to act smoothly as a function of k, then it is very likely that at least one of the parameters that control the damping at high frequencies must present some sort of abrupt change in their dependence on k.

This issue necessitates a deeper and separate study, and we do not intend to suggest here a definitive answer. However, it is interesting to show the behavior of the numerical abscissa, as a function of k, for different values of β, and T/T (Figure 6). We recall that the numerical abscissa measures the derivative of the highest amplification that an initial perturbation can undergo, at time t = 0.

Figure 6.

Figure 6. Numerical abscissa η as a function of k for the following parameters: β = 2 (black), 5 (red), and 10 (green), and for T/T = 0.5 (dot-dashed), 1 (solid line), and 2 (dashed). Note that the isotropic case is marked with solid line.

Standard image High-resolution image

Two features are extremely surprising, as shown in Figure 6. First, the numerical abscissa η follows a power law in wavenumber (i.e., a straight line in the log–log plot). Second, it presents an abrupt steepening of its slope, when T/T ≠ 1. As far as we know, this is the only quantity derived from the linear theory that presents these two features. Therefore, it makes sense to suggest that the numerical abscissa should be included in the set of parameters that control the cascade of turbulent fluctuations at short wavelengths. In some sense, this strengthens our argument that a non-modal approach to the problem might be decisive, and we leave further considerations for the future.

5. DISCUSSION

The linear theory of plasma waves is traditionally carried out as a normal-mode analysis. The modal approach, although useful to define the conditions for which an equilibrium might be unstable, is, in general, unable to describe the time evolution of some arbitrary initial perturbation. In particular, it is known that some linear operators can present transient growth of physical quantities, even in stable conditions. We have addressed the stability of a kinetic plasma modeled through a LF model, via a non-modal approach. The following results have been established.

  • 1.  
    Some small initial perturbations can undergo large amplifications, even if applied to a stable equilibrium (Figure 2). If the plasma is close to marginal stability, this might result in a persistent presence of such fluctuations (Figure 3), which is consistent with observations. This is because the (time-asymptotic) damping rate is so small that the decay of an initial perturbation is appreciable only after a time of the order of ∼106 ion gyroperiods. On the other hand, it is quite unrealistic for a plasma to be undisturbed for such large times.
  • 2.  
    Transient growth of magnetic fluctuations or of high order moments of the distribution function are quite possible also for nearly isotropic plasma, i.e., far from the thresholds of linear anisotropy instabilities (Figures 4 and 5).
  • 3.  
    The numerical abscissa, which is the highest growth rate, at time t = 0, for any possible perturbation follows a power law in k and presents an abrupt change in slope for high wavenumbers (Figure 6).

These results are probably more illustrative than predictive, and the goal of the paper was certainly to present an alternative to the traditional use of linear theory, rather than providing a quantitative tool. However, we believe that these results could help to interpret some space observations in a more consistent way, than the modal theory is able to do. We suggest that dissipation scale fluctuations observed in nominally (marginally) stable conditions might be the result of transient growth. We also suggest that the marginal stability concept should be revisited within a non-modal context, and that the prediction of the evolution of small scale fluctuations in a turbulent plasma purely by their time-asymptotic damping rate is an oversimplification of the physics.

Most of the solar wind protons and electrons are observed to be nearly isotropic. This is in contradiction with fluid theories that predict a continuous anisotropization of the distribution function. We suggest that the increase in anisotropy might be inhibited by the presence of ongoing transient effects that could limit the capability of the plasma to increase the free energy.

As for the dissipation of turbulent fluctuations, it is now acknowledged that the understanding of this complex phenomenon lays in the realm of kinetic plasma physics. However, once again the modal approach does not seem to be very useful here and might be misleading. On the other hand, the power-law dependence of the numerical abscissa, and more importantly the presence of a steepening for high wavenumbers, is very interesting. We have pointed out that, since this is the only quantity (derived within linear theory) that presents such features, it should be taken into account.

5.1. Future Directions

The primary intent of this paper has been to show that a non-modal approach could be helpful in reconciling some observational evidence with plasma linear theory and to suggest a more frequent use of such an approach. Historically, the issues that can arise for a non-normal linear operator have been studied in depth in hydrodynamics, but have not captured the same attention in the space plasma community. This of course does not exclude the use of nonlinear theory and computer simulations that will ultimately be the only way to fully understand the processes in action. However, we have shown that the interpretation of simulations should be made with the awareness of the effects produced by non-normal linear operators. Although the LF model used in this paper is superior to MHD, being able to capture some kinetic effects, it will be desirable in the future to develop a non-modal treatment of the full Vlasov–Maxwell set of equations.

Also, as we have shown, the non-modal approach should be framed in a more statistical theoretical framework.

In conclusion, we believe that the understanding and the ability to address transient behavior will be a crucial ingredient in future modeling of space plasma.

This work was partially supported by STFC grant PP/E001424/1 and by "Programme National Soleil Terre" of CNRS.

APPENDIX: FLR-LANDAU FLUID SET OF EQUATIONS

In the context of this paper, it is useful to specify the form of the model in one space dimension, assuming that all the fields only depend on a coordinate ξ, along a direction of the (x, z)-plane making an angle α with the z-axis defined by the uniform ambient magnetic field (of magnitude B0). The total plasma density field is normalized by ρ0, the magnetic field by B0, the velocities by the Alfvén velocity vA = B0/(4πρ0)1/2, the pressures by the parallel ion pressure p0 = pi(0)||, the heat fluxes by p0vA, and the fourth rank moments by p0v2A. The unit of length is the ion inertial length vAi, and time is normalized to ion gyroperiods. The parameter β = 8πp0/B20 measures the ratio of the (parallel) thermal to the magnetic pressure. Velocities without superscripts refer to the ion velocity. The electron velocity ue is given by

Equation (A1)

Equation (A2)

Equation (A3)

where we define bp = cos α bx − sin α bz. We also define $\bar{u}=\sin \alpha \ u_x+\cos \alpha \ u_z$ which, together with $\nabla \cdot u=\partial _\xi \bar{u}$, take the same form when using the electron velocity. When integrated, the divergenceless condition ∇ · b = 0 rewrites cos α bz + sin α bx = cos α, which allows to write $\nabla \cdot \widehat{b}=\nabla \cdot (b/|b|)=-\cos \alpha \partial _\xi |b|/|b|^2$.

The model involves dynamical equations for the ion density ρ and velocity u, the magnetic field components bp and by, together with, for each species r, the gyrotropic parallel and perpendicular pressures pr|| and pr, the heat fluxes qr|| and qr, and the mixed fourth order cumulants $\tilde{r}_{\Vert \perp }^r$. They read

Equation (A4)

Equation (A5)

Equation (A6)

Equation (A7)

Equation (A8)

Equation (A9)

Equation (A10)

Equation (A11)

Equation (A12)

Equation (A13)

Note that the third component of the magnetic field is calculated by imposing the condition ∇ ·  b = 0. The electric field and flux terms entering the above equations are given by

Equation (A14)

Equation (A15)

Equation (A16)

Equation (A17)

Equation (A18)

Equation (A19)

Equation (A20)

Equation (A21)

The operator ${\cal F}^{-1}$ is the inverse Fourier transform. The second terms on the right-hand-side of Equations (A12) and (A13) involve a Hilbert transform with respect to the z-variable (which in Fourier space reads ${\cal H} \hat{f} = i k/|k| \hat{f}$), signature of the Landau damping. The specific form of the expressions for the gyrotropic fourth rank cumulants $\tilde{r}_{\Vert \Vert }^r$ and $\tilde{r}_{\perp \perp }^r$, of the gyroviscous tensor Π, of the non-gyrotropic contributions to the fourth-rank cumulant RrNG, of the transverse components of the fluxes of parallel and perpendicular heat S||rx and Srx as well as of the coefficients ${\cal R}_{\Vert \perp p}^r$ and function ${\cal C}_1$, are computed from the linear kinetic theory and given in Borgogno et al. (2007). They involve functions Γ0 and Γ1 where Γν(b) is the product of exp(−b) by the modified Bessel function Iν(b), with $b\equiv k_\perp ^2 r_L^2/2=\displaystyle {\frac{\beta }{2}{\overline{T}}_\perp ^i k^2 \sin ^2\alpha }$.

Please wait… references are loading.
10.1088/0004-637X/715/1/260