PROTON, ELECTRON, AND ION HEATING IN THE FAST SOLAR WIND FROM NONLINEAR COUPLING BETWEEN ALFVÉNIC AND FAST-MODE TURBULENCE

We model the plasma properties along an open magnetic flux tube rooted in a polar coronal hole. At solar minimum, large unipolar coronal holes are associated with superradially expanding magnetic fields and the acceleration of the high-speed solar wind. Because we only consider a field line along the polar axis of symmetry, we do not need to include the rotational generation of azimuthal magnetic fields (e.g., the Parker spiral effect; see Weber & Davis 1967; Priest & Pneuman 1974) or other geometrical effects of streamer-like flux tube curvature (Li et al. 2011). We do not distinguish between dense polar plumes and the more tenuous interplume regions between them. The radial dependence of plasma parameters is described as a function of either the heliocentric distance (r) or the height above the solar photosphere (z = r − R_☉).

To specify the radial variation of the time-steady magnetic field strength B₀, mass density ρ₀, and solar wind outflow speed u₀, we used the empirical description of Cranmer & van Ballegooijen (2005). This model combined a broad range of observational constraints with a two-dimensional magnetostatic model of the expansion of thin photospheric flux tubes into a supergranular network canopy. At r = 1 AU in this model, the solar wind outflow speed u₀ is 781 km s⁻¹ and the proton density n_p is 2.56 cm⁻³. This model also specifies the Alfvén speed V_A = B₀/(4πρ₀)^1/2, which decreases from a maximum value of 2890 km s⁻¹ at r = 1.53 R_☉ down to 31 km s⁻¹ at 1 AU. There is a local minimum in V_A at r ≈ 1.02 R_☉ that is the result of the assumed shape of network "funnels" that expand superradially into the corona.

We also need to know the plasma temperature T in order to determine the relative importance of gas pressure versus magnetic pressure. Despite observational evidence for different particle species having different temperatures (and departures from Maxwellian velocity distributions), we generally assume that the majority proton–electron magnetofluid is close enough to thermal equilibrium that strong plasma microinstabilities are not excited (e.g., Gary 1991; Marsch 2006). Thus, we specify a one-fluid temperature T that is assumed to be equal to both the proton temperature T_p and the electron temperature T_e, and we assume temperature isotropy (T_∥ ≈ T_⊥) for both species.

We used the polar coronal hole model of Cranmer et al. (2007) as a starting point to describe T(r), but this model was modified in two ways. First, we moved the sharp transition region (TR) down from a height z of 0.01 to 0.003 solar radii (R_☉) to better match the conditions of semi-empirical models (e.g., Fontenla et al. 1990; Cranmer & van Ballegooijen 2005; Avrett & Loeser 2008). Thus, in the adopted model, at z = 0.01 R_☉ the temperature has risen to 0.48 MK, and it continues to rise to a maximum value of 1.36 MK at z = 0.89 R_☉. We also increased the temperature slightly at distances greater than ∼0.2 AU in order to better agree with the mean of the in situ T_p and T_e measurements of Cranmer et al. (2009). At r = 1 AU, T = 0.17 MK and it declines as T∝r^−0.6. The one-fluid sound speed c_s is defined as c²_s = γk_BT/m_H, where γ = 5/3 is the monatomic ratio of specific heats, k_B is Boltzmann's constant, and m_H is the hydrogen atomic mass.

Figure 1 shows the radial dependence of a selection of the background plasma properties defined above. It also shows the dimensionless plasma beta parameter, which is usually defined as the ratio of gas pressure to magnetic pressure, with

$$\begin{equation} \beta _{0} \, = \, \frac{P_{\rm gas}}{P_{\rm mag}} \, = \, \frac{2}{\gamma } \left(\frac{c_s}{V_{\rm A}} \right)^{2}. \end{equation} \tag{ 1 }$$

However, we will often use a simpler dimensionless parameter β given by

$$\begin{equation} \beta \, = \, \left(\frac{c_s}{V_{\rm A}} \right)^{2} \, = \, \frac{\gamma \beta _{0}}{2}, \end{equation} \tag{ 2 }$$

where β and β₀ differ only by a factor of 1.2 when γ = 5/3. The range of heights shown in Figure 1 extends down into the solar chromosphere, but the wave models discussed below start at a lower boundary condition in the low corona; i.e., they specify the wave and turbulence properties only for z ⩾ 0.01 R_☉.

Zoom In Zoom Out Reset image size

In this section we briefly summarize the dispersion properties of linear MHD waves (i.e., phase and group speeds for the Alfvén mode and the fast and slow magnetosonic modes) and the partitioning between fluctuations in kinetic, magnetic, and thermal energy. In Sections 2.3–2.4, we assume that all three types of MHD waves are present, and we vary their relative strengths arbitrarily in order to match the observations.

The phase speed V_ph = ω/k is defined in terms of the frequency ω and the magnitude of the wavenumber k. In general, V_ph is a function of the Alfvén speed, the sound speed, and the angle θ between the background field direction and the wavevector k. We follow the standard convention of defining a Cartesian coordinate system with the background magnetic field along the z-axis and the k vector having components only in the x–z plane. Also, for now we express ω and k in the frame comoving with the solar wind. For Alfvén waves,

$$\begin{equation} V_{\rm ph}^{2} \, = \, V_{\rm A}^{2} \cos ^{2} \theta \end{equation} \tag{ 3 }$$

and for the magnetosonic modes,

$$\begin{equation} V_{\rm ph}^{2} \, = \, \frac{V_{\rm A}^{2} + c_{s}^{2}}{2} \left(1 \pm \Sigma \right) \end{equation} \tag{ 4 }$$

applies with the upper sign corresponding to the fast mode and the lower sign corresponding to the slow mode, and with

$$\begin{equation} \Sigma \, = \, \sqrt{1 - \sigma \cos ^{2} \theta }, \,\,\,\, \sigma \, = \, \frac{4\beta }{(1+\beta)^{2}} \end{equation} \tag{ 5 }$$

(see, e.g., Whang 1997; Goedbloed & Poedts 2004). In Section 2.3, we also need to know the component of an MHD wave's group velocity in the direction parallel to the background magnetic field. We call this quantity V_gz, and for the Alfvén mode it is identically equal to V_A no matter the value of θ. For the fast and slow modes,

$$\begin{equation} V_{\rm gz} \, = \, V_{\rm ph} \cos \theta \left(1 \mp \frac{\sigma \sin ^{2}\theta }{2 \Sigma \, (1 \pm \Sigma)} \right). \end{equation} \tag{ 6 }$$

MHD waves excite oscillations in the plasma parameters. We denote the root-mean-square (rms) fluctuation amplitudes in velocity as v_x, v_y, v_z, in magnetic field as B_x, B_y, B_z, and in density as δρ. We ignore fluctuations in the electric field because their contribution to the total energy density tends to be negligible when V_A ≪ c. The kinetic, magnetic, and thermal energy densities associated with each type of fluctuation are given as

$$\begin{equation} K_{i} = \frac{\rho _{0} v_{i}^2}{2}, \,\,\,\, M_{i} = \frac{B_{i}^2}{8\pi }, \,\,\,\, \Theta = \beta \frac{B_{0}^{2}}{8\pi } \left(\frac{\delta \rho }{\rho _0} \right)^{2}, \end{equation} \tag{ 7 }$$

respectively, with i = x, y, z. For linear Alfvén waves, the total energy density U_A is divided equally between transverse kinetic and magnetic fluctuations along the y-axis, with U_A = K_y + M_y and

$$\begin{equation} \frac{K_y}{U_{\rm A}} \, = \, \frac{M_y}{U_{\rm A}} \, = \, \frac{1}{2}. \end{equation} \tag{ 8 }$$

For fast- and slow-mode waves,

$$\begin{equation} U_{\rm F,S} \, = \, K_{x} + K_{z} + M_{x} + M_{z} + \Theta \end{equation} \tag{ 9 }$$

and we follow Whang (1997) in expressing the partition fractions as follows:

$$\begin{equation} \frac{K_x}{U_{\rm F,S}} \, = \, f_{n} \sin ^{2}\theta + f_{t} \cos ^{2}\theta, \,\,\,\, \frac{K_z}{U_{\rm F,S}} \, = \, f_{n} \cos ^{2}\theta + f_{t} \sin ^{2}\theta \end{equation} \tag{ 10 }$$

$$\begin{equation} f_{n} \, = \, \frac{V_{\rm ph}^{2} \big(V_{\rm ph}^{2} - V_{\rm A}^{2}\big)}{c_{s}^{2} \Delta }, \,\,\,\, f_{t} \, = \, \frac{V_{\rm A}^{2} \big(V_{\rm ph}^{2} - c_{s}^{2}\big) \cos ^{2}\theta }{V_{\rm ph}^{2} \Delta } \end{equation} \tag{ 11 }$$

$$\begin{equation} \frac{M_x}{U_{\rm F,S}} \, = \, \frac{\big(V_{\rm ph}^{2} - c_{s}^{2}\big) \cos ^{2}\theta }{\Delta }, \,\,\,\, \frac{M_z}{U_{\rm F,S}} \, = \, \frac{\big(V_{\rm ph}^{2} - c_{s}^{2}\big) \sin ^{2}\theta }{\Delta } \end{equation} \tag{ 12 }$$

$$\begin{equation} \frac{\Theta }{U_{\rm F,S}} \, = \, \frac{V_{\rm ph}^{2} - V_{\rm A}^{2}}{\Delta }, \end{equation} \tag{ 13 }$$

where Δ = 4V²_ph − 2V_A² − 2c²_s. The fast and slow velocity fluctuations (K_x + K_z) always occupy exactly half of the total energy density, and the combination of magnetic and thermal fluctuations (M_x + M_z + Θ) takes up the other half.

The energy partition fractions given above are familiar components of plasma physics and MHD textbooks (e.g., Stix 1992; Goedbloed & Poedts 2004). However, it is difficult to see intuitively how these fractions vary throughout the heliosphere from Equations (10)–(13) alone. Thus, in Figure 2 we provide a schematic illustration of the energy partitioning for fast-mode waves. The three columns indicate the variation from low (β ≪ 1) to medium (β = 1) and high (β ≫ 1) beta plasmas. The three rows show the results for purely parallel propagation (θ = 0), an isotropic distribution of wavenumber vectors (see below), and purely perpendicular propagation (θ = π/2). In general, all five terms on the right-hand side of Equation (9) are nonzero, but fractions less than ∼1% are not shown in Figure 2. This diagram can be transformed to show the properties of slow-mode waves by replacing β with 1/β and interchanging the x and z subscripts with one another.

Zoom In Zoom Out Reset image size

In order to determine how the total energy density of a given wave mode evolves with heliocentric distance, we solve equations of wave action conservation that contain multiple sources of wave damping. There have been many discussions of energy conservation for both pure acoustic waves and incompressible Alfvén waves (e.g., Dewar 1970; Isenberg & Hollweg 1982; Velli 1993; Tu & Marsch 1995; Verdini & Velli 2007; Sokolov et al. 2009), but general derivations that can also be applied to fast- and slow-mode waves (for arbitrary θ) are less frequently seen. We utilize the results of Jacques (1977) to write the damped wave action conservation equation as

$$\begin{equation} \frac{\partial }{\partial t} \left(\frac{U_m}{\Omega } \right) + \frac{1}{A_0} \frac{\partial }{\partial r} \left(\frac{\langle u_{0} + V_{{\rm gz,} m} \rangle A_{0} U_m}{\Omega } \right) \, = \, -\frac{Q_m}{\Omega }, \end{equation} \tag{ 14 }$$

where the subscript m can be replaced by A, F, or S for the relevant mode, A₀ is the cross sectional area of the flux tube (i.e., A₀∝1/B₀), and Q_m is the total dissipation rate for the mode in question. The dimensionless factor that takes account of the "stretching" effect of wavelengths in an accelerating reference frame is

$$\begin{equation} \Omega \, = \, \left\langle \frac{V_{{\rm ph,} m}}{u_{0} \cos \theta + V_{{\rm ph,} m}} \right\rangle \end{equation} \tag{ 15 }$$

and the angle brackets denote a weighted average over all angles,

$$\begin{equation} \langle f \rangle \, = \, \frac{\int d\theta \, \sin \theta \, f(\theta)}{\int d\theta \, \sin \theta }, \end{equation} \tag{ 16 }$$

where here we consider outward propagating waves with 0 < θ < π/2. The factor of cos θ in Equation (15) comes from the difference between the wave frequency in the Sun's reference frame (ω₀) and the comoving-frame frequency (ω = ω₀ − k · u₀) that appears in the definition of the wave action; see Section III of Jacques (1977).¹ Equation (14) implicitly assumes that ω₀ remains constant, but it does not require the specification of any given value of ω₀.

