ALFVÉN WAVE REFLECTION AND TURBULENT HEATING IN THE SOLAR WIND FROM 1 SOLAR RADIUS TO 1 AU: AN ANALYTICAL TREATMENT

Beginning with the discovery of large-amplitude Alfvén waves in the interplanetary medium four decades ago (Belcher et al. 1969; Belcher & Davis 1971), spacecraft measurements have demonstrated that turbulent fluctuations in the magnetic field and plasma velocity are ubiquitous in the solar wind (see, e.g., Tu & Marsch 1995; Goldstein et al. 1995; Bruno & Carbone 2005). As in hydrodynamic turbulence, fluctuation energy in turbulent magnetized plasmas cascades from large scales to small scales. Upon reaching sufficiently small scales (typically, of order the proton gyroradius and smaller), the energy is dissipated, thereby heating the ambient plasma (see Schekochihin et al. 2009 for a recent and detailed discussion). Turbulent heating of the solar wind is thus inevitable, and may play an important role in the origin and evolution of the solar wind (see, e.g., Coleman 1968; Barnes 1981; Hollweg 1983; Tu et al. 1984; Hollweg 1986; Tu 1987, 1988; Hollweg & Johnson 1988; Velli et al. 1989, 1990; Grappin & Velli 1991; Grappin et al. 1993; Marsch & Tu 1996; Matthaeus et al. 1999).

The bulk of the fluctuation energy in solar-wind turbulence is at length scales greatly exceeding the proton gyroradius, is non-compressive, and to a good approximation can be treated within the approximation of reduced magnetohydrodynamics (Kadomtsev & Pogutse 1974; Strauss 1976; Zank & Matthaeus 1992; Schekochihin et al. 2009). The cascade of energy to smaller scales in reduced magnetohydrodynamic (MHD) turbulence arises only from interactions between oppositely propagating Alfvén waves. It has thus been known for some time that the launching of outward-propagating Alfvén waves by the Sun is not sufficient to cause turbulent heating, unless some additional process generates Alfvén waves that propagate toward the Sun (in the solar-wind frame). Different mechanisms for producing Sunward-propagating waves have been considered, including the parametric instability (Galeev & Oraevskii 1963; Viñas & Goldstein 1991; Del Zanna et al. 2001) and instabilities driven by velocity shear (Roberts et al. 1992; Breech et al., 2008). In this paper, we focus on a third mechanism: the non-WKB reflection of Alfvén waves.

When the non-uniformity of the background plasma is taken into account, purely outward-propagating waves are no longer a solution to the reduced MHD equations. Instead, outward and inward waves are coupled through linear terms that are proportional to the gradients in the Alfvén speed (Velli 1993; Hollweg & Isenberg 2007). This coupling gives rise to wave reflection. At sufficiently high frequencies, linear Alfvén waves propagating in a hypothetical, steady solar wind are insensitive to radial variations in the Alfvén speed and undergo little reflection. On the other hand, for wave periods exceeding ∼20 minutes, the transmission coefficient for outward-propagating Alfvén waves launched from the Sun drops significantly below unity (Velli 1993). A number of authors (e.g., Velli et al. 1989; Matthaeus et al. 1999; Dmitruk et al. 2001; Dmitruk & Matthaeus 2003; Chandran et al. 2009; Verdini et al. 2009) have discussed the propagation of Alfvénic fluctuations originating at the Sun, emphasizing the self-consistent generation of reflected waves and the resulting cascade and dissipation of wave energy. In particular, Dmitruk et al. (2002, hereafter D02) developed a simple phenomenological model describing reflections and turbulence and obtained an analytical approximation for the turbulent heating rate valid close to the Sun in the limit of strong turbulent dissipation. Elements of their approach were incorporated into detailed solar wind models by Cranmer & van Ballegooijen (2005, hereafter CvB05), Cranmer et al. (2007), and Verdini & Velli (2007, hereafter VV07).

In view of the potential usefulness of D02's analytical approximation, our purpose here is to extend their approximation out to the Alfvén critical point and beyond, i.e., throughout the region where the principal solar wind heating and acceleration occur, and to obtain analytical formulas for the radial profiles of the turbulent heating rate and wave amplitudes. In Section 2, we focus on waves with periods of hours and longer. In Section 3, we extend our model to account for shorter-period waves.

We begin with the usual equations of ideal MHD,

$$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\bf v}) = 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \rho \left(\frac{\partial }{\partial t} + {\bf v} \cdot \nabla \right) {\bf v} = - \nabla p_{\rm tot} + \frac{{\bf B} \cdot \nabla {\bf B}}{4\pi } - \rho \nabla \Phi, \end{equation} \tag{ 2 }$$

and

$$\begin{equation} \frac{\partial {\bf B}}{\partial t} = \nabla \times ({\bf v} \times {\bf B}), \end{equation} \tag{ 3 }$$

where ρ, v, and B are the mass density, velocity, and magnetic field, Φ is the gravitational potential, p_tot = p + B²/8π is the total pressure, and p is the plasma pressure. We set

$$\begin{equation} {\bf v} = {\bf U} + \delta {\bf v} \end{equation} \tag{ 4 }$$

and

$$\begin{equation} {\bf B} = {\bf B}_0 + \delta {\bf B}, \end{equation} \tag{ 5 }$$

neglect the Sun's rotation, and take the background flow velocity U to be aligned with the background magnetic field B₀. We assume that δv and δB are perpendicular to B₀ and non-compressive. We neglect density fluctuations¹ and take ρ, U, and B₀ to be steady-state solutions of Equations (1) through (3). The vector Alfvén velocity is given by

$$\begin{equation} {\bf v}_{\rm A} = \frac{{\bf B}_0}{\sqrt{4\pi \rho }}, \end{equation} \tag{ 6 }$$

and the Elsässer variables are defined as

$$\begin{equation} {\bf z}^\pm = \delta {\bf v} \mp \delta {\bf b}, \end{equation} \tag{ 7 }$$

where $\delta {\bf b} = \delta {\bf B}/\sqrt{4\pi \rho }$. Re-writing Equations (2) and (3) in terms of z^±, we obtain (Velli 1993; VV07)

$$\frac{\partial {\bf z}^\pm }{\partial t} +\left({\bf U} \pm {\bf v}_{\rm A}\right)\cdot \nabla {\bf z}^\pm + {\bf z}^\mp \cdot \nabla \left({\bf U} \mp {\bf v}_A\right)$$

$$\begin{equation} + \frac{1}{2} \left({\bf z}^- - {\bf z}^+\right) \left(\nabla \cdot {\bf v}_{\rm A} \mp \frac{1}{2}\nabla \cdot {\bf U} \right) = - \!\left({\bf z}^\mp \cdot \nabla {\bf z}^\pm + \frac{\nabla p_{\rm tot}}{\rho }\right). \end{equation} \tag{ 8 }$$

The ∇p_tot term cancels the compressible part of the remaining terms in Equation (8) to maintain the incompressibility condition ∇ · z^± = 0.

To proceed further, we assume that the background magnetic field B₀ possesses a field line that is purely radial. This field line can be at any heliographic latitude. We then consider a cylindrical coordinate system (R, ϕ, z) whose z-axis is aligned with this radial field line, and we restrict our analysis to the region close to the z-axis, in which we assume that v_Aϕ = U_ϕ = 0 and ∂U/∂ϕ = ∂v_A/∂ϕ = 0. Near the z-axis, the components of U and v_A can be expanded in powers of R,

$$\begin{eqnarray} U_R & = & a_1(z) R + \cdots, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} v_{{\rm A}, R} & = & a_2(z) R + \cdots. \end{eqnarray} \tag{ 10 }$$

Keeping just the first terms in these Taylor series, we find that ∂U_R/∂R = U_R/R, ∂v_A,R/∂R = v_A,R/R, and

$$\begin{equation} {\bf z}^\mp \cdot \nabla \left({\bf U} \mp {\bf v}_{\rm A}\right) = \frac{{\bf z}^\mp (U_R \mp v_{{\rm A},R})}{R}. \end{equation} \tag{ 11 }$$

We define

$$\begin{equation} \frac{\partial }{\partial s} = ({\bf B}_0/B_0) \cdot \nabla, \end{equation} \tag{ 12 }$$

so that s is the distance along a background magnetic field line. Since $\nabla R = \hat{R}$,

$$\begin{equation} \frac{U_R \mp v_{{\rm A},R}}{R} = \frac{(U\mp v_{\rm A})}{R} \;\frac{\partial R(s)}{\partial s}, \end{equation} \tag{ 13 }$$

where R(s) is to be interpreted as the distance from the z-axis to a given background magnetic field line.

We define vector versions of the Heinemann & Olbert (1980) variables f and g as follows:

$$\begin{eqnarray} {\bf z}^+ & = & \frac{{\bf g} \,\eta ^{1/4}}{1 + \eta ^{1/2}} \end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray} {\bf z}^- & = & \frac{{\bf f} \,\eta ^{1/4}}{1 - \eta ^{1/2}}, \end{eqnarray} \tag{ 15 }$$

where

$$\begin{equation} \eta = \rho /\rho _a, \end{equation} \tag{ 16 }$$

and ρ_a is the value of ρ at the Alfvén critical point (on the z-axis). Mass conservation requires that ∂(ρU/B₀)/∂s = 0, and thus η^1/2 = v_A/U. With the use of Equations (11) and (13), Equation (8) can be re-written as (VV07)²

$$\begin{eqnarray} &&\dsty \frac{\partial {\bf g}}{\partial t} + (U + v_{\rm A}) \frac{\partial {\bf g}}{\partial s} + (U+v_{\rm A}) {\bf f} \frac{\partial L}{\partial s}\nonumber\\ &&\quad\dsty = {-} \!\left(\frac{1+\eta ^{1/2}}{\eta ^{1/4}}\right) \left({\bf z}^- \cdot \nabla {\bf z}^+ + \frac{\nabla p_{\rm tot}}{\rho }\right) \end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray} &&\dsty \frac{\partial {\bf f}}{\partial t} + (U - v_{\rm A}) \frac{\partial {\bf f}}{\partial s} + (U-v_{\rm A}) {\bf g} \frac{\partial L}{\partial s}\nonumber\\ &&\quad\dsty= \left(\frac{\eta ^{1/2} - 1}{\eta ^{1/4}}\right) \left({\bf z}^+ \cdot \nabla {\bf z}^- + \frac{\nabla p_{\rm tot}}{\rho }\right), \end{eqnarray} \tag{ 18 }$$