Our use of weighted averages over θ is derived from the assumption that wave power is distributed isotropically in three-dimensional k-space. In Appendix A we discuss the motivations for assuming such an isotropic distribution of wavenumber vectors (specifically for the fast-mode waves). For Alfvén waves, this assumption has no impact on solving Equation (14), since the arguments of both angle-bracketed quantities given above are independent of θ (see also Hollweg 1974). Thus, one obtains the same result for Alfvén waves whether one assumes a single value of θ or the isotropic distribution. For fast- and slow-mode MHD waves, some quantities depend strongly on θ and others do not. For example, slow-mode waves in low-beta plasmas have values of the angle-dependent quantity u₀ + V_gz,m that are always nearly equal to u₀ + c_s. Figure 1 shows the radial dependence of 〈u₀ + V_{gz, F}〉 for the isotropic distribution of fast-mode waves.

We solve Equation (14) for the energy densities of the three MHD modes (U_A, U_F, U_S), and we compute the dispersion and energy partition properties of all three wave types as given in Section 2.2. At this stage, we neglect couplings between multiple modes and other nonlinear effects. This is an approximation that is likely to break down wherever the wave amplitudes become large (e.g., Chin & Wentzel 1972; Wentzel 1974; Goldstein 1978; Lacombe & Mangeney 1980; Poedts et al. 1998; Vasquez & Hollweg 1999; Del Zanna et al. 2001; Gogoberidze et al. 2007). In Section 4, we discuss the likelihood of rapid coupling between the high-wavenumber tails of the Alfvén and fast-mode power spectra. However, we continue to assume that the total energy densities are given by the solution of the individual transport equations.

To specify the dissipation rates Q_m, we include both linear collisional effects (e.g., viscosity, thermal conductivity, and electrical resistivity) for all three modes and nonlinear turbulent damping for the Alfvén and fast mode. Thus, we use

$$\begin{equation} Q_{\rm A} = \tilde{Q}_{\rm A} + 2 \gamma _{\rm A} U_{\rm A}, \,\,\,\, Q_{\rm F} = \tilde{Q}_{\rm F} + 2 \gamma _{\rm F} U_{\rm F}, \,\,\,\, Q_{\rm S} = 2 \gamma _{\rm S} U_{\rm S}. \end{equation} \tag{ 17 }$$

We give the amplitude damping rates γ_m, which include an approximation for the transition from strongly collisional to collisionless regimes, in Appendix B. The turbulent damping rates $\tilde{Q}_{\rm A}$ and $\tilde{Q}_{\rm F}$ are described in more detail below. In general, these rates depend on the parallel and perpendicular components of the wavenumber (k_∥, k_⊥). For the purposes of evaluating these rates in the global wave transport equations, we assumed that k_⊥ = 1/λ_⊥ for all three modes, where λ_⊥ is the turbulent correlation length described below. For the fast and slow modes, our assumption of an isotropic distribution of wavenumbers is consistent with also assuming k_∥ = k_⊥. For the Alfvén mode, we found that γ_A never depended on the assumed value of k_∥ at all, but for completeness we used the critical balance condition (introduced in Appendix A) to specify k_∥.

We adopt phenomenological forms for the turbulent dissipation rates that are equivalent to the total energy fluxes that cascade from large to small scales. Thus, $\tilde{Q}_{\rm A}$ and $\tilde{Q}_{\rm F}$ are constrained only by the properties of fluctuations at the largest scales, and they do not specify the exact kinetic means of dissipation once the energy reaches the smallest scales (but see, however, Section 5). Dimensionally, these are similar to the rate of cascading energy flux derived by von Kármán & Howarth (1938) for isotropic hydrodynamic turbulence. For the nonlinear dissipation of Alfvénic fluctuations, we use

$$\begin{equation} \tilde{Q}_{\rm A} \, = \, \rho _{0} \, \tilde{\alpha }_{\rm A} {\cal E}_{\rm turb} \, \frac{Z_{-}^{2} Z_{+} + Z_{+}^{2} Z_{-}}{4 \lambda _{\perp }} \end{equation} \tag{ 18 }$$

(see also Hossain et al. 1995; Zhou & Matthaeus 1990a; Matthaeus et al. 1999; Dmitruk et al. 2001, 2002; Breech et al. 2008). For the fast-mode waves, we use

$$\begin{equation} \tilde{Q}_{\rm F} \, = \, \rho _{0} \, \tilde{\alpha }_{\rm F} \, \frac{\big(v_{x}^{2} + v_{z}^{2}\big)^2}{V_{\rm A} \lambda _{\perp }}, \end{equation} \tag{ 19 }$$

where the quantity (v²_x + v_z²) collects together the total kinetic energy in fast-mode velocity fluctuations (Chandran 2005; Suzuki et al. 2007). Many of the terms introduced in Equations (18) and (19) are defined throughout the remainder of this subsection.

Equation (18) depends on the magnitudes of the Elsasser (1950) variables, Z_± = v_y ± B_y/(4πρ₀)^1/2, which specify the power in outward (Z₋) and inward (Z₊) propagating Alfvénic fluctuations. Alfvénic turbulent heating occurs only when there is energy in both modes. In practice, we compute an effective reflection coefficient ${\cal R} = |Z_{+}| / |Z_{-}|$ whose magnitude is always less than unity, and thus we express the Elsasser variables in terms of the Alfvénic energy density as

$$\begin{equation} Z_{-} \, = \, \sqrt{\frac{4 U_{\rm A}}{\rho _{0} (1 + {\cal R}^{2})}}, \,\,\,\, Z_{+} \, = \, {\cal R} Z_{-}. \end{equation} \tag{ 20 }$$

An accurate solution for Z_± requires the integration of non-WKB equations of Alfvén wave reflection (e.g., Heinemann & Olbert 1980; Verdini & Velli 2007). However, our assumption that the total power U_A varies in accord with straightforward wave action conservation has been shown to be reasonable, even in environments where ${\cal R}$ is not small such as the chromosphere (van Ballegooijen et al. 2011) and interplanetary space (Zank et al. 1996; Cranmer & van Ballegooijen 2005).

We estimate the reflection coefficient ${\cal R}$ using a modification of the low-frequency approximation of Chandran & Hollweg (2009). Specifically, we examine the magnitudes of terms in the transport equation for the inward Elsasser variable,

$$\begin{equation} \frac{\partial Z_{+}}{\partial t} + (u_{0} - V_{\rm A}) \frac{\partial Z_{+}}{\partial r} \, = \, (u_{0} + V_{\rm A}) \left(\frac{Z_{+}}{4 H_{\rm D}} + \frac{Z_{-}}{2 H_{\rm A}} \right) - \frac{Z_{+} Z_{-}}{2 \lambda _{\perp }}, \end{equation} \tag{ 21 }$$

where

$$\begin{equation} H_{\rm A} = \frac{V_{\rm A}}{\partial V_{\rm A} / \partial r}, \,\,\,\, H_{\rm D} = \frac{\rho _0}{\partial \rho _{0} / \partial r}. \end{equation} \tag{ 22 }$$

Chandran & Hollweg (2009) neglected both terms on the left-hand side of Equation (21) as well as the term containing H_D, and thus were able to solve for Z₊ straightforwardly. However, in cases of strong reflection, the term containing H_D may have a magnitude comparable to the other dominant terms. Thus, we keep all three terms on the right-hand side and solve for

$$\begin{equation} {\cal R} \, \approx \, \frac{2 h / | H_{\rm A} |}{1 + (h/ | H_{\rm D} |)}, \end{equation} \tag{ 23 }$$

where

$$\begin{equation} h \, = \, \frac{\lambda _{\perp } (u_{0} + V_{\rm A})}{2 Z_{-}}. \end{equation} \tag{ 24 }$$

Equation (22) of Chandran & Hollweg (2009) is recovered in the limit of h ≪ |H_D|, with ${\cal R} \approx 2h/ |H_{\rm A}|$. In the case of purely linear reflection, Cranmer (2010) found that the most accurate local estimates for ${\cal R}$ were obtained when H_A was replaced with the positive-definite quantity

$$\begin{equation} \tilde{H}_{\rm A} \, = \, V_{\rm A} t_{\rm ref} \, = \, (r+R_{\odot }) \left(1 - \frac{R_{\odot }}{r} \right). \end{equation} \tag{ 25 }$$

We used $\tilde{H}_{\rm A}$ instead of H_A in Equations (23) and (24) to compute ${\cal R}$.

The definitions of the turbulent dissipation rates contain the perpendicular length scale λ_⊥, which is an effective transverse correlation length of the turbulence for the largest "outer-scale" eddies. For simplicity we use the same correlation length for both the Alfvénic and fast-mode fluctuations, but this may not be universally valid (e.g., Suzuki et al. 2007). In previous papers, we assumed that λ_⊥ scales with the transverse width of the magnetic flux tube; i.e., that λ_⊥∝B^−1/2₀ (Hollweg 1986). Here we describe the evolution of the transverse correlation length λ_⊥ with the following transport equation,

$$\begin{equation} \frac{\partial \lambda _{\perp }}{\partial r} \, = \, \frac{\lambda _{\perp }}{2A_0} \frac{\partial A_0}{\partial r} + \frac{\tilde{\beta }_{\rm A}}{u_{0} + V_{\rm A}} \left(\frac{Z_{-}^{2} Z_{+} + Z_{+}^{2} Z_{-}}{Z_{-}^{2} + Z_{+}^{2}} \right), \end{equation} \tag{ 26 }$$

where $\tilde{\beta }_{\rm A}$ is a dimensionless constant that is often assumed to be equal to $\tilde{\alpha }_{\rm A}/2$ (e.g., Hossain et al. 1995). The first term on the right-hand side of Equation (26) drives the correlation length to expand linearly with the perpendicular flux-tube cross section (Hollweg 1986). The second term takes account of the nonlinear coupling between the fluctuations and the background plasma properties. It is given in a form suggested initially by Matthaeus et al. (1994) and later generalized to nonzero cross-helicity turbulence by Breech et al. (2008) and others. Our transport equation attempts to bridge together the effects of the two terms. In the lower solar atmosphere (between the photosphere and the chosen lower boundary of 0.01 R_☉ for the wave transport models) we assumed that the first term in Equation (26) is dominant, and thus λ_⊥∝A^1/2₀.

The turbulent dissipation rates also depend on dimensionless Kolmogorov-type constants $\tilde{\alpha }_{\rm A}$ and $\tilde{\alpha }_{\rm F}$ that are often assumed to have values of order unity. For example, Hossain et al. (1995) and Breech et al. (2009) found that $\tilde{\alpha }_{\rm A} \approx 0.5$ gives rise to dissipation rates that agree well with both numerical simulations and heliospheric observations. In our case, we used this value as a starting point, but we also varied $\tilde{\alpha }_{\rm A}$ as a free parameter in order to produce the best match to the well-constrained Alfvénic fluctuations. On the other hand, the properties of heliospheric fast-mode turbulence are not known nearly as well as the Alfvén wave turbulence. We thus relied on the independent wave-kinetic simulations of P. Pongkitiwanichakul & B. D. G. Chandran (2012, in preparation) to fix $\tilde{\alpha }_{\rm F}$ at a value of 2.3.

The Alfvénic cascade rate contains an efficiency factor ${\cal E}_{\rm turb}$ that attempts to account for regions where the turbulent cascade may not have time to develop before the fluctuations are carried away by the wind. Cranmer et al. (2007) estimated this efficiency factor to scale as

$$\begin{equation} {\cal E}_{\rm turb} \, = \, \frac{1}{1 + (t_{\rm eddy} / t_{\rm ref})}, \end{equation} \tag{ 27 }$$

where the two timescales above are t_eddy, a nonlinear eddy cascade time, and t_ref, a timescale for large-scale Alfvén wave reflection (see also Dmitruk & Matthaeus 2003; Oughton et al. 2006). The reflection time is often defined as t_ref = 1/|∇ · V_A|, but we solved Equation (25) for t_ref in order to remain consistent with the adopted model for ${\cal R}$. The eddy cascade time is given by

$$\begin{equation} t_{\rm eddy} \, = \, \frac{\lambda _{\perp } \sqrt{3\pi }}{(1 + M_{\rm A}) \, v_{y}}, \end{equation} \tag{ 28 }$$

where the Alfvén Mach number M_A = u₀/V_A and the numerical factor of 3π comes from the normalization of an assumed shape of the turbulence spectrum (see Appendix C of Cranmer & van Ballegooijen 2005). The two limiting cases of ${\cal E}_{\rm turb} \ll 1$ and ${\cal E}_{\rm turb} \approx 1$ are roughly equivalent to the "weak" and "strong" cascade phenomenologies discussed in Section 3, but they are not precisely the same.

We solved the transport equations given in Section 2.3 by numerically integrating upward from a specified set of lower boundary conditions at z = 0.01 R_☉ and assuming time-steady conditions (i.e., ∂U_m/∂t = 0). We used a logarithmic grid of 500 radial zones in z that expands out to a maximum distance of 860 R_☉ ≈ 4 AU. The transport equations were solved with straightforward first-order Euler steps. The values of the Elsasser variables Z_± in each zone were determined by iteration, since Equations (20) and (23) do not give a simple closed-form solution for Z₊ and Z₋ by themselves.

There are a number of free parameters in this model whose values were not easily obtained from either theoretical calculations or observations. In addition to the lower boundary conditions on the wave energy densities U_A, U_F, and U_S, there is also the lower boundary condition on the correlation length λ_⊥ and the values of the two von Kármán constants $\tilde{\alpha }_{\rm A}$ and $\tilde{\beta }_{\rm A}$. Initially, we varied these six parameters randomly in order to build up a large Monte Carlo ensemble of trial solutions. For each model, we synthesized the radial variation of observable plasma fluctuations such as the rms parallel and perpendicular fluctuation speeds

$$\begin{equation} v_{\parallel } \, = \, v_{z}, \,\,\,\,\,\,\, v_{\perp } \, = \, \big(v_{x}^{2} + v_{y}^{2}\big)^{1/2}, \end{equation} \tag{ 29 }$$

the Elsasser variables Z_±, and the rms fractional density fluctuation amplitude δρ/ρ₀. The velocity amplitudes v_∥ and v_⊥ contain contributions from all three MHD wave types. In nearly all models produced here, v_∥ is dominated by the fast mode and v_⊥ is dominated by the Alfvén mode. Observations of these quantities are discussed below.

There was no single set of parameter values that gave rise to perfect agreement between all of the synthesized and observed fluctuation quantities. This is not surprising, since the models are certainly incomplete and there are significant uncertainties in the observations and their interpretation. Also, even though we aimed to restrict ourselves to measurements made in "quiet" high-latitude fast-wind streams, sometimes only low-latitude data were available. Thus, in Table 1 we give a set of optimized parameters that were chosen because they produce adequate agreement with the full set of observed quantities. There were other combinations of the six parameters that gave better agreement on any single observation, but in most of these cases the agreement became worse for other observations. Although the ratio of the two von Kármán constants $\tilde{\alpha }_{\rm A}/\tilde{\beta }_{\rm A}$ was allowed to vary freely, the optimal value was nonetheless found to be close to the commonly used value of 2 (Hossain et al. 1995). The best photospheric value of λ_⊥ ≈ 120 km is intermediate between the values of 75 km (Cranmer et al. 2007) and 300 km (Cranmer & van Ballegooijen 2005) found from earlier models.

Table 1. Standard Model Parameters for Coronal Hole MHD Wave Transport

Parameter	Value
$\tilde{\alpha }_{\rm A}$	0.60
$\tilde{\beta }_{\rm A}$	0.31
$\tilde{\alpha }_{\rm F}$	2.3
(UA/ρ0)1/2 (at z = 0.01 R☉)	29.0 km s−1
(UF/ρ0)1/2 (at z = 0.01 R☉)	24.3 km s−1
(US/ρ0)1/2 (at z = 0.01 R☉)	9.17 km s−1
λ⊥ (at photosphere)	120 km

Download table as:  ASCIITypeset image

Figure 3 shows the comparison between synthesized and observed fluctuation quantities for the model parameters given in Table 1. The observational constraints on v_⊥ at z ≲ 0.1 R_☉ are a combination of the off-limb nonthermal emission line widths given by Banerjee et al. (1998) and Landi & Cranmer (2009). The observations shown between 0.3 and 1 R_☉ are from Esser et al. (1999). At larger heights, v_⊥ becomes approximately equal to Z₋/2, so we truncate the v_⊥ curve in favor of showing the radial dependence of Z₊ and Z₋ more clearly. The latter are compared directly with high-speed wind data from Helios and Ulysses (Bavassano et al. 2000). Observations of longitudinal velocity fluctuations are more difficult to find, and we show only the on-disk nonthermal line width velocities of Chae et al. (1998) as a way to compare with the modeled values of v_∥.

Zoom In Zoom Out Reset image size

Figure 3(a) shows how the modeled density fluctuation amplitude δρ/ρ₀ is dominated by slow-mode waves in the low corona (z ≲ 0.1 R_☉) and by fast-mode waves in the extended corona and solar wind (z ≳ 1 R_☉). The low-corona observations are drawn as an approximate boundary region around the polar plume data given by Ofman et al. (1999). The intermediate data point at z = 4 R_☉ is an empirical value of δρ/ρ₀ estimated from radio sounding data (Coles & Harmon 1989; Spangler 2002; Harmon & Coles 2005; Chandran et al. 2009), but it is still unclear what fraction of the measured density fluctuations is due to anything even close to ideal MHD waves. At larger distances, we show approximate ranges of density fluctuations as reported by Marsch & Tu (1990) (blue rectangles at z < 200 R_☉), Tu & Marsch (1994) (open rectangle), and Issautier et al. (1998) (blue rectangle at z > 300 R_☉).

Figure 3(b) compares the result of solving Equation (26) for λ_⊥ with the simpler approximation of λ_⊥∝B^−1/2₀. The plot also shows fast-wind estimates of λ_⊥ between 1.4 and 5 AU from Ulysses (Breech et al. 2008). Figure 3(b) also compares the total heating rate Q_tot = Q_A + Q_F + Q_S with observational constraints and with the modeled coronal heating rate from Cranmer et al. (2007). The shaded area between 0.2 and 5 R_☉ is an envelope surrounding a collection of empirical and theoretical heating curves from Wang (1994), Hansteen & Leer (1995), and Allen et al. (1998). These rates illustrate what is needed to produce the observed coronal heating and solar wind acceleration. The area shown at larger distances (z > 60 R_☉) is a representation of the range of total (proton and electron) empirical heating rates estimated by Cranmer et al. (2009). Note that the turbulent heating rates $\tilde{Q}_{\rm A}$ and $\tilde{Q}_{\rm F}$ dominate the total heating rate, with approximately 70% of the total coming from $\tilde{Q}_{\rm A}$ and 20% from $\tilde{Q}_{\rm F}$. Less than 10% of Q_tot comes from the linear damping terms.

There are additional measurement techniques that may be used to further constrain the model parameters, and in future work we will incorporate as many of these as possible. For example, Hollweg et al. (2010) argued that radio measurements of Faraday rotation fluctuations may put unique empirical constraints on the value of λ_⊥ in the corona. Also, Sahraoui et al. (2010) used multi-spacecraft data to tease out new details of the wavenumber anisotropy of MHD fluctuations, which may lead to better limits on, e.g., v_∥/v_⊥ in the heliosphere. Unfortunately, the vast majority of these measurements have been made for the slow solar wind and not the much less structured fast wind associated with polar coronal holes. Nearer to the Sun, Kitagawa et al. (2010) used the dispersive and energy partition properties of thin-tube MHD waves to diagnose the presence and strengths of various modes in active regions. These techniques may be useful in open-field regions as well.

Although we did not include any explicit multi-mode coupling in the transport equations of Section 2.3, there is some feedback between the modes. For example, the correlation length λ_⊥ is used in both the Alfvénic and fast-mode turbulent heating expressions, and it is also used to set the wavenumbers k_∥ and k_⊥ in the linear dissipation rates γ_m. Thus, the choice of the lower boundary condition on λ_⊥ can have a significant impact on the radial evolution of all three wave types. Figure 4 illustrates this by varying the photospheric value of λ_⊥ between 30 and 300 km and using the other standard parameters from Table 1. The integrated energy densities are plotted in velocity units as (U_m/ρ₀)^1/2. The power in the Alfvén waves changes by only a small amount because the damping is never a strong contributor to the U_A transport equation. However, damping is a major effect for the fast and slow modes, and thus small changes in the damping rate's normalization can have large relative impacts on the resulting energy densities.

Zoom In Zoom Out Reset image size

Figure 4 shows that, no matter the choice of normalization for λ_⊥, it seems unlikely for the slow-mode waves to have significantly large amplitudes anywhere but in the lowest few tenths of a solar radius. This appears to be consistent with models of slow-mode shock formation and dissipation in polar plumes (Cuntz & Suess 2001). Therefore, in the remainder of this paper, our models of turbulence in the fast solar wind ignore the slow-mode waves altogether. We also note that Figure 4 suggests that the actual fast-mode wave properties in the high-speed solar wind may be more highly variable than the Alfvén wave properties. Our use of a "standard" model for the fast-mode waves (using the parameters given in Table 1) is thus presented as an example case and not a definitive prediction.

We model the MHD fluctuations as time-steady Fourier distributions of wave power in three-dimensional wavenumber space. Although additional information about the physics of turbulence can be found in more complex statistical measures of the system (e.g., higher-order structure functions), we limit ourselves to describing the power spectrum because that is the basic quantity needed to compute the quasilinear particle heating rates.

Because of the simplified flux-tube geometry discussed in Section 2.1, we assume the background magnetic field is parallel to the bulk flow velocity, and thus the system has only one preferred spatial direction (see, however, Narita et al. 2010). The random turbulent motions create a statistical equivalence between the x- and y-directions transverse to the background field, so that we can describe the power spectra as two-dimensional functions of k_∥ and k_⊥ only. By convention, we define the full three-dimensional power spectrum E_m in effective velocity-squared units; i.e., when integrated over the full volume of wavenumber space, the spectrum gives the fluctuation energy density per unit mass, or

$$\begin{equation} \frac{U_m}{\rho _0} \, = \, \int d^{3} {\bf k} \,\, E_{m} ({\bf k}). \end{equation} \tag{ 30 }$$

In Appendix A, we review some of the basic physical processes that determine the shape of the spectrum for Alfvénic (m = A) and fast-mode (m = F) fluctuations.

We describe the driven turbulent cascade as a combination of advection and diffusion in wavenumber space. At first, it may appear that a smooth and continuous description of the spectral "spreading" of a cascade ignores too much of the inherently stochastic and nonlocal nature of turbulence. However, Chandrasekhar (1943) showed that such a model can be made to capture the essential statistics of a large ensemble of random-walk-like (i.e., Brownian) processes. Specific models of turbulent wavenumber transport using diffusion or advection equations include those of Pao (1965), Leith (1967), Tu et al. (1984), Tu (1988), Zhou & Matthaeus (1990b), Miller et al. (1996), Stawicki et al. (2001), Chandran (2008b), Matthaeus et al. (2009), Jiang et al. (2009), and Galtier & Buchlin (2010). For the cascade of Alfvénic fluctuations, we generally follow the approach taken by Cranmer & van Ballegooijen (2003). The general forms of these equations are given as

$$\begin{eqnarray} \frac{\partial E_{\rm A}}{\partial t} \, &=& \, \frac{1}{k_{\perp }} \frac{\partial }{\partial k_{\perp }} \left\lbrace D_{{\rm A}\perp } \left[ \frac{\alpha _{\perp }}{k_{\perp }} \frac{\partial }{\partial k_{\perp }} \left(k_{\perp }^{2} E_{\rm A} \right) - \mu _{\perp } E_{\rm A} \right] \right\rbrace\nonumber\\ && + \alpha _{\parallel } \frac{\partial }{\partial k_{\parallel }} \left(D_{{\rm A}\parallel } \frac{\partial E_{\rm A}}{\partial k_{\parallel }} \right) + S_{\rm A} - 2 \gamma _{\rm A} E_{\rm A} + C_{\rm AF} \end{eqnarray} \tag{ 31 }$$

$$\begin{equation} \frac{\partial E_{\rm F}}{\partial t} \, = \, \frac{\alpha _{\rm F}}{k^2} \frac{\partial }{\partial k} \left(k^{2} D_{\rm F} \frac{\partial E_{\rm F}}{\partial k} \right) + S_{\rm F} - 2 \gamma _{\rm F} E_{\rm F} - C_{\rm AF} \end{equation} \tag{ 32 }$$

and the terms on the right-hand sides of Equations (31) and (32) are defined throughout the remainder of this subsection. The mode coupling term C_AF is described further in Section 4, and the dissipation rates γ_A and γ_F are described in Section 5.

The perpendicular Alfvénic cascade is described by the first term on the right-hand side of Equation (31), and we assume an arbitrary linear combination of advection and diffusion. Cranmer & van Ballegooijen (2003) found that many key properties of the turbulence do not depend on whether the cascade is modeled as advection, diffusion, or both, so we retain all terms for maximum generality. For both the parallel Alfvénic spectral transport and the isotropic fast-mode transport, a more standard diffusion coefficient is assumed. The dimensionless multipliers to the E_A diffusion coefficients are denoted α_⊥ and α_∥, to correspond roughly to $\tilde{\alpha }_{\rm A}$ in Equation (18), and the dimensionless multiplier for the wavenumber advection coefficient is denoted μ_⊥.

For the Alfvénic cascade, the overall behavior of wavenumber transport in the perpendicular and parallel directions is specified by the diffusion-like coefficients

$$\begin{equation} D_{{\rm A}\perp } \, = \, \frac{k_{\perp }^2}{\tau _{\rm A}}, \,\,\,\, D_{{\rm A}\parallel } \, = \, \left(\frac{v_{\perp }}{V_{\rm A}} \right)^{2} D_{{\rm A}\perp }, \end{equation} \tag{ 33 }$$