where L = ln [R(s)η^1/4/R_ref], and R_ref is an arbitrary constant. Close to the z-axis, ∂(B₀R²)/∂s = 0, and R_ref can be chosen so that

$$\begin{equation} L = - \frac{1}{2} \ln \left(\frac{v_{\rm A}}{v_{{\rm A}a}}\right), \end{equation} \tag{ 19 }$$

where v_Aa is the Alfvén speed at the Alfvén critical point on the z-axis.

Dmitruk et al. (2002) numerically solved Equation (8) in the limit U → 0. They found that the numerical results could to a large extent be understood in terms of a phenomenological model in which z^∓ · ∇z^± is replaced by |z^∓|z^±/2λ(s). We employ the same phenomenological model (see Chandran et al. 2009 for a more detailed discussion), replacing the nonlinear terms z^∓ · ∇z^± + ρ⁻¹∇p_tot in Equations (17) and (18) with |z^∓|z^±/2λ(s). Equations (17) and (18) then become

$$\begin{equation} \frac{\partial {\bf g}}{\partial t} + (U + v_{\rm A}) \frac{\partial {\bf g}}{\partial s} + (U+v_{\rm A}) {\bf f} \frac{\partial L}{\partial s} = - \frac{\eta ^{1/4}}{2\lambda } \left|\frac{f}{1 - \eta ^{1/2}}\right| {\bf g} \end{equation} \tag{ 20 }$$

and

$$\begin{equation} \frac{\partial {\bf f}}{\partial t} + (U - v_{\rm A}) \frac{\partial {\bf f}}{\partial s} + (U-v_{\rm A}) {\bf g} \frac{\partial L}{\partial s} = - \frac{\eta ^{1/4}|g|{\bf f}}{2\lambda (1 + \eta ^{1/2})}. \end{equation} \tag{ 21 }$$

Observations (e.g., Belcher & Davis 1971) show that most of the Alfvénic power in the solar wind is at long periods, of the order of hours, as seen in the spacecraft frame. In this section, we thus restrict our attention to waves with periods of hours and longer in the frame of the Sun. We further assume that the amplitude of the outward-propagating waves is sufficiently large, and λ is sufficiently small, that the timescale for the cascade and dissipation of z⁻ fluctuations, λ/z⁺, is shorter than the wave period. In this case, we can neglect the first and second terms on the left-hand side of Equation (21) and solve Equation (21) to obtain

$$\begin{equation} {\bf f} = \frac{ 2\lambda (\eta - 1) U}{\eta ^{1/4}} \;\frac{\partial L}{\partial s} \;\hat{g}, \end{equation} \tag{ 22 }$$

where $\hat{g} = {\bf g}/g$. We assume that the rapid cascading of z⁻ energy leads to the inequality

$$\begin{equation} z^- \ll z^+, \end{equation} \tag{ 23 }$$

so that we can neglect the third term on the left-hand side of Equation (20). Taking the dot product of Equation (20) with g, and then taking the time average of the resulting equation, we find that

$$\begin{equation} \frac{\partial }{\partial s} \langle g^2 \rangle = - 2 \langle g^2 \rangle \left| \frac{\partial L}{\partial s}\right|, \end{equation} \tag{ 24 }$$

where 〈 ⋅ ⋅ ⋅ 〉 denotes a time average. Since we restrict our analysis to the vicinity of our z-axis, s ≃ r, where r is the radial coordinate in spherical coordinates centered on the Sun. We take L to be a decreasing function of r near the Sun, to have a minimum at r = r_m (corresponding to the maximum of v_A), and to increase with r at r>r_m. We take r_m to be less than the radius r_a of the Alfvén critical point. Integrating Equation (24), we find that

$$\begin{equation} \langle g^2 \rangle ^{1/2} = g_a \left(\frac{v_{\rm A}}{v_{{\rm A}a}}\right)^{1/2} \end{equation} \tag{ 25 }$$

for r>r_m, where g_a is the rms value of g at r = r_a. For r < r_m, we find that

$$\begin{equation} \langle g^2 \rangle ^{1/2} = g_a \left(\frac{v_{{\rm A}m}^2}{v_{{\rm A}a}v_{\rm A}}\right)^{1/2}, \end{equation} \tag{ 26 }$$

where v_Am is the Alfvén speed at r = r_m.

The energy density of the Alfvénic fluctuations is ρ[(z⁺)² + (z⁻)²]/4. We assume that the energy that is drained from z⁺ and z⁻ fluctuations by the nonlinear "damping" terms in our model is converted to thermal energy. The turbulent heating rate, for z⁻ ≪ z⁺, is thus