where τ_A is the cascade timescale defined below, and v_⊥ is the k_⊥-dependent velocity response of the waves. Note that D_A∥ is independent of k_∥, so it can be pulled out of the derivative in Equation (31). Cranmer & van Ballegooijen (2003) showed that the above form for the diffusion coefficients tends to reproduce the Goldreich & Sridhar (1995) critical balance, and Matthaeus et al. (2009) derived similar functional forms for the coefficients. When specifying the properties of the wavenumber cascade, we apply the scalings for "balanced" turbulence (i.e., zero cross helicity, or Z₊ = Z₋), which is more straightforward to implement but is formally inconsistent with the large-scale transport model of Section 2.

For ideal MHD Alfvénic fluctuations, v²_⊥ is equal to b²_⊥, the latter representing the transverse magnetic variance spectrum divided by 4πρ₀ to convert it to units of velocity squared. Following the usual convention, the power spectrum E_A tracks the magnetic fluctuations, so the reduced spectra are defined formally as

$$\begin{equation} b_{\perp }^{2} \, = \, k_{\perp }^{2} \int dk_{\parallel } \, E_{\rm A}, \,\,\, v_{\perp }^{2} \, = \, \phi b_{\perp }^{2}. \end{equation} \tag{ 34 }$$

The dimensionless factor ϕ describes the departure from ideal MHD energy equipartition. For small values of k_⊥, we assume ϕ ≈ 1. However, as k_⊥ increases into the regime of kinetic Alfvén waves (KAWs), ϕ can become much larger than 1. Hollweg (1999) described how the main difference between v_⊥ and b_⊥ in the KAW regime comes from an enhanced response of the electron velocity distribution to the electric and magnetic fluctuations. For simplicity, we use an approximate analytic expression

$$\begin{equation} \phi \, = \, \frac{\omega ^2}{k_{\parallel }^{2} V_{\rm A}^2} \, \approx \, \frac{1 + k_{\perp }^{2} \rho _{p}^2}{1 + k_{\perp }^{2} \rho _{p}^{2} m_{e} / (\beta m_{p})}, \end{equation} \tag{ 35 }$$

where ρ_p = w_p/Ω_p is the proton thermal gyroradius, with the proton most-probable speed given by w_p = (2k_BT_p/m_p)^1/2 and the proton cyclotron frequency by Ω_p = eB/m_pc. Our term ϕ is equivalent to α² as defined by Howes et al. (2008).

Inspired by Equation (A6), we define the Alfvénic spectral transport timescale as

$$\begin{equation} \tau _{\rm A} \, = \, \frac{1 + \chi _0}{k_{\perp } v_{\perp }}, \end{equation} \tag{ 36 }$$

where we chose to replace the general critical balance parameter χ by its value at the outer-scale parallel wavenumber k_0∥. Thus,

$$\begin{equation} \chi _{0} \, = \, \frac{\omega _0}{k_{\perp } v_{\perp }} \, \approx \, \frac{k_{0 \parallel } V_{\rm A}}{k_{\perp } b_{\perp }}, \end{equation} \tag{ 37 }$$

and χ₀ is the appropriate critical balance parameter to use when solving for the properties of the dominant low-frequency cascade. From Equation (35) we see that KAW outer-scale frequency ω₀ ≈ ϕ^1/2k_0∥V_A, so that a factor of ϕ^1/2 cancels out of both the numerator and denominator to give the final approximate expression above. The wavenumber k_0∥ specifies the spatial scale along the field at which energy is injected in the source term S_A (see below). Because k_0∥ is assumed to be constant (at a given heliocentric distance r), the parameters χ₀ and τ_A are both functions of k_⊥ and not k_∥. The above form for Equation (36) was motivated by the analysis of Zhou & Matthaeus (1990b), Chandran (2008b), and Howes et al. (2011), who described how the cascade and wavenumber anisotropy change when the system transitions from weak (χ₀ ≫ 1) to strong (χ₀ ≪ 1) turbulence.

As mentioned above, our expressions for τ_A, D_A⊥, and D_A∥ assume zero cross helicity (i.e., ${\cal R} = 1$). There is still no agreement about how to generalize these terms when inefficient wave reflection gives rise to nonzero cross helicity. Lithwick et al. (2007) found that the cascade timescales for outward and inward wave modes are different from one another when ${\cal R} \ne 1$, but their parallel spatial scales are the same. However, Beresnyak & Lazarian (2008, 2009) found that k_∥ for the outward mode should be larger than k_∥ for the inward mode, and thus the Goldreich & Sridhar (1995) critical balance must be modified (see Equation (45) below). Chandran (2008b) outlined a method for setting up the advection–diffusion equations in the case of ${\cal R} \ne 1$, but we defer a full implementation of that approach to future work.

Putting aside the issue of imbalanced turbulence, the dominant perpendicular nature of the Alfvénic cascade allows us to define a reduced transport equation that follows the evolution of the spectrum as a function of k_⊥ only. If we ignore the mode coupling term C_AF for now, we can multiply Equation (31) by k²_⊥ and integrate over k_∥ to obtain

$$\begin{equation} \frac{\partial b_{\perp }^2}{\partial t} \, = \, k_{\perp } \frac{\partial }{\partial k_{\perp }} \left[ \frac{1}{\tau _{\rm A}} \left(\alpha _{\perp } k_{\perp } \frac{\partial b_{\perp }^2}{\partial k_{\perp }} - \mu _{\perp } b_{\perp }^{2} \right) \right] + \tilde{S}_{\rm A} - 2 \tilde{\gamma }_{\rm A} b_{\perp }^{2}. \end{equation} \tag{ 38 }$$

This is essentially the same as Equation (11) of Cranmer & van Ballegooijen (2003). The reduced source term $\tilde{S}_{\rm A}$ and dissipation rate $\tilde{\gamma }_{\rm A}$ are defined similarly to the corresponding terms in Equation (31), but they are weighted toward the low-k_∥ regions of wavenumber space that are "filled" by the cascade. In Appendices C.1–C.3, we derive analytic solutions for the time-steady Alfvén wave power spectrum in various limiting cases.

The cascade of fast-mode waves, described by Equation (32), appears to be conceptually simpler than the strongly anisotropic Alfvén wave cascade. The diffusion coefficient is given by D_F = k²/τ_F, where τ_F is related to the IK-like cascade time given by Equation (A2) with p = 1. There is increasing evidence (e.g., Markovskii et al. 2010) that a fast-mode cascade is more rapid in the directions perpendicular to the field than along the field. However, the cascade does appear to proceed outward "radially" in the direction of increasing k. Thus, it makes the most sense to use an isotropic diffusion formalism as in Equation (32), but scale the magnitude of the diffusion timescale with θ. Following the weak turbulence model of Chandran (2005), we adopt

$$\begin{equation} \tau _{\rm F} \, = \, \frac{V_{\rm A}}{k v_{k}^{2} \sin \theta }, \end{equation} \tag{ 39 }$$

which implies that

$$\begin{equation} D_{\rm F} \, = \, \frac{k^{3} v_{k}^{2} \sin \theta }{V_{\rm A}} \, = \, \frac{4\pi k^{6} E_{\rm F} \sin \theta }{V_{\rm A}}. \end{equation} \tag{ 40 }$$

Chandran (2005) showed that the sin θ dependence in the denominator of τ_F is consistent with an isotropic energy flux for the cascade, but it does not guarantee an isotropic wavenumber spectrum E_F(k). More information about how we chose to implement the fast-mode cascade is given in Appendices C.4 and C.5.

In order to fully describe the cascade in the advection–diffusion equations, four dimensionless spectral transport constants (α_⊥, α_∥, μ_⊥, α_F) need to be specified. Matthaeus et al. (2009) summarized the results of many MHD turbulence models and found that α_⊥ often takes on values between 0.2 and 0.5, and α_∥ ≈ 0.43α_⊥ seems to be a useful parameterization (see Equation (13) of Matthaeus et al. 2009). Zhou & Matthaeus (1990b) and Matthaeus et al. (2009) made a case for a classical form of the diffusion operator that implies μ_⊥ = 2α_⊥. Alternately, van Ballegooijen (1986) found that a cascade of random-walk-like displacements of magnetic flux tubes is described well by μ_⊥ = α_⊥. Howes et al. (2008) and Chandran (2008b) used a straightforward advection equation to model an Alfvénic cascade, which sets α_⊥ = 0 and assumes μ_⊥ ≠ 0. For this type of model, Howes et al. (2008) derived μ_⊥  ≈  0.2.

In our models, we are constrained by the values of the cascade constants $\tilde{\alpha }_{\rm A}$ and $\tilde{\alpha }_{\rm F}$ used in the global transport equations of Section 2. We related these constants to the ones defined above by integrating the cascade advection–diffusion terms over wavenumber to find ∂U_m/∂t. By demanding this quantity be equal to the heating rate $\tilde{Q}_{m}$, we obtained

$$\begin{equation} \frac{2 \alpha _{\perp }}{3} + \mu _{\perp } \, = \, \frac{3 \sqrt{6\pi }}{4} \, \tilde{\alpha }_{\rm A}, \end{equation} \tag{ 41 }$$

which assumes that the perpendicular cascade is dominant and that ${\cal R} \approx 1$, and

$$\begin{equation} \alpha _{\rm F} \, = \, \frac{32}{7\pi } \, \tilde{\alpha }_{\rm F}. \end{equation} \tag{ 42 }$$

We keep the ratio s = μ_⊥/α_⊥ as a free parameter and we explore the ramifications of varying it below. Note, however, that if we used s = 2 (as assumed by Zhou & Matthaeus 1990b), then Equation (41) gives α_⊥ ≈ 0.73 and μ_⊥ ≈ 1.47. These are roughly consistent with the constants given by Zhou & Matthaeus (1990b) and Matthaeus et al. (2009). To complete the system of cascade constants, we adopt the Matthaeus et al. (2009) choice for α_∥ = 0.43α_⊥, but we compute this quantity using the Matthaeus et al. (2009) assumption of s = 2.

The source terms, S_A in Equation (31) and S_F in Equation (32), describe the outer-scale injection of fluctuation energy. The global energy balance of the waves is already described by the radial transport model of Section 2. Thus, we specify the magnitudes of S_A and S_F by demanding that the time-steady total energy densities U_A and U_F be maintained at their known values at a given distance r. From a physical standpoint, however, it is unclear whether the passive propagation of waves dominates the source terms, or whether there is significant local "stirring" that converts large-scale dynamical motions (e.g., velocity shears in evolved corotating streams) into new fluctuations.

We adopt specific functional forms for S_A(k_∥, k_⊥) and S_F(k) that are described in detail in Appendix C. Generally, the source terms are nonzero only at the lowest wavenumbers, at which the fluctuations are driven. For the Alfvén waves, we continue to use the assumption from Section 2.3 that the perpendicular driving scale is set by the turbulence correlation length; i.e., k_0⊥ = 1/λ_⊥. For the fast-mode fluctuations, we assume their outer-scale wavenumber magnitude k_0F is also equal to k_0⊥, since the largest-scale transverse stirring motions are likely to be common to both Alfvénic and fast-mode waves. There are several ways that one could imagine defining the parallel outer-scale Alfvén wavenumber k_0∥.

• 1.  
  Monochromatic Alfvén waves that propagate up from the corona retain a constant frequency ω₀ in the Sun's inertial frame. However, because the phase speed varies with distance, the corresponding wavelength undergoes "stretching" commensurate with the dispersion relation

  $$\begin{equation} k_{0 \parallel } \, = \, \frac{\omega _0}{u_{0} + V_{\rm A}}. \end{equation} \tag{ 43 }$$
• 2.  
  The fluctuations propagating up from the Sun may already be fully turbulent (see, e.g., van Ballegooijen et al. 2011). Thus, the outer-scale parallel wavenumber may be coupled continuously to the perpendicular wavenumber via critical balance (Goldreich & Sridhar 1995), with

  $$\begin{equation} k_{0 \parallel } \, \approx \, \frac{k_{0 \perp }}{V_{\rm A}} \sqrt{\frac{U_{\rm A}}{\rho _0}}. \end{equation} \tag{ 44 }$$
• 3.  
  In flux tubes with nonzero cross helicity (i.e., ${\cal R} < 1$), Beresnyak & Lazarian (2008, 2009) found that the inward waves should obey the Goldreich & Sridhar (1995) critical balance, but the outward waves (which are generally what we intend to model) obey a modified version of critical balance, which we approximate as

  $$\begin{equation} k_{0 \parallel } \, \approx \, \frac{k_{0 \perp }}{V_{\rm A}} \sqrt{\frac{U_{\rm A}}{\rho _0}} \frac{1}{\cal R}. \end{equation} \tag{ 45 }$$
• 4.  
  In some cases we assume that the dimensionless ratio k_0∥/k_0⊥ remains fixed at a constant specified value. Many studies of MHD turbulence assume isotropic forcing at the outer scale, which is consistent with the fixed ratio k_0∥/k_0⊥ = 1. The lack of a physical justification for this approximation is offset by its simplicity.

Figure 5 illustrates the ratio k_0∥/k_0⊥ for several of the above methods of setting the parallel outer scale. For example, it shows the result of evaluating Equation (43) for a range of wave periods P = 2π/ω₀ between 1 and 100 minutes. Constant assumed values of k_0∥/k_0⊥ would correspond to horizontal lines in Figure 5.

Zoom In Zoom Out Reset image size

Here we present some example results for the power spectra E_A(k_∥, k_⊥) and E_F(k_∥, k_⊥). These spectra are computed from Equations (31), (32), and (38) in the limiting cases of time independence and no mode coupling (C_AF = 0). The Alfvénic spectrum was first computed in its reduced form using the solutions for b_⊥(k_⊥) given in Appendices C.1 and C.2, and then it was expanded into full wavenumber space by using the results of Appendix C.3. The shape of the fast-mode spectrum was determined from the analytic solutions given in Appendices C.4 and C.5.

To illustrate the wavenumber dependence of the power spectra, we chose a single coronal height z = 10 R_☉ at which β ≈ 0.04. We typically plot the wavenumbers in terms of dimensionless quantities k_∥V_A/Ω_p and k_⊥ρ_p. Dissipative wave–particle interactions tend to become important when these quantities reach order-unity values, and ideal MHD conditions apply when these quantities are small. Typically, the driving scale for Alfvénic turbulence occurs at k_0⊥ρ_p ≈ 10⁻⁶ to 10⁻⁴, with the larger values generally occurring at larger heliocentric distances.

In Figure 6 we show the time-steady k_⊥ dependence for the Alfvénic b_⊥ and v_⊥ fluctuations, both with and without KAW dissipation. To set the cascade properties, we utilized the values of the constants given in Section 3.1, and we also assumed s = μ_⊥/α_⊥ = 2 and k_0∥/k_0⊥ = 0.1. The KAW damping ratio $\tilde{\gamma } / \omega$ appropriate for the assumed value of β, which was used in Equation (C11), is also shown in green (see also Section 5). At the outer scale, the peak value of v_⊥ is 42 km s⁻¹. We caution that this value should not be assumed to be equivalent to the full rms velocity amplitude. In this case (U_A/ρ₀)^1/2 = 196 km s⁻¹, which is almost a factor of five larger than the maximum value of v_⊥ at this height.

Zoom In Zoom Out Reset image size

The damped spectra shown in Figure 6 have several features that resemble those of measured KAWs in the solar wind. Using the conventional form of the reduced energy spectrum (e_A ≈ b²_⊥/k_⊥), we found that the magnetic fluctuation power made a transition from a Kolmogorov-like power law k^−5/3_⊥ to a steeper spectrum with k^−2.5_⊥ at k_⊥ρ_p ≈ 1. The spectrum becomes shallower again around k_⊥ρ_p ≈ 40 because the wavenumber dependence of ϕ flattens out at low values of β. This behavior is reminiscent of that predicted by Voitenko & De Keyser (2011). At larger radial distances where β ≳ 1, the KAW dispersion relation (Equation (35)) gives rise to a more sustained increase in ϕ with increasing k_⊥. This in turn produces spectra that remain steep, with e_A∝k^−2.5_⊥ persisting over several orders of magnitude of k_⊥ in agreement with both measurements (Smith et al. 2006; Sahraoui et al. 2010) and other models (Howes et al. 2008). We note that the predicted undamped KAW power-law decline of k^−7/3_⊥ (see Appendix C.1) was not seen for any sustained range of k_⊥.

Figure 7 shows the result of varying the normalization of the parallel outer-scale wavenumber k_0∥ on the shape of b_⊥(k_⊥). We kept the same value of s = 2 that was used in Figure 6, but we varied the constant ratio k_0∥/k_0⊥ over five orders of magnitude. For the lowest values of k_0∥ the outer-scale critical balance ratio χ₀ always remains much smaller than unity. This means that the stirring or forcing takes place well within the "filled" region of wavenumber space, and thus strong turbulence occurs. In this case, b_⊥∝k^−1/3_⊥ and thus e_A∝k^−5/3_⊥. The opposite extreme case of large k_0∥ corresponds to χ₀ ≫ 1 and weak turbulence with less anisotropic driving. In that limit, the inertial range spectra are given by b_⊥∝k^−1/2_⊥ and e_A∝k⁻²_⊥. Our model shows the gradual transition between these two extreme cases.

Zoom In Zoom Out Reset image size

In Figure 8 we compare the Alfvén and fast-mode spectra with one another. As above, we used the background conditions at a coronal height of z = 10 R_☉ and we assumed k_0∥/k_0⊥ = 0.1. We illustrate the most extreme case of a lack of high-frequency Alfvénic power by showing the contours of E_A(k_∥, k_⊥) for the case s → ∞. In this limit, Equation (C16) describes an exponential decrease of power with increasing χ. Other comparable examples of this kind of spectrum can be found in Figure 4(b) of Cranmer & van Ballegooijen (2003) and Figures 1 and 2 of Jiang et al. (2009). We computed the Alfvénic and fast-mode spectra with the kinetic sources of damping that were described in Section 5. Note that E_F experiences the strongest damping at intermediate values of θ. For θ ≲ 10° or θ ≳ 85°, the transit-time damping described by Equation (58) is relatively weak.

Zoom In Zoom Out Reset image size

There are several ways that the ideal linear MHD wave modes can become coupled to one another in the corona and solar wind.

• 1.  
  Inhomogeneities in the background plasma can blur the definitions of the individual modes. For example, linear reflection due to radial variations in V_A (Ferraro & Plumpton 1958; Heinemann & Olbert 1980) may produce not only incoming Alfvén waves (i.e., $0 < {\cal R} < 1$), but also fast and slow magnetosonic waves (e.g., Stein 1971; McDougall & Hood 2007). In addition, large-scale bends in the background magnetic field B (Frisch 1964; Wentzel 1974), density variations between flux tubes (Valley 1974; Markovskii 2001; Mecheri & Marsch 2008), or velocity shears (Poedts et al. 1998; Gogoberidze et al. 2007) can drive instabilities that partially convert Alfvén waves into other modes.
• 2.  
  Even in a homogeneous medium, the MHD waves begin to lose their ideal linear character when their amplitudes become large. Nonlinear Alfvén waves naturally drive second-order fluctuations in v_∥ and δρ that mimic the properties of both slow and fast magnetosonic waves (Hollweg 1971; Spangler 1989; Vasquez & Hollweg 1999). Large-amplitude waves also excite a range of wave–wave interactions that can often be characterized either as two modes giving birth to a third, or one mode splitting into several others (e.g., Chin & Wentzel 1972; Goldstein 1978; Del Zanna et al. 2001; Sharma & Kumar 2010). Models of weak turbulence, in which the wave–wave interactions describe the cascade process (Chandran 2005, 2008a; Luo & Melrose 2006; Yoon & Fang 2008) also create this kind of coupling.
• 3.  
  Although not strictly a multi-mode coupling, when k_⊥ρ_p ≳ 1 the Alfvén mode begins to exhibit oscillations in density, parallel electron velocity, and the parallel electric and magnetic fields (Hasegawa & Chen 1976; Hollweg 1999). Observationally, it has proved difficult to separate such dispersive KAW density fluctuations from those arising from independent sources of fast or slow MHD waves (e.g., Harmon & Coles 2005; Chandran et al. 2009).

In this paper, we take account of one particular nonlinear effect from the second entry in the above list. Specifically, Chandran (2005) suggested that weak turbulence couplings between Alfvén and fast-mode fluctuations may provide enough power at high k_∥ to induce substantial ion cyclotron heating. Suzuki et al. (2007) argued that this effect may be relatively unimportant because the fast-mode cascade timescale τ_F is long in comparison to the Alfvén cascade timescale τ_A. This may be the case in the low-frequency regime of wavenumber space where χ ≪ 1, but at the cyclotron resonant frequencies of interest (k_∥ ∼ Ω_p/V_A) the Alfvénic cascade is quenched because χ ≫ 1. The fast-mode cascade may in fact even be faster than any intrinsic Alfvénic spectral transfer in this region of wavenumber space. Therefore, we proceed using the Chandran (2005) results for Alfvén/fast-mode coupling.

We express the coupling term in Equations (31) and (32) as

$$\begin{equation} C_{\rm AF} \, = \, \frac{E_{\rm F} - E_{\rm A}}{\tau _{\rm AF}} \end{equation} \tag{ 46 }$$

such that, in the absence of other processes, the power spectra at a given wavenumber k are driven toward a common value over a coupling timescale τ_AF(k). The weak turbulence model of Chandran (2005) gave an approximate value for this timescale of

$$\begin{equation} \tau _{\rm AF} \, \approx \, \frac{15}{23 \pi ^2} \tau _{\rm F} \, \sin ^{2} \theta \end{equation} \tag{ 47 }$$

which holds in the limiting cases of E_F > E_A and nearly parallel propagation (θ ≪ 1). In the opposite case of E_A ≫ E_F, it is likely that τ_AF would no longer depend linearly on τ_F, and may scale instead with τ_A. However, the region of wavenumber space with which we are most concerned is the high-k_∥, low-θ ion cyclotron regime. At those wavenumbers, we know that in the absence of coupling the condition E_F ≫ E_A is likely to be satisfied, and the coupling will be a transfer of energy from the dominant fast-mode spectrum to the much less intense Alfvén mode.

The wave–wave conditions of frequency and wavenumber matching (e.g., Sagdeev & Galeev 1969) confirm that the most rapid coupling should occur when the dispersive properties of the Alfvén and fast-mode waves are the most similar to one another; i.e., at θ → 0. Note that Equation (39) gave τ_F∝1/sin θ, so the combined dependence for the coupling time is τ_AF∝sin θ. In practice, however, we found that using this ideal expression for τ_AF could lead to an unphysical singularity at θ = 0. We removed this singularity by replacing θ in Equation (47) by θ + δθ. We set δθ to a constant value of 0.01 to avoid having an infinitely fast coupling rate at parallel propagation.³ To retain the most generality, we chose to reparameterize the coupling timescale as

$$\begin{equation} \tau _{\rm AF} \, = \, \frac{1}{\Phi } \, \tau _{\rm F} \sin ^{2} (\theta + \delta \theta), \end{equation} \tag{ 48 }$$

where we find it useful to vary the constant coupling strength Φ up or down from the value of 23π²/15 ≈ 15.1 derived by Chandran (2005). The case Φ = 0 corresponds to ignoring the coupling altogether.

Note that the above form for the coupling timescale implies that τ_AF∝k_⊥/k^3/2, so that the coupling is rapid at wavenumbers corresponding to ion cyclotron resonance (large k_∥, small k_⊥). The coupling is much slower at KAW wavenumbers favored by the pure Alfvénic cascade (small k_∥, large k_⊥). Thus, the bulk of the Alfvénic spectrum at χ ≪ 1 is likely to be more or less unaffected by the coupling. This seems to be consistent with our assumption that the integrated energy densities U_A and U_F also remain uncoupled from one another. We realize that this may be a severe underestimate of the degree of energy transfer between Alfvén and magnetosonic modes in the corona and solar wind. However, one main purpose of this paper is to investigate how much can be accomplished with only this small degree of coupling in the high-k_∥ tails of the power spectra.

The exact solutions to Equations (31) and (32) with coupling (C_AF ≠ 0) must be found numerically. Here we present an approximate solution that is both (1) likely to reflect the proper behavior of more rigorous numerical solutions in many limiting regimes of parameter space, and (2) efficient to implement on a large grid of model spectra spanning a wide range of heliocentric distances. We begin by approaching the problem iteratively; i.e., we solve Equation (31) for E_A under the assumption that E_F is known, and we then solve Equation (32) for E_F under the assumption that E_A is known. The analytic solutions derived below suggest a natural way to terminate this iteration after only one round.

When solving the advection–diffusion equation for Alfvénic fluctuations, let us temporarily ignore the outer-scale source term S_A and the dissipation-range damping term that depends on γ_A. Since we are most concerned with the generation and transport of wave power in the high-k_∥ regions that undergo ion cyclotron resonance, we consider the weak turbulence regime of χ ≫ 1, in which the transport of energy is mainly from low to high k_⊥ and there is negligible parallel spreading (see also Oughton et al. 2006). Thus, we solve the advection–diffusion equation for discrete, non-interacting "strips" of wavenumber space each having constant k_∥. The nonlinear coupling supplies wave energy locally, and the Alfvénic cascade takes it from low to high k_⊥. If we simplify further by assuming pure advection (i.e., α_⊥ = 0), the time-steady version of Equation (31) becomes