$$\begin{equation} Q = \frac{\rho |z^-| (z^+)^2}{4\lambda }. \end{equation} \tag{ 27 }$$

With the use of Equations (22), (25), and (26), we find that the average heating rate is

$$\begin{equation} Q = \frac{\rho \eta ^{1/2} g_a^2 U}{4 (1 + \eta ^{1/2}) v_{{\rm A }a}} \; \left| \frac{dv_{\rm A}}{dr}\right| \end{equation} \tag{ 28 }$$

for r>r_m, and

$$\begin{equation} Q = \frac{\rho \eta ^{1/2} g_a^2 U}{4 (1 + \eta ^{1/2})}\left(\frac{v_{{\rm A}m}^2}{v_{{\rm A}a} v_{\rm A}^2}\right) \frac{dv_{\rm A}}{dr} \end{equation} \tag{ 29 }$$

for r < r_m.

To check the validity of some of our approximations and compare to observations, we consider a model solar wind in which the proton number density is given by Equation (4) of Feldman et al. (1997), which describes coronal holes out to several solar radii, plus an additional r⁻² component:

$$\begin{equation} n(r) = \left(\frac{3.23 \times 10^8}{x^{15.6}} + \frac{2.51 \times 10^6}{x^{3.76}} + \frac{1.85 \times 10^5}{x^2}\right)\; {\rm cm}^{-3}, \end{equation} \tag{ 30 }$$

where x = r/R_☉ and R_☉ = 1 solar radius. This leads to n = 4 cm⁻³ at 1 AU. We take the magnetic field strength to be (Hollweg & Isenberg 2002)

$$\begin{equation} B_0 = \left[\frac{1.5(f_{\rm max} - 1)}{x^6} + \frac{1.5}{x^2}\right] \;{\rm Gauss}, \end{equation} \tag{ 31 }$$

where f_max is the usual super-radial expansion factor, which we set equal to 5. The solar wind speed U is given by flux conservation

$$\begin{equation} U = 9.25 \times 10^{12} \;\frac{\tilde{B}}{\tilde{n}} \;\mbox{cm}\;\mbox{s}^{-1}, \end{equation} \tag{ 32 }$$

where $\tilde{B}$ is B₀ in Gauss, and $\tilde{n}$ is n in units of cm⁻³. Equation (32) implies a flow speed of 750 km s⁻¹ and a proton flux of 3 × 10⁸ cm⁻² s⁻¹ at 1 AU. The Alfvén critical point in this model is at r_a = 11.1 R_☉, and the maximum of v_A is at r_m = 1.60 R_☉. The solar wind speed and Alfvén speed are plotted in Figure 1.

Zoom In Zoom Out Reset image size

We set

$$\begin{equation} \lambda = \lambda _0 \left[\frac{B(R_{\odot })}{B}\right]^{1/2}, \end{equation} \tag{ 33 }$$

where B(R_☉) = 7.5 G from Equation (31), and λ₀ is a constant. Previous studies suggest a range of possible values for λ₀. One choice is to set λ₀ equal to the mean spacing between photospheric flux tubes. Following Spruit (1981), we take this mean spacing to be $R_p \sqrt{\pi B_p/B_{\rm avg}}$, where R_p = 200 km and B_p = 1500 G are the typical radius and magnetic field strength of a photospheric flux tube, and B_avg is the average magnetic field at the coronal base above the flux-tube merging height. If we set B_avg equal to the value of B₀ in Equation (31) as x → 1 (which is 7.5 G), then λ₀ = 5000 km. On the other hand, D02 used λ₀ = 3 × 10⁴ km (comparable to the average size of supergranules), CvB05 considered values of λ₀ in the range 1–30 × 10³ km, and the quantity corresponding to our λ₀ in VV07 was 1.7 × 10⁴ km. We will consider several values of λ₀ in the range spanned by the above values.

To check whether we are justified in neglecting the terms that we have dropped in Equations (20) and (21), it would be desirable to compare our approximate solutions to exact, time-dependent solutions of Equations (20) and (21). However, such time-dependent solutions are beyond the scope of this study. Instead, we numerically solve the time-independent versions of Equations (20) and (21), taking f = 0 at the Alfvén critical point, setting λ₀ = 5000 km and $g_a = 8 \times 10^7\;\rm cm\;s^{-1}$, and assuming that the directions of f and g are the same and do not change with s. This steady-state numerical solution is still a useful reference point since we are considering only long-period oscillations in this section. In Figure 2, we plot the numerical solutions for |f| and g, as well as our approximate analytical solutions from Equations (22), (25), and (26). For r>2 R_☉, our approximate solutions are close to the numerical solutions. Closer to the Sun, our analytical expressions overestimate the values of |f| and g. In particular, our approximate value of g is 23% larger than the numerical solution at r = R_☉. If λ₀ is increased to 10⁴ km (1.5 × 10⁴ km), the approximate solution is 43% (57%) larger than the numerical solution at r = R_☉, although the approximate solutions remain reasonably accurate at r>2 R_☉. Increasing λ₀ reduces the nonlinear terms in Equations (20) and (21) and degrades the accuracy of our approximations.

Zoom In Zoom Out Reset image size