$$\begin{equation} \frac{\mu _{\perp }}{k_{\perp }} \frac{\partial }{\partial k_{\perp }} \left(\frac{k_{\perp }^{2} E_{\rm A}}{\tau _{\rm A}} \right) \, = \, \frac{E_{\rm F} - E_{\rm A}}{\tau _{\rm AF}}, \end{equation} \tag{ 49 }$$

where we use Equation (A6) to give the timescale τ_A ≈ χ/(k_⊥v_⊥) in the weak turbulence regime, and we use Equation (48) for τ_AF.

The advection–coupling equation above can be rewritten as a first-order ordinary differential equation,

$$\begin{equation} \frac{\partial E_{\rm A}}{\partial k_{\perp }} + \left(\frac{10}{3k_{\perp }} + \frac{f_0}{k_{\perp }^{10/3}} \right) E_{\rm A} \, = \, \frac{f_{0} E_{\rm F}}{k_{\perp }^{10/3}}, \end{equation} \tag{ 50 }$$

where

$$\begin{equation} f_{0} \, = \, \frac{\Phi }{\mu _{\perp }} \left(\frac{v_{0 {\rm F}}}{v_{0 \perp }} \right)^{2} \frac{k_{0 {\rm F}}^{1/2} k_{\parallel }^{5/2}}{k_{0 \perp }^{2/3}}. \end{equation} \tag{ 51 }$$

To obtain Equation (50), we made several power-law assumptions for the timescales τ_A and τ_AF, which depend on the velocity spectra v_⊥ (for Alfvén waves) and v_k (for fast-mode waves), respectively, with

$$\begin{equation} v_{\perp } \, = \, v_{0 \perp } \left(\frac{k_{\perp }}{k_{0\perp }} \right)^{-1/3}, \,\,\,\, v_{k} \, = \, v_{0 {\rm F}} \left(\frac{k}{k_{0 {\rm F}}} \right)^{-1/4}. \end{equation} \tag{ 52 }$$

We also assumed that we are solving for E_A mainly in the small-θ region of wavenumber space in which k ≈ k_∥.

With the above assumptions taken into account, Equation (50) can be solved by means of an integrating factor. We first define the dimensionless independent variable

$$\begin{equation} y \, = \, \frac{3f_0}{7 k_{\perp }^{7/3}} \, = \, \left(\frac{k_c}{k_{\perp }} \right)^{7/3} \end{equation} \tag{ 53 }$$

which is a measure of the relative strength of the nonlinear coupling. When y ≫ 1 (or k_⊥ ≪ k_c) the coupling is strong and we should expect E_A ≈ E_F. When y ≪ 1 (or k_⊥ ≫ k_c) the coupling is weak in comparison to the cascade and we expect E_A ≪ E_F. Note also that y depends much more sensitively on θ than on the magnitude k. Working through the integrating factor method and choosing an integration constant of zero (to avoid the solution diverging to infinity when y ≫ 1), we obtain

$$\begin{equation} E_{\rm A} \, = \, \frac{7y}{3} \left[ 1 - e^{y} y^{3/7} \, \Gamma \left(\frac{4}{7}, \, y \right) \right] E_{\rm F}, \end{equation} \tag{ 54 }$$

where Γ(a, y) is the incomplete gamma function. This function behaves as expected in the limits of strong and weak coupling as discussed above.

Next, we solve the coupled fast-mode advection–diffusion equation for E_F under the assumption that E_A is known. Making use of many of the same simplifications that were used to solve the E_A equation, we include only the cascade and coupling terms, with

$$\begin{equation} \frac{\alpha _{\rm F}}{k^2} \frac{\partial }{\partial k} \left(k^{2} D_{\rm F} \frac{\partial E_{\rm F}}{\partial k} \right) \, = \, C_{\rm AF} \, = \, \left[ \frac{1}{\tau _{\rm AF}} \left(1 - \frac{E_{\rm A}}{E_{\rm F}} \right) \right] E_{\rm F}. \end{equation} \tag{ 55 }$$

Noticing that the quantity in square brackets above is an effective damping rate γ_eff, we use Equation (54) to write the ratio E_A/E_F as a known function of k_∥ and k_⊥. After substituting in the wavenumber dependence for τ_AF, we found that γ_eff∝(k/k_0F)^1/2. The analytic solution of E_F(k) for this special case is given in Equation (C32), and the constant c_γ in that expression is specified here to be

$$\begin{equation} c_{\gamma } \, = \, \frac{8 \Phi }{49 \, \alpha _{\rm F} \sin ^{2}\theta } \left(1 - \frac{E_{\rm A}}{E_{\rm F}} \right). \end{equation} \tag{ 56 }$$

The solution of Equation (C32) is applied only for k ⩾ k_0F, and the uncoupled/undamped fast-mode power spectrum E_0F is used for k < k_0F.

Since our solution for the ratio E_A/E_F depends only on wavenumber and not on any prior solutions of E_A or E_F, we found that there is no need for further iteration. We solve first for E_F as described above, using Equation (54) for the ratio E_A/E_F, and then we use this ratio to solve for E_A. Note that the complete solution for E_F must take account of both coupling and transit-time damping (i.e., the damping rate given by Equation (58)). In practice, we apply both types of damping separately to the uncoupled and undamped fast-mode power spectrum E_0F and we use the solution that gives rise to stronger local damping at any given wavenumber. At high enough values of k_∥, the complete solution for E_A must take into account the effects of ion cyclotron damping. We use the approximate prescription given by Equation (C20) to implement this damping.

If the original uncoupled spectra obey E_0A ⩽ E_0F, then the coupled spectra follow

$$\begin{displaymath} E_{0 {\rm A}} \le E_{\rm A} \le E_{\rm F} \le E_{0 {\rm F}} \end{displaymath}$$

at wavenumbers in the high-k_∥ regime where the coupling is applied. Usually, the relative increase in E_A from its uncoupled solution is greater than the relative decrease in E_F from its uncoupled solution. In all cases, however, we found that the variations in the spectra introduced by the coupling do not significantly affect the total wavenumber-integrated power in either E_A or E_F.

Figure 9 illustrates the effects of including coupling on E_A. As in earlier plots of spectrum results, we used the representative height z = 10 R_☉ and we assumed k_0∥ = k_0⊥/10. In order to show that the coupling can be efficient even when the uncoupled Alfvén wave power E_0A is negligibly small, we assumed the extreme limiting case of s → ∞. In Figure 9(a) we show the k_∥ dependence of the spectra along a slice taken at a constant value of k_⊥ = k_0⊥. We varied the parameter Φ between 10⁻⁶ and 10⁺³. Even if the coupling is several orders of magnitude weaker than estimated by Chandran (2005), it is still likely to be efficient at generating some Alfvénic wave power at k_∥ ≈ Ω_p/V_A. However, if the coupling constant Φ is significantly smaller than ∼10⁻³, the ion cyclotron damping at k_∥ ≈ Ω_p/V_A is likely to overwhelm the "local supply" of wave energy from the coupling and give rise to a low level of resonant wave power.

Zoom In Zoom Out Reset image size

Figure 9(b) shows how the power at a given wavenumber (k_⊥ = k_0⊥ and k_∥V_A/Ω_p = 10⁻³) varies as a function of Φ. The fast-mode power decreases monotonically as Φ is increased, which confirms our treatment of the coupling in Equation (55) as an effective damping. The Alfvénic power generally increases (from its uncoupled value far below the lower edge of the plot) with increasing Φ, but there is some nonmonotonicity around Φ ≈ 10⁻². This gives rise to a slightly counterintuitive result that there may be more E_A power at high k_∥ (and thus more proton and ion heating) at some values of Φ than in the Φ → ∞ limit.

An example of the full wavenumber dependence of the coupled E_A(k_∥, k_⊥) spectrum is shown in Figure 10(a) for a radial distance of r = 10 R_☉. This model has the same parameters as the one shown in Figure 8, except that we set Φ = 10. Despite the appearance of substantial wave power at large values of k_∥, most of the power is still contained within the critical balance locus of χ ≲ 1. This is illustrated in another way by Figure 10(b), in which we show the radial dependence of the spectrum-averaged angle Θ_Bk between the background field direction and the wavenumber vector k. We used a definition for the spectrum-averaged wavevector anisotropy that is similar to that of Gary et al. (2010),

$$\begin{equation} \tan ^{2} \langle \Theta _{{\rm B}k} \rangle \, = \, \frac{\int d^{3}{\bf k} \,\, E_{\rm A}({\bf k}) \, k_{\perp }^2}{\int d^{3}{\bf k} \,\, E_{\rm A}({\bf k}) \, k_{\parallel }^2}. \end{equation} \tag{ 57 }$$

Note that the model result at r = 1 AU (895) is reasonably close to the value of ∼88° measured by Sahraoui et al. (2010) from the four Cluster satellites at 1 AU. It is evident that a strongly perpendicular ("quasi-two-dimensional") sense of wavenumber anisotropy is not incompatible with the existence of high-frequency ion cyclotron resonant wave power.

Zoom In Zoom Out Reset image size

The cascade of energy from large to small eddies was first described in the context of isotropic hydrodynamic turbulence (von Kármán & Howarth 1938; Kolmogorov 1941; Obukhov 1941; Batchelor 1953). The spectral transport timescale for energy to be transferred down to the next order of magnitude of eddy size is estimated generally as τ_s ≈ (kv_k)⁻¹, where k is the magnitude of the local wavevector k and v_k is the local eddy velocity at this value of k. For isotropic fluctuations that depend only on k and not its direction, we can define the reduced one-dimensional spectrum e_m(k) = v²_k/k. Thus, since

$$\begin{equation} \frac{U_m}{\rho _0} \, = \, \int dk \,\, e_{m} (k), \end{equation} \tag{ A1 }$$

we relate the eddy velocity to the full three-dimensional spectrum via v²_k = 4πk³E_m. The cascade rate is estimated as ε ∼ v²_k/τ_s. Assuming that ε is constant in the inertial range leads to the time-steady Kolmogorov–Obukhov spectrum e_m∝k^−5/3, or E_m∝k^−11/3.

When the background magnetic field becomes strong, other physical processes become important. Iroshnikov (1963) and Kraichnan (1965) (hereafter IK) realized that the "eddy" description of hydrodynamic turbulence could be generalized by referring to colliding MHD wave packets, and that the Alfvén speed V_A introduces a new absolute scale into the problem. If one continues to treat the cascade isotropically in k-space, a more generalized spectral transport time can be defined as

$$\begin{equation} \tau _{s} \, = \, \frac{1}{k v_{k}} \left(\frac{V_{\rm A}}{v_{k}} \right)^{p}, \end{equation} \tag{ A2 }$$

where p = 0 gives the Kolmogorov–Obukhov limit and p = 1 is the result of the IK analysis. Using the same assumption above that ε is constant, we obtain a more general one-dimensional power spectrum e_m∝k^{−(p + 5)/(p + 3)}. For the IK value of p = 1, the spectrum is e_m∝k^−3/2 (see also Boldyrev 2005).

It has been known for several decades that a cascade of Alfvén-wave-like fluctuations does not lead to an isotropic distribution of power in wavenumber space (Strauss 1976; Montgomery & Turner 1981; Shebalin et al. 1983; Higdon 1984). The dominant energy cascade takes place mainly in the two-dimensional plane perpendicular to the background field. For the Alfvénic fluctuations, we can define the local eddy velocity as v_⊥ being mainly a function of k_⊥. The one-dimensional spectrum in this case is given by e_A = v²_⊥/k_⊥ and the integration over wavenumber space is best done in cylindrical coordinates with

$$\begin{equation} \frac{U_{\rm A}}{2\pi \rho _0} \, = \, \int dk_{\parallel } \int dk_{\perp } \, k_{\perp } E_{\rm A} \, = \, \int dk_{\perp } \, e_{\rm A}. \end{equation} \tag{ A3 }$$

Taking into account the spectral anisotropy (k_∥ ≠ k_⊥) we can also write an even more general perpendicular transport time for the Alfvén waves as

$$\begin{equation} \tau _{\rm A} \, = \, \frac{1}{k_{\perp } v_{\perp }} \left(\frac{V_{\rm A}}{v_{\perp }} \right)^{p} \left(\frac{k_{\parallel }}{k_{\perp }} \right)^{q}. \end{equation} \tag{ A4 }$$

A perpendicular generalization of the IK model is given by p = 1 and q = 0, which gives e_A∝k^−3/2_⊥ (see also Nakayama 1999, 2001; Boldyrev 2006; Podesta 2011). Weak three-wave couplings have been shown to give rise to the case p = q = 1, which yields e_A∝k⁻²_⊥ (e.g., Galtier et al. 2000; Bhattacharjee & Ng 2001; Boldyrev & Perez 2009). However, in that case nonlinear effects grow in magnitude as k_⊥ gets larger, so it is generally believed that a weakly turbulent inertial range must eventually become strongly turbulent (see also Goldreich & Sridhar 1997).

Goldreich & Sridhar (1995) described strong Alfvénic turbulence with a spectral transfer time given by p = q = 0, and thus e_A∝k^−5/3_⊥ reminiscent of the original Kolmogorov–Obukhov model. In this case of strong mixing between the turbulent motions (perpendicular to the field) and the flow of Alfvén wave packets (parallel to the field) there is a so-called critical balance that couples k_⊥ and k_∥ to one another. We define a critical balance parameter

$$\begin{equation} \chi \, \approx \, \frac{k_{\parallel } V_{\rm A}}{k_{\perp } v_{\perp }} \end{equation} \tag{ A5 }$$

such that the Goldreich & Sridhar (1995) strong cascade is consistent with the condition χ ≈ 1. Combining this with the velocity scaling v_⊥∝k^−1/3_⊥ yields the wavenumber anisotropy scaling k_∥∝k^2/3_⊥. Note that assuming p = q in Equation (A4) is equivalent to τ_A being given by χ^p/(k_⊥v_⊥); see also Galtier et al. (2005). An alternate way of describing the cascade was given by Zhou & Matthaeus (1990b), who defined a triple correlation timescale equivalent to

$$\begin{equation} \tau _{\rm A} \, \approx \, \frac{1 + \chi }{k_{\perp } v_{\perp }}. \end{equation} \tag{ A6 }$$

This expression naturally bridges the strong (χ ≲ 1) and weak (χ ≫ 1) turbulence scaling limits, and we use a similar relation in Section 3.1.

The cascade of compressible fast-mode waves has received less attention than that of the incompressible Alfvén waves. Because fast-mode waves propagate at a roughly constant phase speed no matter the direction angle θ, the fast-mode cascade has been suspected to resemble isotropic hydrodynamic turbulence. In fact, numerical simulations do tend to find that fast-mode waves produce a more isotropic spectrum than do Alfvén waves (Cho & Lazarian 2003; Svidzinski et al. 2009).⁴ The rate of the cascade is generally assumed to follow the weak IK-type scaling of Equation (A2) with p = 1 (see, e.g., Chandran 2005; Suzuki et al. 2007). Thus, because in most cases we expect v_k ≪ V_A, the fast-mode cascade timescale τ_F is likely to be significantly longer than the Alfvénic timescale τ_A.

It is important to emphasize that there is still no agreement concerning the most realistic and universal way to describe MHD turbulence. There remains controversy about the applicability of the various power-law exponents (especially 5/3 versus 3/2) for the Alfvénic inertial range (Beresnyak 2011; Forman et al. 2011; Mason et al. 2011; Podesta 2011). Simulations have not been able to accurately pin down the amount of slow "leakage" of power to the high-k_∥ region of the spectrum where χ ≫ 1. Furthermore, the observed steepening of the spectrum at high values of k_⊥ is still not well understood (Leamon et al. 1998; Stawicki et al. 2001; Howes et al. 2008). In many models, the precise scalings depend on the degree of cross helicity of the fluctuations (i.e., on the imbalance between Z₊ and Z₋) and on whether the turbulence is driven or decaying (e.g., Lithwick et al. 2007; Chandran 2008b; Chen et al. 2011). In this paper, we attempt to identify the most controversial aspects of the models and discuss how they can be modified once such issues are resolved.

Equation (38) is a reduced one-dimensional version of the full advection–diffusion equation for Alfvénic fluctuations. The Appendix of Cranmer & van Ballegooijen (2003) presented one method of solving this equation in the low-wavenumber, strong turbulence (χ₀ ≪ 1) limit. Here we derive a more general case for arbitrary χ₀. Ignoring both wave damping and mode coupling, and assuming a steady state (i.e., ∂b²_⊥/∂t = 0), Equation (38) can be simplified to

$$\begin{equation} \frac{\partial \varepsilon }{\partial x} \, = \, \tilde{S}_{\rm A}, \end{equation} \tag{ C1 }$$

where here we define x = ln k_⊥ and we write the cascade rate as

$$\begin{eqnarray} \varepsilon & = & \frac{1}{\tau _{\rm A}} \left(\mu _{\perp } b_{\perp }^{2} - \alpha _{\perp } \frac{\partial b_{\perp }^{2}}{\partial x} \right) \end{eqnarray} \tag{ C2 }$$

$$\begin{eqnarray} \hspace{7pt}& = &\! - \frac{\alpha _{\perp } k_{\perp }^{1+s}}{\tau _{\rm A}} \frac{\partial }{\partial k_{\perp }}(b_{\perp }^{2} k_{\perp }^{-s}), \end{eqnarray} \tag{ C3 }$$

where s = μ_⊥/α_⊥. Note that the ratio s was called β/γ by Cranmer & van Ballegooijen (2003) and Landi & Cranmer (2009). The second form for ε given in Equation (C3) helps to show that the power-law spectrum for b_⊥ in the inertial range (i.e., where $\tilde{S}_{\rm A} = 0$ and ε is constant) should be independent of the value of s. In the limiting cases of strong (χ₀ ≪ 1) and weak (χ₀ ≫ 1) turbulence, we use Equation (36) to find that b_⊥ is proportional to k^−1/3_⊥ and k^−1/2_⊥, respectively.

In regions of wavenumber space where the source term is nonzero, ε is not constant and the simple inertial-range scalings do not apply. If we assume that most of the fluctuation power is injected near x ≈ x₀ = ln k_0⊥, then it makes sense to use a compact Gaussian shape for the source term,

$$\begin{equation} \tilde{S}_{\rm A}(x) \, = \, \frac{\varepsilon _0}{\pi ^{1/2} \sigma _0} \exp \left[ - \left(\frac{x - x_0}{\sigma _0} \right)^{2} \right], \end{equation} \tag{ C4 }$$

where the dimensionless width of the Gaussian is specified by σ₀ = 1 in our models. The constant ε₀ is varied arbitrarily to produce the desired total fluctuation energy density U_A. Cranmer & van Ballegooijen (2003) showed how the above form for the source function integrates to a cascade rate

$$\begin{equation} \varepsilon (x) \, = \, \frac{\varepsilon _0}{2} \left[ 1 + \mbox{erf} \left(\frac{x - x_0}{\sigma _0} \right) \right]. \end{equation} \tag{ C5 }$$

Finally, we define an auxiliary parameter q = b²_⊥k_⊥^−s and integrate Equation (C3) to obtain

$$\begin{equation} q^{3/2} \, = \, \frac{3}{2 \alpha _{\perp }} \int _{k_{\perp }}^{\infty } \frac{d\kappa \, (1 + \chi _{0}) \, \varepsilon (\kappa)}{\phi ^{1/2} \kappa ^{2 + (3s/2)}}. \end{equation} \tag{ C6 }$$

In practice, we integrate this equation numerically and use the definition of q to obtain b_⊥(k_⊥). In the energy containing range (x ≪ x₀), we see that ε → 0, and thus q is constant and b²_⊥∝k_⊥^s. The low-k_⊥ region of wavenumber space stands in contrast to the inertial range because here the shape of the fluctuation spectrum does depend on the value of s.

In the MHD strong-turbulence inertial range (where ε ≈ ε₀, ϕ ≈ 1, and χ₀ ≪ 1), we obtain the standard solution b_⊥∝v_⊥∝k^−1/3_⊥. However, when k_⊥ρ_p increases past unity into the KAW dispersive range, we see from Equation (35) that for β ≫ 1 there exists a sizable "dispersion range" in which ϕ∝k²_⊥ and we obtain

$$\begin{equation} b_{\perp } \propto k_{\perp }^{-2/3}, \,\,\,\, v_{\perp } \propto k_{\perp }^{+1/3}. \end{equation} \tag{ C7 }$$

Converting these into the more commonly used one-dimensional spectra (see Appendix A), the MHD inertial range has e_A∝k^−5/3_⊥, and the KAW inertial range has e_A∝k^−7/3_⊥ for the magnetic fluctuations and e_A∝k^−1/3_⊥ for the (electron) velocity fluctuations. These scalings have been described by, e.g., Biskamp et al. (1996), Cranmer & van Ballegooijen (2003), and Howes et al. (2008). Note, however, that strong damping also begins to occur in the KAW regime, so the above power laws may not be evident in the final modeled spectra.

We include the effects of high-k_⊥ dissipation in the perpendicular cascade by assuming the damping acts only at wavenumbers significantly above those where the source term is dominant. Thus, we assume that $\tilde{S}_{\rm A} = 0$ and we continue to ignore mode coupling. In this limiting case, the time-steady advection–diffusion equation becomes

$$\begin{equation} k_{\perp } \frac{\partial \varepsilon }{\partial k_{\perp }} \, = \, - 2 \tilde{\gamma }_{\rm A} b_{\perp }^{2}. \end{equation} \tag{ C8 }$$

We note that Howes et al. (2008) found it is important to include KAW damping when solving for the steady-state wave power (and thus the proton/electron energy partitioning) at large values of k_⊥.

In Appendix C.1 we showed that the properties of the inertial range spectra should be independent of the value of s. In order to produce a closed-form solution, here we follow Howes et al. (2008) and assume that s → ∞ (i.e., the cascade proceeds purely by wavenumber advection). Thus, by assuming that α_⊥ = 0, we can use Equation (C2) to write the b²_⊥ term on the right-hand side of Equation (C8) as

$$\begin{equation} b_{\perp }^{2} \, \approx \, \frac{\varepsilon }{\mu _{\perp } k_{\perp } v_{\perp }} \end{equation} \tag{ C9 }$$

if we also assume χ₀ ≪ 1 in the high-k_⊥ limit. To retain the most generality in cases when s is not infinitely large, we can use Equation (41) to replace μ_⊥ by μ_⊥* = μ_⊥ + 2α_⊥/3. This may not be completely accurate for the cases of weakest advection (i.e., s ≈ 0), but it is an improvement on ignoring the influence of the α_⊥ diffusion term altogether.

To simplify the modified transport equation, we take a further cue from Howes et al. (2008) and use the critical balance condition ω ≈ k_⊥v_⊥ to rewrite the equation as

$$\begin{equation} \frac{k_{\perp }}{\varepsilon } \frac{\partial \varepsilon }{\partial k_{\perp }} \, = \, - \frac{2}{\mu _{\perp \ast }} \left(\frac{\tilde{\gamma }}{\omega } \right)_{\rm A}, \end{equation} \tag{ C10 }$$

where the ratio $(\tilde{\gamma }/ \omega)_{\rm A}$ is the output of the Vlasov–Maxwell dispersion analysis discussed in Section 5. For the KAW domain, this ratio is largely independent of k_∥ and thus it can be treated as a function of k_⊥ only. When we solve this equation numerically, we start the integration at a low enough value of k_⊥ that the damping is negligibly small (i.e., at which ε = ε₀). Thus, we integrate upward in k_⊥ with

$$\begin{equation} \varepsilon \, = \, \varepsilon _{0} \, \exp \left[ - \frac{2}{\mu _{\perp \ast }} \int \frac{dk_{\perp }}{k_{\perp }} \left(\frac{\tilde{\gamma }}{\omega } \right)_{\rm A} \right]. \end{equation} \tag{ C11 }$$

Finally, the damped solution for ε(k_⊥) is used in Equation (C6) to obtain the damped power spectrum.

The previous sections described the cascade as a function of k_⊥ and ignored the behavior of the full power spectrum E_A(k_∥, k_⊥). To first order, the strong predicted anisotropy of MHD turbulence justifies this approach, but we are also concerned with the possible leakage of power to high values of k_∥ and thus to high frequencies. Here we mainly follow the analysis of Cranmer & van Ballegooijen (2003), but we also include the possible effects of weak turbulence when χ₀ ≫ 1. Goldreich & Sridhar (1995) wrote the full power spectrum as a separable function of two variables, k_⊥ and χ, with

$$\begin{equation} E_{\rm A} (k_{\parallel }, k_{\perp }) \, = \, \frac{V_{\rm A} b_{\perp }(k_{\perp })}{k_{\perp }^3} g(\chi) \end{equation} \tag{ C12 }$$

and χ = k_∥V_A/(k_⊥b_⊥). This definition allows the dimensionless function g(χ) to be normalized to unity,

$$\begin{equation} \int _{-\infty }^{+\infty } d\chi \, g(\chi) \, = \, 1, \end{equation} \tag{ C13 }$$

but in practice we usually calculate the normalization for g(χ) from the condition that the total power over all wavenumber space integrates properly to U_A.

It has been known for some time that the dominant contribution to the integral in Equation (C13) should come from the region where |χ| ≲ 1. For values of k_∥ at which |χ| ≲ 1, the Goldreich & Sridhar (1995) solution for b_⊥∝k^−1/3_⊥ gives a dominant perpendicular dependence for the inertial range of E_A∝k^−10/3_⊥. We also expect g(χ) to grow negligibly small for |χ| ≫ 1, and thus for these large-k_∥ regions of wavenumber space there should be very little Alfvénic wave power. Cho et al. (2002) found that numerical simulations of anisotropic MHD turbulence were consistent with g(χ) being fitted reasonably well with either a simple exponential function (g ∼ e^−χ) or a Castaing function (a convolution of multiple exponentials). Cranmer & van Ballegooijen (2003) derived an analytic solution to a simplified version of Equation (31), with