The loss of accuracy in our analytical approximations as λ₀ increases is further illustrated in Figure 3, which is analogous to Figure 5 of D02. The dashed line in this figure shows Q/ρ from Equations (28) and (29). The three solid lines plot the total heating rate Q_total = ρ[(z⁺)²|z⁻| + (z⁻)²|z⁺|]/4λ in the numerical solutions of the time-independent versions of Equations (20) and (21) for $g_{\rm a} = 8.0 \times 10^7 \;\rm cm\;s^{-1}$ and for three different values of λ₀. Going from the top of the figure to the bottom, these values are 2500 km, 5000 km, and 10⁴ km. At r ⩾ 2 R_☉, Equation (28) overestimates the numerical heating rate by ≲50% for the chosen values of λ₀. At r < 1.6 R_☉, Equation (29) overestimates the numerical heating rate by almost an order of magnitude in the worst case (λ₀ = 10⁴ km).

Zoom In Zoom Out Reset image size

In Figure 4, we use our time-independent solutions of Equations (20) and (21) with $g_{\rm a} = 8.0 \times 10^7 \;\rm cm\;s^{-1}$ and λ₀ = 5000 km to plot each term in Equation (21), divided by (1 − η^1/2) to make the plot easier to read. Although we dropped ∂f/∂t to obtain these time-independent solutions, we estimate the approximate magnitude of the ∂f/∂t term by multiplying our numerical solution for f by 2π/P, where P is the wave period, which we take to be three hours. Our estimate of the ∂f/∂t term is smaller than the terms that we have kept inside the Alfvén critical point, but ∂f/∂t increases relative to the dominant terms as r increases toward r_a. We have compared the different terms in Equation (21) at larger r, and find that our neglect of the ∂f/∂t term becomes unjustified for r ≳ 20 R_☉, because the nonlinear timescale becomes larger than the wave period (measured in the Sun's frame) due to the increase in λ and the decrease in z⁺ as r increases. Figure 4 shows that the ∂f/∂s term, which we have neglected in Equation (21), is negligible except in a small region near r = r_m, where our solution for f passes through zero. Although it is not plotted, we have also carried out the same comparison for the different terms in Equation (20), using our time-independent numerical solution. We find that the ∂L/∂s term that we neglect in Equation (20) is indeed smaller than the other terms out to well beyond 20 R_☉, because f ≪ g.

Zoom In Zoom Out Reset image size

In Figure 5, we plot the turbulent heating rate in Equations (28) and (29) (thin solid line). Although the approximations leading to Equations (28) and (29) break down at r ≳ 20 R_☉, we plot Equation (28) out to beyond 100 R_☉ to illustrate how the analytical formulas extrapolate to larger radii. The heating rate in our model depends on a single free parameter, g_a, which we set equal to $8.0 \times 10^7 \; \rm cm\;s^{-1}$ in order to roughly match two previous empirically constrained models of plasma heating in the fast solar wind (Allen et al. 1998; Cranmer et al. 2009). We note that for z⁻ ≪ z⁺, this value of g_a corresponds to a fluctuating velocity of 200 km s⁻¹ at the Alfvén critical point. The dashed line in Figure 5 is the total (electron plus proton) heating rate Q_A98 in model SW2 of Allen et al. (1998). The thick solid line is the proton heating rate in the model of Cranmer et al. (2009) (their Equation (14), which applies at r>0.3 AU), multiplied by 4/3 to convert the proton heating rate to an approximate total heating rate, denoted as Q_C09. (The proton heating rate is typically 60%–90% of the total heating rate in Cranmer et al.'s 2009 model).

Zoom In Zoom Out Reset image size

Model SW2 of Allen et al. (1998) does not provide a tight constraint on the heating rate, since widely different heating functions and momentum-deposition profiles (their models SW2, SW3, and SW4) were shown to be consistent with in situ observations of the fast wind. We include Q_A98 in Figure 5 to show that our results for the heating rate at 2 R_☉ < r < 7 R_☉ are reasonably close to at least one empirically constrained model for a plausible value of g_a. On the other hand, our heating rate greatly exceeds Q_A98 at r < 1.5 R_☉. This disagreement may be due in part to the effects described in Section 4 and also to the point illustrated in Figure 3 that we substantially overestimate the heating rate close to the Sun when λ₀>5000 km. We emphasize that although Figure 5 represents a rough first attempt to compare our model to observational constraints on the heating rate, this comparison must be viewed with some caution, because the density profiles in the three models are different, and Q depends upon ρ in Equations (28) and (29).

In Figure 6, we plot the rms amplitude of the fluctuating velocity in our low-wave-frequency model as a function of r (dashed line), using Equations (25) and (26) and the relation δv_rms ≃ 〈(z⁺)²〉^1/2/2, which is valid when z⁻ ≪ z⁺. We also plot lower and upper limits from Ultraviolet-Coronagraph-Spectrometer (UVCS) observations (Esser et al. 1999; CvB05) and in situ measurements from Helios and Ulysses (Bavassano et al. 2000). (The solid line in this figure is from the extended model described in Section 3.) For $g_a = 8 \times 10^7\; \rm cm\;s^{-1}$, our low-wave-frequency model overpredicts the rms amplitude of velocity fluctuations at 1.5 R_☉ < r < 2 R_☉. As noted above, the approximations leading to Equations (25) and (26) break down at r ≳ 20 R_☉. Nevertheless, we plot Equation (25) out to larger r to illustrate how the analytical formulas extrapolate to larger radii. This extrapolation shows that for $g_a = 8 \times 10^7 \; \rm cm\;s^{-1}$, Equations (25) and (26) underpredict δv_rms at r>70 R_☉.

Zoom In Zoom Out Reset image size

In Section 2, we took the waves launched by the Sun to have periods of hours and longer. On the other hand, CvB05 suggested that a significant fraction of the outward wave power near the Sun is in waves with periods of 1 minute to 1 hr.³ In this section, we extend our model to account for such shorter-period waves in an approximate manner. The approximations we use are very simplistic, but are not limited to r ≲ 20 R_☉.

In the spirit of our simple, phenomenological modeling, we begin by replacing each g and f in Equations (20) and (21) with the (time-independent) rms values 〈g²〉^1/2 and 〈f²〉^1/2, respectively. As shown by Velli (1993) and others, wave reflection becomes less efficient as the wave period decreases below ∼1 hr, since higher-frequency waves with shorter parallel wavelengths are less affected by the radial gradient in the Alfvén speed. To account for this effect, we set

$$\begin{equation} \frac{\partial L}{\partial s} \rightarrow \chi \frac{\partial L}{\partial s} \end{equation} \tag{ 34 }$$

in Equations (20) and (21), where χ is a dimensionless constant. In the low-wave-frequency limit of Section 2, χ → 1. To model the mix of wave frequencies launched from the Sun, we leave χ as a free parameter between 0 and 1. We take 〈g²〉^1/2 and 〈f²〉^1/2 to have scale lengths in the direction of B₀ of order r. We continue to assume that f ≪ g, and set (η − 1)∂L/∂s → |(η − 1)∂L/∂s| in Equation (21) to maintain the positivity of 〈f²〉, so that Equation (21) becomes

$$\begin{equation} \langle f^2\rangle ^{1/2} = \frac{ 2\chi \lambda |\eta - 1| U}{\eta ^{1/4}} \;\left|\frac{\partial L}{\partial s}\right|. \end{equation} \tag{ 35 }$$

Dropping the third term on the left-hand side of Equation (20) and using Equation (35), we obtain

$$\begin{equation} \frac{\partial }{\partial s} \langle g^2 \rangle = - 2\chi \langle g^2 \rangle \left| \frac{\partial L}{\partial s}\right|, \end{equation} \tag{ 36 }$$

which we integrate to obtain

$$\begin{equation} \langle g^2 \rangle ^{1/2} = g_a \left(\frac{v_{\rm A}}{v_{{\rm A}a}}\right)^{\chi /2} \end{equation} \tag{ 37 }$$

for r>r_m and

$$\begin{equation} \langle g^2 \rangle ^{1/2} = g_a \left(\frac{v_{{\rm A}m}^2}{v_{{\rm A}a}v_{\rm A}}\right)^{\chi /2} \end{equation} \tag{ 38 }$$

for r < r_m. From Equation (27), we then find that

$$\begin{equation} Q = \frac{\chi \rho \eta ^{1/2} g_a^2 U}{4v_{\rm A} (1 + \eta ^{1/2})} \left(\frac{v_{\rm A}}{v_{{\rm A}a}}\right)^\chi \; \left|\frac{dv_{\rm A}}{dr}\right| \end{equation} \tag{ 39 }$$

for r>r_m, and

$$\begin{equation} Q = \frac{\chi \rho \eta ^{1/2} g_a^2 U}{4v_{\rm A} (1 + \eta ^{1/2})}\left(\frac{v_{{\rm A}m}^2}{v_{{\rm A}a} v_{\rm A}}\right)^\chi \frac{dv_{\rm A}}{dr} \end{equation} \tag{ 40 }$$

for r < r_m. The results for 〈f²〉^1/2, 〈g²〉^1/2 and Q in our low-frequency model of Section 2 are recovered by setting χ = 1 in Equations (35) through (40).

When z⁻  ≪   z⁺, the rms amplitude of the fluctuating velocity δv_rms is approximately 〈(z⁺)²〉^1/2/2, and there are two free parameters that determine δv_rms: g_a and χ. We set g_a = 7.2 × 10⁷ cm s⁻¹ and χ = 0.65 in order to match the observational constraints on the fluctuating velocity, shown in Figure 6. The solid line in Figure 6 is the resulting value of δv_rms obtained from Equations (37) and (38). In Figure 7, we plot the heating rate from Equations (39) and (40) for these values of g_a and χ, along with the empirically constrained heating rates Q_A98 and Q_C09 from the models of Allen et al. (1998) and Cranmer et al. (2009) described in Section 2. We plot the rms amplitudes of z⁺ and z⁻ resulting from Equations (35), (37), and (38) in Figure 8, as well as in situ measurements from Helios and Ulysses (Bavassano et al. 2000). The rms amplitude of z⁻ that results from Equation (35) depends upon both χ and λ₀. To match the in situ data, we set λ₀ = 1.6 × 10⁴ km in this figure, keeping χ = 0.65 as in Figures 6 and 7. As Figure 8 shows, the approximation z⁻ ≪ z⁺ fails at r ≳ 1 AU.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall, the model matches the observations and previous models fairly well. However, as in our low-wave-frequency model, the heating rate in our extended model is much larger than Q_A98 at r < 1.5 R_☉. This discrepancy may arise in part for the reasons discussed in Section 4. In addition, our extended model is likely similar to our low-frequency model in significantly overestimating Q close to the Sun when λ₀ ⩾ 5000 km (Figure 3). A second difference between the heating rates plotted in Figure 7 is that Q becomes smaller than Q_C09 at r>1 AU. This discrepancy may arise because of the increasing importance at large r of sources of inward waves (such as velocity shear) that are not included in our model. We note, however, that some caution is warranted when comparing our heating rate to Q_A98 and Q_C09, because the density profiles in the three models are different, and Q depends upon ρ in Equations (39) and (40).

Compared to our model for low-frequency waves in Section 2, our extended model accounting for shorter-period waves does a notably better job of reproducing the observational constraints on the velocity fluctuations in the solar wind. This is not surprising, given that there is an additional free parameter in the extended model and the approximations underlying our low-wave-frequency model break down at r ≳ 20 R_☉. There is, however, an additional reason for the closer agreement between the extended model and the data. As can be seen in Figure 6, the value of δv_rms in our low-wave-frequency model drops off more rapidly with radius than do the data. When the efficiency of reflection is reduced, the rate at which the outward wave energy cascades and dissipates is reduced, and so |(d/dr)δv_rms| decreases.

If the Sun launches an outward-propagating (z⁺) Alfvén wave packet with a perpendicular scale length λ₀ at the coronal base, with λ₀ greatly exceeding the dissipation scale, then this wave packet must propagate some distance into the corona before its energy cascades to small perpendicular scales and dissipates. We estimate this distance using approximations and results from our analytical models. As in Equation (20), the rate at which z⁻ waves shear z⁺ waves is

$$\begin{equation} \gamma _{\rm nl}^- = \frac{|z^-|}{2\lambda }. \end{equation} \tag{ 41 }$$

The quantity

$$\begin{equation} \Gamma = \int _{r_b}^r \frac{|z^-| dr_1}{2\lambda (U+v_{\rm A})}. \end{equation} \tag{ 42 }$$

is the time integral of γ⁻_nl[r(t), t] evaluated at a position r(t) that moves outward with the z⁺ wave packet at speed U + v_A, starting from the time the wave packet leaves the coronal base at r = r_b ≃ 1 R_☉ until the time the wave packet reaches radius r. Roughly speaking, the z⁺ wave packet must propagate out to a distance at which Γ ≳ 1 before its energy cascades and dissipates. We evaluate the right-hand side of Equation (42) in our low-frequency model using Equation (22). To evaluate Equation (42) using our extended-model results, we replace z⁻ with 〈(z⁻)²〉^1/2 and use Equation (35). The integral in Equation (42) can then be evaluated analytically to yield

$$\begin{equation} \Gamma = \left\lbrace \begin{array}{@{}ll} (\chi /2)\ln (v_{\rm A}/v_{{\rm A}b}) & {\rm if}\;r<r_m \\ \ms (\chi /2)\ln [v_{{\rm A}m}^2/(v_{{\rm A}b}v_{\rm A})] & {\rm if}\;r>r_m, \end{array} \right. \end{equation} \tag{ 43 }$$

where the result for our low-frequency model is obtained by setting χ = 1. The value of Γ is plotted in Figure 9, assuming χ = 1 and using the profiles for n, B₀, and U given in Equations (30) through (32). The condition Γ>1 only holds for r ≳ 7 R_☉, suggesting that the turbulent heating rate may be smaller than in our model within a few R_☉ of the solar surface. We note, however, that the reduction in Q due to the condition Γ < 1 may itself be mitigated by a further complication, that turbulence within the chromosphere (where z⁺ ∼ z⁻ (CvB05)) and photosphere may lead to the launching of Alfvén waves into the corona with a broad range of perpendicular length scales, extending to very small values (Chandran 2008). Further work is needed to explore how these effects modify the radial profile of the turbulent heating rate.

Zoom In Zoom Out Reset image size

In one sense, Figure 9 provides a rough consistency check on models in which the solar wind is heated and accelerated by low-frequency Alfvén-wave turbulence. If Γ were ≪1 at r = r_a, then only a tiny fraction of the Alfvén-wave energy launched from the base of the Sun would dissipate in the region where most of the heating and acceleration of the solar wind take place. A small fraction of the Alfvén-wave energy would be unable to generate the solar wind, because the Alfvén wave energy flux measured in situ at r>0.3 AU is much smaller than the total solar-wind energy flux, even allowing for loss of Alfvén-wave energy due to work done on the flow by the ponderomotive force. On the other hand, if Γ rose from 0 at r = r_b to a value ≫1 at say r = 7 R_☉, then most of the wave energy would dissipate very close to the Sun (within the radius at which Γ reached a value of a few). If the heating rate were to fall too rapidly with r, then the heating profile would become inconsistent with UVCS observations that show ion temperatures staying flat or increasing with radius out to at least 3.5 R_☉ despite adiabatic cooling (Kohl et al. 1998; Antonucci et al. 2000).

In this paper, we consider Alfvén-wave reflection and turbulent heating in the solar wind using the Heinemann–Olbert equations, which describe non-compressive MHD fluctuations in an inhomogeneous medium with a background flow parallel to the background magnetic field. We approximate the nonlinear terms in these equations using a simple phenomenological model from D02. Wave reflection plays an essential role in our calculation of the turbulent heating rate, because only outward-propagating waves are generated by the Sun, and turbulent interactions arise only from interactions between oppositely propagating waves (for the low-frequency, non-compressive Alfvén waves that we consider).

In Section 2, we restrict our attention to waves with periods of hours and longer. Our focus on low-frequency waves in this section is motivated by in situ observations showing that most of the wave power is at periods of hours in the spacecraft frame. We assume that the energy in inward-propagating waves is much less than the energy in outward-propagating waves, i.e., z⁻ ≪ z⁺, an approximation that is appropriate for coronal holes and the fast solar wind at r < 1 AU. We also assume that the wave period is longer than the energy cascade time for Sunward-propagating waves, ∼λ/z⁺. With the use of these approximations, we solve our model equations analytically to obtain expressions for the rms amplitudes of outward and inward waves as a function of r. We also obtain an analytical expression for the radial profile of the turbulent heating rate. There is one free parameter in the heating rate in our low-wave-frequency model: the amplitude of the outward-propagating waves at some fixed reference point. We take this reference point to be the Alfvén critical point r = r_A, which is at 11.1 R_☉ in our model solar wind. Interestingly, as in the model of D02, the heating rate does not depend upon the choice of the dominant perpendicular length scale of the turbulence, λ. Ordinarily, the heating rate is inversely proportional to λ. However, in our model, the amplitude of the inward-propagating waves is proportional to λ, because small λ leads to rapid turbulent dissipation of the inward waves. As a result, Q ≃ ρ〈(z⁺)²〉〈(z⁻)²〉^1/2/4λ becomes independent of λ, as in Equation (5) of D02. Our model can be thought of as generalizing the phenomenological model of D02 by accounting for the solar wind velocity, so that the model can be applied all the way from the coronal base out past the Alfvén critical point.

As discussed in Section 2, one of the approximations in our low-wave-frequency model—that the energy cascade time of Sunward-propagating waves is shorter than the wave period measured in the Sun's frame—breaks down for r ≳ 20 R_☉ for a wave period of ∼3 hr. In addition, our low-wave-frequency model does not account for waves with periods much shorter than one hour, which may make a significant contribution to the outward wave flux from the Sun (CvB05; but see VV07). To overcome these limitations, we develop an "extended" model in Section 3 that approximately accounts for waves with shorter periods and does not require the nonlinear time scale to be shorter than the wave period in the reference frame of the Sun. The main additional ingredient in this extended model is a new free parameter, χ, that models the decrease in the efficiency of wave reflection at shorter wave periods. For appropriate choices of the two free parameters in the model (χ = 0.65 and $g_a = 7.2 \times 10^7\;\rm cm\;s^{-1}$), our extended model is successful at matching observations of velocity fluctuations and the Elsässer fields.

Both of the models we have developed provide analytical expressions for the turbulent heating rate that can be used to incorporate Alfvén-wave reflection and turbulent heating into fluid models of the solar wind, both inside and outside the Alfvén critical point. We conclude by summarizing the relative advantages and disadvantages of the two models:

• 1.  
  The low-wave-frequency model described in Section 2 involves a simplistic phenomenological approximation of only the nonlinear terms in the basic equations. The heating rate in the model involves only a single free parameter, g_a, which determines the amplitude of the outward waves at the Alfvén critical point. On the other hand, one of the assumptions of the model breaks down at r ≳ 20 R_☉, where the nonlinear timescale becomes longer than the wave period measured in the frame of the Sun (for wave periods of ∼3 hr). In addition, the model does not account for waves with periods much shorter than 1 hr, which may contribute significantly to the outward wave flux (CvB05).
• 2.  
  The extended model described in Section 3 does a better job of matching observational constraints on the wave amplitudes, accounts (in a very approximate way) for both long-period waves and waves with periods much less than one hour, and can be applied all the way out to 1 AU. On the other hand, this model involves simple phenomenological modeling of both the nonlinear terms and the linear reflection term. The heating rate in the model also involves two free parameters instead of one.

We thank the anonymous referee for a very helpful report. This work was supported in part by NSF grant ATM-0851005, by DOE grant DE-FG02-07-ER46372, by NSF-DOE grant AST-0613622, and by NASA grants NNX07AP65G and NNX08AH52G.

Note added in proof—a closely related study has been carried out independently by Cranmer (2009).