$$\begin{equation} g (\chi) \, = \, \frac{2\Gamma (n)}{3\Gamma (n - 0.5) \sqrt{\pi }} \left(1 + \frac{4 \alpha _{\perp } \chi ^{2}}{9 \alpha _{\parallel }} \right)^{-n} \end{equation} \tag{ C14 }$$

and

$$\begin{equation} n \, = \, 1 + \frac{3 s}{4}. \end{equation} \tag{ C15 }$$

These expressions are appropriate for the MHD inertial range where v_⊥∝k^−1/3_⊥, but we use them for the entire range of modeled wavenumbers. The above form for g(χ) resembles a generalized Lorentzian, or kappa distribution (e.g., Vasyliunas 1968; Pierrard & Lazar 2010) that is Gaussian for small arguments and evolves to a power-law tail for large arguments. We can simplify the argument of the power-law term by using the values of the cascade parameters discussed in Section 3.1, with α_⊥/α_∥ ≈ 18.6/(3s + 2).

An example choice for the dimensionless advection–diffusion ratio (s = μ_⊥/α_⊥ = 2) gives rise to a power-law exponent n = 5/2, and the large-k_∥ behavior of the Alfvénic power spectrum is E_A∝k⁻⁵_∥. Smaller values of s give shallower power-law slopes. In fact, Cranmer & van Ballegooijen (2003) and Landi & Cranmer (2009) found that if s could be maintained at small values of order 0.1–0.3, there would be sufficient high-k_∥ power to heat protons and heavy ions in the corona via ion cyclotron resonance. In the opposite limit of pure advection (i.e., s → ∞ or α_⊥ → 0, with α_∥ and μ_⊥ remaining finite) Equation (C14) becomes a Gaussian,

$$\begin{equation} g(\chi) \, \propto \, \exp \left(- \frac{\mu _{\perp }}{3 \alpha _{\parallel }} \, \chi ^{2} \right). \end{equation} \tag{ C16 }$$

Chandran (2008b) also obtained a similar Gaussian solution for the parallel spectrum under the assumption of pure advection. Using the values of the cascade parameters discussed in Section 3.1, we can set μ_⊥ ≈ 1.95 in the limit of s → ∞. Thus, the ratio μ_⊥/α_∥ ≈ 6.2 and we can write $g \sim e^{-2 \chi ^{2}}$ in the pure-advection limit.

In this paper, we modify the analysis described above in one additional way. Instead of using the usual critical balance parameter χ as the argument of g(χ), we use

$$\begin{equation} \chi _{\rm eff} \, = \, \frac{\chi }{\sqrt{1 + \chi _{0}^{2}}}, \end{equation} \tag{ C17 }$$

where χ₀ is defined in Equation (37). In the strong turbulence regime (χ₀ ≪ 1) this modification makes no difference. In the weak turbulence regime (χ₀ ≫ 1) this has the effect of extending the "filled" region of the spectrum (i.e., g(χ_eff) ≈ 1) up through all wavenumbers with k_∥ ≲ k_0∥.

Finally, we must adjust the highest frequency part of the spectrum to account self-consistently for the effects of ion cyclotron damping. Because the cyclotron resonance at high k_∥ has a rapid onset (see Figure 11(a)), we need only model its effects over a limited range of wavenumber space. We truncate the calculation of the spectrum at a maximum parallel wavenumber

$$\begin{equation} \frac{k_{\parallel {\rm max}} V_{\rm A}}{\Omega _p} \, = \, \frac{0.72}{\beta ^{0.43}} \end{equation} \tag{ C18 }$$

at which |γ_A/Ω_p| ≈ 1. Above this wavenumber, we found that slowly varying solutions to the Vlasov–Maxwell dispersion relation cease to exist (see also Stix 1992). Between 0.1k_∥max and k_∥max, we include the time-steady effect of resonant damping by assuming that the local Alfvénic wave power is produced solely by the nonlinear coupling with the fast mode. If only the coupling and damping are present, the time-steady transport equation simplifies to

$$\begin{equation} \frac{\partial E_{\rm A}}{\partial t} \, \approx \, \frac{E_{\rm F} - E_{\rm A}}{\tau _{\rm AF}} - 2 \gamma _{\rm A} E_{\rm A} \, = \, 0, \end{equation} \tag{ C19 }$$

which can be solved analytically for E_A. However, we note that we already have a time-steady solution for E_A in the presence of nonlinear coupling, but it does not take into account the ion cyclotron damping. Equation (54) gives that solution, which we now call E_0A. Thus, we insert it in place of E_F in Equation (C19) above, since that is the solution toward which the coupling will drive the system in the absence of damping. We then use the analytic solution

$$\begin{equation} E_{\rm A} \, \approx \, \frac{E_{0 {\rm A}}}{1 + 2 \gamma _{\rm A} \tau _{\rm AF}} \end{equation} \tag{ C20 }$$

to account for ion cyclotron dissipation at high k_∥. This solution gives rise to a significant reduction in the power spectrum when the damping rate γ_A exceeds the rate at which power is supplied from the nonlinear coupling.

This section is conceptually similar to Appendix C.1 in that we ignore both damping and mode coupling, assume time-steady conditions, and model the spectral transport of fast-mode waves as a balance between diffusive cascade and the outer-scale source term. In that limiting case, Equation (32) becomes

$$\begin{equation} \frac{\partial \varepsilon }{\partial k} \, = \, k^{2} S_{\rm F}, \end{equation} \tag{ C21 }$$

where the cascade rate is defined here as

$$\begin{equation} \varepsilon \, = \, - \frac{4\pi \alpha _{\rm F} k^{8} \sin \theta }{V_{\rm A}} E_{\rm F} \frac{\partial E_{\rm F}}{\partial k}. \end{equation} \tag{ C22 }$$

When S_F = 0, the cascade rate is constant along radial rays of constant θ, and thus E_F∝k^−7/2. We follow Chandran (2005) and others by assuming a Gaussian shape for the fast-mode source term, but we also add a corresponding sin θ angle dependence, with

$$\begin{equation} S_{\rm F} (k) \, = \, S_{0} \, \exp \left[ - \left(\frac{k}{k_{0 {\rm F}}} \right)^{2} \right] \, \sin \theta. \end{equation} \tag{ C23 }$$

We eventually set S₀ at the level required to maintain the spectrum at the known total energy density U_F. Our choice to constrain k_0F to be equal to k_0⊥ = 1/λ_⊥ is discussed in Section 3.1.

Because both the left and right sides of Equation (C21) depend identically on sin θ, the resulting time-steady solution for E_F(k) becomes independent of θ. This outcome was motivated by simulation results that show that the fast-mode power spectrum is largely isotropic in wavenumber space (Cho & Lazarian 2003; Svidzinski et al. 2009). It is also possible that additional isotropization of the fast-mode spectrum can come from couplings with slow-mode waves (Chandran 2008a) or from multi-scale "wandering" of the magnetic field that gives rise to a continuously varying direction for θ = 0 (Shalchi & Kourakis 2007). In future work, our method of artificially forcing isotropy via the source term should be replaced with a more realistic description.

In any case, we cancel out both instances of sin θ and integrate Equation (C21) to obtain

$$\begin{equation} \varepsilon (k) \, = \, k_{0 {\rm F}}^{3} S_{0} \left(\frac{\sqrt{\pi }}{4} \mbox{erf} \, x - \frac{x e^{-x^{2}}}{2} \right), \end{equation} \tag{ C24 }$$

where here x = k/k_0F. In the limit x ≫ 1, the term in parentheses above approaches a constant value of $\sqrt{\pi } / 4$. In the limit x ≪ 1, the term in parentheses is approximately equal to x³/3. Finally, we integrate the definition of ε to get the time-steady spectrum,

$$\begin{equation} \frac{E_{\rm F}^{2}}{2} \, = \, \frac{V_{\rm A}}{4\pi \alpha _{\rm F}} \int _{k}^{\infty } d\kappa \, \frac{\varepsilon (\kappa)}{\kappa ^8}, \end{equation} \tag{ C25 }$$

which we evaluate numerically using Equation (C24) for the cascade rate in the integrand. In the energy containing range (k ≪ k_0F), this yields E_F∝k⁻², and thus v_k∝k^+1/2.

If we consider high values of k above those affected by the outer-scale source term, we can solve for the transition from the inertial range to the dissipation range in the fast-mode power spectrum. An analytic solution becomes possible if we rewrite the fast-mode spectral transport time τ_F as a function of k and θ only. This can be done by using the time-steady inertial range scaling for v_k, with

$$\begin{equation} v_{k}^{2} \, = \, v_{0}^{2} \left(\frac{k}{k_0} \right)^{-1/2}, \,\,\,\, \tau _{\rm F} \, = \frac{V_{\rm A}}{v_{0}^{2} \sin \theta \sqrt{k \, k_0}}. \end{equation} \tag{ C26 }$$

Note that Equation (8) of Suzuki et al. (2007) gave the same result for the fast-mode cascade timescale (but without the sin θ term). The normalization wavenumber k₀ is defined arbitrarily here; it needs to be set well below the regime of strong damping, but well above the outer-scale wavenumber k_0F so that we can justify ignoring the source term.

The above approximation gives D_F∝k^5/2sin θ. It is also straightforward to model the fast-mode damping rate as being proportional to a constant power of the wavenumber, and thus we assume γ_F = γ₀(k/k₀)^z. Note that the normalizing constant γ₀ may depend on the angle θ as well. The time-steady version of Equation (32) becomes

$$\begin{equation} \frac{\alpha _{\rm F} v_{0}^{2} k_{0}^{1/2} \sin \theta }{k^{2} V_{\rm A}} \frac{\partial }{\partial k} \left(k^{9/2} \frac{\partial E_{\rm F}}{\partial k} \right) \, = \, 2 \gamma _{0} \left(\frac{k}{k_0} \right)^{z} E_{\rm F}. \end{equation} \tag{ C27 }$$

Defining the auxiliary variable

$$\begin{equation} x \, = \, \frac{2}{7} \left(\frac{k_0}{k} \right)^{7/2} \end{equation} \tag{ C28 }$$

helps to greatly simplify the differential equation. Thus,

$$\begin{equation} \frac{\partial ^{2} E_{\rm F}}{\partial x^2} \, = \, \frac{c_{\gamma } E_{\rm F}}{x^{(2z + 13)/7}}, \end{equation} \tag{ C29 }$$

where

$$\begin{equation} c_{\gamma } \, = \, \frac{2 \gamma _{0} V_{\rm A}}{\alpha _{\rm F} v_{0}^{2} k_{0} \sin \theta } \left(\frac{2}{7} \right)^{(2z + 13)/7} \end{equation} \tag{ C30 }$$

is a constant that is essentially the ratio of the damping rate to the cascade rate at the normalization wavenumber k₀.

For z ⩾ 1, Equation (C29) is solved analytically with two linearly independent terms proportional to the two types of modified Bessel function (I_n and K_n). Knowing that the only physically realistic solution is one that decreases monotonically with increasing k (or with decreasing x), we then use only one of those terms, which is given by

$$\begin{equation} E_{\rm F}(k) \, \propto \, k^{-7/4} \, K_{\zeta } \left[ 2\zeta \left(\frac{2}{7} \right)^{1/(2\zeta)} \sqrt{c_{\gamma }} \left(\frac{k}{k_0} \right)^{7/(4\zeta)} \right], \end{equation} \tag{ C31 }$$

where ζ = 7/(2z − 1). Note that the transit-time damping rate of Equation (58) gives z = 1 and ζ = 7. In the limiting case that the modified Bessel function of the second kind has a small argument, we have K_ζ(x) ∼ x^−ζ and thus E_F∝k^−7/2, independent of the value of ζ. This is the proper inertial-range solution in the case of either low wavenumber (k ≪ k₀) or weak damping (c_γ ≪ 1). The opposite case of a large argument gives exponentially steep dissipation in the limit of large k and/or large γ₀. This kind of solution was also derived by Hunana & Zank (2010).

Another useful special case for the damping exponent is z = 1/2. For this value of the exponent, Equation (C29) is solved with two linearly independent power-law terms. As above, we keep only the solution that does not diverge as x → 0 (i.e., as k → ∞), and the time-steady spectrum is given by

$$\begin{equation} E_{\rm F}(k) \, \propto \, k^{-7(1 + \sqrt{1 + 4 c_{\gamma }})/4}. \end{equation} \tag{ C32 }$$

The weak-damping limit of c_γ ≪ 1 gives the proper inertial-range solution E_F∝k^−7/2, but the presence of a nonzero value of c_γ makes the spectrum steeper. This is one (possibly rare) case in which a physically motivated source of damping gives rise to a power-law "dissipation range."