Parametric Instability, Inverse Cascade, and the $1/f$ Range of Solar-Wind Turbulence

In this paper, weak turbulence theory is used to investigate the nonlinear evolution of the parametric instability in 3D low-$\beta$ plasmas at wavelengths much greater than the ion inertial length under the assumption that slow magnetosonic waves are strongly damped. It is shown analytically that the parametric instability leads to an inverse cascade of Alfv\'en wave quanta, and several exact solutions to the wave kinetic equations are presented. The main results of the paper concern the parametric decay of Alfv\'en waves that initially satisfy $e^+ \gg e^-$, where $e^+$ and $e^-$ are the frequency ($\,f$) spectra of Alfv\'en waves propagating in opposite directions along the magnetic field lines. If $e^+$ initially has a peak frequency $f_0$ (at which $f e^+$ is maximized) and an"infrared"scaling $f^p$ at smaller $f$ with $-1<p<1$, then $e^+$ acquires an $f^{-1}$ scaling throughout a range of frequencies that spreads out in both directions from $f_0$. At the same time, $e^-$ acquires an $f^{-2}$ scaling within this same frequency range. If the plasma parameters and infrared $e^+$ spectrum are chosen to match conditions in the fast solar wind at a heliocentric distance of 0.3 astronomical units (AU), then the nonlinear evolution of the parametric instability leads to an $e^+$ spectrum that matches fast-wind measurements from the Helios spacecraft at 0.3 AU, including the observed $f^{-1}$ scaling at $f \gtrsim 3 \times 10^{-4} \mbox{ Hz}$. The results of this paper suggest that the $f^{-1}$ spectrum seen by Helios in the fast solar wind at $f \gtrsim 3\times 10^{-4} \mbox{ Hz}$ is produced in situ by parametric decay and that the $f^{-1}$ range of $e^+$ extends over an increasingly narrow range of frequencies as $r$ decreases below 0.3 AU. This prediction will be tested by measurements from the Parker Solar Probe.


Introduction
The origin of the solar wind is a long-standing problem (Parker 1958) that continues to receive considerable attention. A leading model for the origin of the fast solar wind appeals to Alfvén waves (AWs) that are launched by photospheric motions. As these AWs propagate away from the Sun, they undergo partial reflection due to the radial variation of the Alfvén speed (Heinemann & Olbert 1980). Nonlinear interactions between counter-propagating AWs then cause AW energy to cascade to small scales and dissipate, heating the plasma (Velli et al. 1989;Zhou & Matthaeus 1989;Cranmer & van Ballegooijen 2005;Verdini et al. 2012;Perez & Chandran 2013;van Ballegooijen & Asgari-Targhi 2017). This heating increases the plasma pressure, which, in conjunction with the wave pressure, accelerates the plasma to high speeds (Suzuki & Inutsuka 2005;Cranmer et al. 2007;Verdini et al. 2010;Chandran et al. 2011;van der Holst et al. 2014).
Although non-compressive AWs are the primary mechanism for energizing the solar wind in this model, a number of considerations indicate that compressive fluctuations have a significant impact on the dynamics of turbulence in the corona and solar wind.
Observations of the tail of Comet-Lovejoy reveal that the background plasma density ρ 0 at r = 1.2R ⊙ (where R ⊙ is the radius of the Sun) varies by a factor of ∼ 6 over distances of a few thousand km measured perpendicular to the background magnetic field B 0 (Raymond et al. 2014). These density variations (denoted δρ) lead to phase mixing of AWs, which transports AW energy to smaller scales measured perpendicular to B 0 (Heyvaerts & Priest 1983). Farther from the Sun, where δρ/ρ 0 is significantly smaller than |δB|/B 0 (Tu & Marsch 1995;Hollweg et al. 2010), AWs still couple to slow magnetosonic waves ("slow waves") through the parametric instability, in which outwardpropagating AWs decay into outward-propagating slow waves and inward-propagating AWs † (Galeev & Oraevskii 1963;Sagdeev & Galeev 1969;Goldstein 1978;Spangler 1986Spangler , 1989Spangler , 1990Hollweg 1994;Dorfman & Carter 2016). This instability and its nonlinear evolution are the focus of the present work.
A number of studies have investigated the parametric instability in the solar wind within the framework of magnetohydrodynamics (MHD) (e.g., Malara et al. 2000;Del Zanna et al. 2001;Shi et al. 2017), while others have gone beyond MHD to account for temperature anisotropy (Tenerani et al. 2017) or kinetic effects such as the Landau damping of slow waves (e.g . Inhester 1990;Vasquez 1995;Araneda et al. 2008;Maneva et al. 2013). Cohen & Dewar (1974), for example, derived the growth rate of the parametric instability in the presence of strong slow-wave damping and randomly phased, parallel-propagating AWs. Terasawa et al. (1986) carried out 1D hybrid simulations and found that Landau damping reduces the growth rate of the parametric instability and that the parametric instability leads to an inverse cascade of AWs to smaller frequencies.
In this paper, weak turbulence theory is used to investigate the nonlinear evolution of the parametric instability assuming a randomly phased collection of AWs at wavelengths much greater than the proton inertial length d i in a low-β plasma, where β is the ratio of plasma pressure to magnetic pressure. The fluctuating fields are taken to depend on all three spatial coordinates, but the wave kinetic equations are integrated over the perpendicular (to B 0 ) wave-vector components, yielding equations for the 1D power spectra that depend only on the parallel wavenumber and time. The starting point of the analysis is the theory of weak compressible MHD turbulence. Collisionless damping of slow waves is incorporated in a very approximate manner analogous to the approach of Cohen & Dewar (1974), by dropping terms containing the slow-wave energy density in the wave kinetic equations that describe the evolution of the AW power spectra.
The remainder of the paper is organized as follows. Section 2 reviews results from the theory of weak compressible MHD turbulence, and Section 3 uses the weak-turbulence wave kinetic equations to recover the results of Cohen & Dewar (1974) in the linear regime. Section 4 shows how the wave kinetic equations imply that AW quanta undergo an inverse cascade towards smaller parallel wavenumbers, and Section 5 presents several exact solutions to the wave kinetic equations. The main results of the paper appear in Section 6, which uses a numerical solution and an approximate analytic solution to the wave kinetic equations to investigate the parametric decay of an initial population of randomly phased AWs propagating in the same direction with negligible initial power in counter-propagating AWs. The numerical results are compared with observations from the Helios spacecraft at a heliocentric distance of 0.3 AU. Section 7 critically revisits the main assumptions of the analysis and the relevance of the analysis to the solar wind. † The terms outward-propagating and inward-propagating refer to the propagation direction in the plasma rest frame. Beyond the Alfvén critical point, all AWs propagate outward in the rest frame of the Sun.
Section 8 summarizes the key findings of the paper, including predictions that will be tested by NASA's Parker Solar Probe.

The Wave Kinetic Equations for Alfvén Waves Undergoing Parametric Decay
In weak turbulence theory, the quantity ω nl /ω linear is treated as a small parameter, where ω nl is the inverse of the timescale on which nonlinear interactions modify the fluctuations, and ω linear is the linear wave frequency. Because the fluctuations can be viewed as waves to a good approximation. The governing equations lead to a hierarchy of equations for the moments of various fluctuating quantities, in which the time derivatives of the second moments (or second-order correlation functions) depend upon the third moments, and the time derivatives of the third moments depend upon the fourth moments, and so on. This system of equations is closed via the randomphase approximation, which allows the fourth-order correlation functions to be expressed as products of second-order correlation functions (see, e.g., Galtier et al. 2000).
The strongest nonlinear interactions in weak MHD turbulence are resonant three-wave interactions. These interactions occur when the frequency and wavenumber of the beat wave produced by two waves is identical to the frequency and wavenumber of some third wave, which enables the beat wave to drive the third wave coherently in time. If the three waves have wavenumbers p, q, and k and frequencies ω p , ω q , and ω k , respectively, then a three-wave resonance requires that and An alternative interpretation of Equations (2.2) and (2.3) arises from viewing the wave fields as a collection of wave quanta at different wavenumbers and frequencies, restricting the frequencies to positive values, and assigning a wave quantum at wavenumber k and frequency ω k the momentum k and energy ω k . Equations (2.2) and (2.3) then correspond to the momentum-conservation and energy-conservation relations that arise when either one wave quantum decays into two new wave quanta or two wave quanta merge to produce a new wave quantum. In the parametric instability in a low-β plasma, a parent AW (or AW quantum) at wavenumber k decays into a slow wave at wavenumber p propagating in the same direction and an AW at wavenumber q propagating in the opposite direction. Regardless of the direction of the wave vector, the group velocity of an AW is either parallel or anti-parallel to the background magnetic field (2.4) and the same is true for slow waves when which is henceforth assumed. At low β slow waves travel along field lines at the sound speed c s , which is roughly β 1/2 times the Alfvén speed v A . Thus, regardless of the perpendicular components of k, p, and q, the frequency-matching condition (Equation (2.3)) for the parametric instability is Combining the z component of Equation (2.2) with Equation (2.6) and taking c s ≪ v A yields and (2.8) Equation (2.8) implies that the frequency |q z v A | of the daughter AW is slightly smaller than the frequency |k z v A | of the parent AW (Sagdeev & Galeev 1969). Thus, the energy of the daughter AW is slightly smaller than the energy of the parent AW. This reduction in AW energy is offset by an increase in slow-wave energy. Chandran (2008) derived the wave kinetic equations for weakly turbulent AWs, slow waves, and fast magnetosonic waves ("fast waves") in the low-β limit. The resulting equations were expanded in powers of β, and only the first two orders in the expansion (proportional to β −1 and β 0 , respectively) were retained. Slow waves are strongly damped in collisionless low-β plasmas (Barnes 1966). Chandran (2008) neglected collisionless damping during the derivation of the wave kinetic equations, but incorporated it afterward in an ad hoc manner by assuming that the slow-wave power spectrum S ± k was small and discarding terms ∝ S ± k unless they were also proportional to β −1 . † (The ± sign in S ± k indicates slow waves propagating parallel (+) or anti-parallel (−) to B 0 .) In the present paper, the wave kinetic equations derived by Chandran (2008) are used to investigate the nonlinear evolution of the parametric instability. It is assumed that slowwave damping is sufficiently strong that all terms ∝ S ± k , even those ∝ β −1 , can be safely discarded. All other types of nonlinear interactions are neglected, including resonant interactions between three AWs, phase mixing, and resonant interactions involving fast waves. Given these approximations, Equation (8) of Chandran (2008) becomes is the 3D wavenumber spectrum of AWs propagating parallel (antiparallel) to B 0 , δ(x) is the Dirac delta function, and the integral over each Cartesian component of p and q extends from −∞ to +∞. The 3D AW power spectra depend upon all three wave-vector components and time. The δ(k − p − q) term enforces the wavenumber-resonance condition (Equation (2.2)), and the δ(q z + k z ) term enforces the frequency-resonance condition (Equation (2.8)) to leading order in β. The integral over the components of p in Equation (2.9) can be carried out immediately, thereby annihilating the first delta function. Equation (2.9) can be further simplified by introducing the 1D wavenumber spectra (2.10) † The one exception to this rule was that Chandran (2008)   . The mathematical expressions next to the arrows represent the contributions to ∂E + (kz2)/∂t from the parametric decay of AWs at kz3, which acts to increase E + (kz2), and the parametric decay of AWs at kz2, which acts to decrease E + (kz2).
In these expressions, and integrating Equation (2.9) over k x and k y , which yields (2.11) Equation (2.11) describes how the 1D (parallel) power spectra E ± evolve and forms the basis for much of the discussion to follow. Given the aforementioned assumptions, the evolution of the 1D power spectra E ± is not influenced by the way that A ± depends on k x and k y . For future reference, the normalization of the power spectra is such that where δv AW and δB AW are the velocity and magnetic-field fluctuations associated with AWs, and . . . indicates an average over space and time (Chandran 2008). Figure 1 offers a way of understanding Equation (2.11). The horizontal color bars in this figure represent the spectra of outward-propagating and inward-propagating AWs, with red representing longer-wavelength waves and violet representing shorter-wavelength waves. AWs propagating in the −B 0 direction at |k z | = k z3 decay into slow waves propagating anti-parallel to B 0 at |k z | ≃ 2k z3 and AWs propagating parallel to B 0 at |k z | = k z2 . AWs propagating parallel to B 0 at |k z | = k z2 decay into slow waves propagating parallel to B 0 at |k z | ≃ 2k z2 and AWs propagating anti-parallel to B 0 at |k z | = k z1 . Equation (2.11) is approximately equivalent to the statement that the rate at which E + 2 = E + (k z2 ) increases via the decay of AWs at |k z | = k z3 is

Physical Interpretation of the Wave Kinetic Equation
, while the rate at which E + 2 decreases via the decay of AWs at and k z1 E − 1 about k z2 in Equation (2.15) thus allows this equation to be rewritten as which is the same as Equation (2.11) to within a factor of order unity.
To be clear, no independent derivation is being presented for Equations (2.13) and (2.14). The foregoing discussion merely points out that Equations (2.13) and (2.14) are equivalent (up to a factor of order unity) to Equation (2.11), which is derived on the basis of weak turbulence theory. It is worth pointing out, however, that several features of Equations (2.13) and (2.14) make sense on a qualitative level. If either E + = 0 or E − = 0, then R 3→2 = R 2→1 = 0, because the parametric instability is a stimulated decay, which ceases if initially all the AWs travel in the same direction. For fixed E + and E − , R 3→2 and R 2→1 vanish as v A → ∞, since the fractional nonlinearities vanish in this limit. Also, R 3→2 and R 2→1 are proportional to β −1/2 (when S ± k is negligibly small, as assumed) because the parametric-decay contribution to ∂A ± k /∂t is an integral (over p and q) of third-order correlation functions such as δv k ·δB q δn p , where δv k and δB q are the velocity and magnetic-field fluctuations associated with AWs at wave vectors k and q, and δn p is the density fluctuation associated with the slow waves at wave vector p that are driven by the beating of the AWs at wave vectors k and q. For fixed AW amplitudes and fixed B 0 and v A , this driven density fluctuation is proportional to β −1/2 , because as β decreases the thermal pressure is less able to resist the compression along B 0 resulting from the Lorentz force that arises from the beating of the AWs.

Linear Growth of the Parametric Instability
In the linear regime of the parametric instability, the spectrum of AWs propagating in one direction, say E + , is taken to be fixed, and E − ≪ E + . Equation (2.11) then implies that E − increases exponentially in time with growth rate Cohen & Dewar (1974) given the different normalizations of the AW power spectra in the two equations. For example, Equation (2.12) implies that can be compared with the un-numbered but displayed equation under Equation (9) of Cohen & Dewar (1974). As in the present paper, Cohen & Dewar (1974) assumed that slow waves are strongly damped and that the AWs satisfy the randomphase approximation. The present paper builds upon the results of Cohen & Dewar (1974) by investigating the coupled nonlinear evolution of E + and E − . Also, whereas Cohen & Dewar (1974) took the wave vectors to be parallel or anti-parallel to B 0 , the derivation of Equation (2.11) in the present paper allows for obliquely propagating waves.

Conservation of Wave Quanta and Inverse Cascade
To simplify the presentation, it is assumed that (4.1) No generality is lost, because E ± is an even function of k z , and thus it is sufficient to solve for the spectra at positive k z values. Equation (2.11) can be rewritten as the two equations ∂N ∂t is the number of wave quanta per unit k z per unit mass and is the flux of wave quanta in k z -space. Equation (4.2) implies that the number of wave quanta per unit mass, is conserved. The fact that Γ is negative indicates that there is an inverse cascade of wave quanta from large k z to small k z (c.f. Terasawa et al. 1986). The wavenumber drift velocity of the wave quanta, is determined primarily by the smaller of E + and E − .

Exact Solutions to the Wave Kinetic Equations
In this section, several exact solutions to Equation (2.11) are presented under the assumption that k z > 0. The spectra at negative k z follow from the relation E ± (−k z , t) = E ± (k z , t).

Decaying, Balanced Turbulence
One family of exact solutions to Equation (2.11) follows from setting is the Heaviside function. When Equation (5.1) is substituted into Equation (2.11), each side of Equation (2.11) becomes the sum of terms proportional to δ(k z − b) and terms that contain no delta function. By separately equating the two groups of terms, one can show that Equation (5.1) is a solution to Equation (2.11) if . In Appendix A it is shown that Equation (5.4) can be recovered by adding a small amount of nonlinear diffusion to Equation (2.11) and replacing the discontinuous jump in the spectrum at k z = b(t) with a boundary layer. Equation (5.4) implies that, for solutions of the form given in Equation (5.1), the mean-square amplitudes of forward and backward-propagating AWs must be equal just above the break wavenumber b. An exact solution to Equations (5.3) and (5.4) corresponding to decaying turbulence is where a 0 and b 0 are the values of a and b at t = 0. This solution can be further truncated at large k z by setting where q 0 is the value of q at t = 0, which is taken to exceed b 0 . Equations (5.5) through (5.9) can be recovered numerically by solving Equation (2.11) for freely decaying AWs. Whether the spectra satisfy Equations (5.1) and (5.5) through (5.7) or, alternatively, Equations (5.6) through (5.9), the number of wave quanta N tot defined in Equation (4.6) is finite and independent of time.

Forced, Balanced Turbulence
An exact solution to Equations (5.3) and (5.4) corresponding to forced turbulence is where c is a constant and b 0 is the value of b at t = 0. In this solution, the number of wave quanta N tot is not constant, because there is a nonzero influx of wave quanta from infinity. A version of this solution can be realized in a numerical solution of Equation (2.11) by holding E ± fixed at some wavenumber k f , which mimics the effects of energy input from external forcing. In this case, the numerical solution at k z < k f is described by Equations (5.1), (5.10), and (5.11), with b(t) < k f . The solution in Equations (5.10) and (5.11) can be truncated at large k z in a manner analogous to Equation (5.8), but with q = [(πct/2v A ) + (1/q 0 )] −1 , where q 0 is the value of q at t = 0. In this solution, N tot is independent of time. Numerical solutions of Equation (2.11) show, however, that this solution is unstable. If the spectra initially satisfy , then they evolve towards the solution described by Equations (5.5) through (5.9).

Exact Solutions Extending over All k z
In addition to the truncated solutions described in Sections 5.1 and 5.2, Equation (2.11) possesses several exact solutions that extend over all k z . These solutions are unphysical, because they correspond to infinite AW energy and neglect dissipation (which becomes important at sufficiently large k z ) and finite system size (which becomes important at sufficiently small k z ). However, they illustrate several features of the nonlinear evolution of the parametric instability, which are summarized at the end of this section.
The simplest solution to Equation (2.11) spanning all k z is where c ± is a constant. It follows from Equation (4.5) that Equation (5.12) corresponds to a constant flux of AW quanta to smaller k z . In contrast to the truncated E ± ∝ k −1 z forcedturbulence solution in Section 5.2, E + and E − need not be equal in Equation (5.12).
A second, non-truncated, exact solution to Equation (2.11) is given by where a + 0 and a − 0 are the initial values of a + and a − . In this solution, If a + 0 > a − 0 , then E − decays faster than E + , and, after a long time has passed, E − decays to zero while a + decays to the value a + 0 − a − 0 . Conversely, if a − 0 > a + 0 , then E + decays faster than E − , and the turbulence decays to a state in which E + = 0. In the limit that where a 0 = a + 0 = a − 0 . Equations (5.13) and (5.16) are a non-truncated version of the decaying-turbulence solution presented in Section 5.1.
Equations (5.12) and (5.13) can be combined into a more general class of solution, where d + and d − are constants and a ± (t) is given by Equation (5.14). Another type of solution combining k −1 z and k −2 z scalings is where c 0 , c 1 , and c 2 are constants. The exact solutions presented in this section illustrate three properties of the nonlinear evolution of the parametric instability at low β when slow waves are strongly damped. First, when E ± ∝ k −1 z , ∂E ∓ /∂t vanishes. Second, if E ± ∝ k −2 z , then (∂/∂t) ln E ∓ is negative and independent of k z , and E ∓ (k z , t) can be written as the product of a function of k z and a (decreasing) function of time. (More general principles describing the evolution of E ± are summarized in Figure 3 and Equation (6.11).) Third, the parametric instability does not necessarily saturate with E + = E − . For example, in Equations (5.13) and (5.14), when a + 0 = a − 0 , the AWs decay to a maximally aligned state reminiscent of the final state of decaying cross-helical incompressible MHD turbulence (Dobrowolny et al. 1980).

Nonlinear Evolution of the Parametric Instability When Most of the AWs Initially Propagate in the Same Direction
This section describes a numerical solution to Equation (2.11) in which, initially, As in Section 5, k z is taken to be positive, and the spectra at negative k z can be inferred from the fact that E ± (−k z ) = E ± (k z ). The spectra are advanced forward in time using a second-order Runge-Kutta algorithm on a logarithmic wavenumber grid consisting of 2000 grid points. To prevent the growth of numerical instabilities, a nonlinear diffusion term is added to the right-hand side of Equation (2.11), where ν is a constant. The value of ν is chosen as small as possible subject to the constraint that the diffusion term suppress instabilities at the grid scale.
To represent the solution in a way that can be readily compared with spacecraft measurements of solar-wind turbulence, the wavenumber spectra are converted into frequency spectra, where U is the solar-wind velocity, and is the frequency in the spacecraft frame that, according to Taylor's (1938) hypothesis, corresponds to wavenumber k z when the background magnetic field is aligned with the nearly radial solar-wind velocity. The Alfvén speed is taken to be the approximate average of the observed values of v A in three fast-solar-wind streams at r = 0.3 AU (see Table 1 of Marsch et al. (1982) and Table 1a of Marsch & Tu (1990)), v A = 150 km/s. (6.5) In order to compare directly with Figure 2-2c of Tu & Marsch (1995), the solar-wind velocity is taken to be U = 733 km/s. (6.6) The power spectra are initialized to the values and  Tu & Marsch (1995) in the aforementioned fast-solar-wind stream at 10 −5 Hz < f < 10 −4 Hz. The numerical results shown below suggest that the parametric instability has little effect on e + at these frequencies at r = 0.3 AU. The observed f −0.5 scaling in this frequency range is thus presumably inherited directly from the spectrum of AWs launched by the Sun. Like the scaling e + ∝ f −0.5 , the value of σ + is chosen to match the observed spectrum of outward-propagating AWs at 0.3 AU at small f . The reason for the f −2 scaling in e + at large f is that a (parallel) k −2 z spectrum is observed in the solar wind Podesta 2009;Forman et al. 2011) and predicted by the theory of critically balanced MHD turbulence (see, e.g., Goldreich & Sridhar 1995;Mallet et al. 2015). The value of σ − is set equal to a minuscule value (10 −12 σ + ), so that the only source of dynamically important, inward-propagating AWs is the parametric decay of outward-propagating AWs. Figure 2 summarizes the results of the calculation. Between t = 0 and t = 4 hr, e + changes little while e − grows rapidly between roughly 2 and 5 mHz, where the growth rate γ − given in Equation (3.1) peaks. Between t = 4 hr and t = 8 hr, e + develops a broad ∼ 1/f scaling between f = 3 × 10 −4 Hz and f = 3 × 10 −2 Hz, which shuts off the growth of e − at these frequencies. At the same time, e − acquires an ∼ f −2 scaling over much of this same frequency range. Between t = 8 hr and t = 16 hr, the low-frequency limit of the 1/f range of e + decreases to ∼ 10 −4 Hz, and the high-frequency limit of the 1/f range of e + increases to ∼ 0.1 Hz.
The dotted lines in the upper left corner of each panel in Figure 2 show the tracks followed by the values of e + and e − at the low-frequency end of the frequency range in which e + ∝ f −1 in the approximate analytic solution to Equation (2.11) that is described in Appendix B. In this solution, E + and E − are expanded in negative powers of k z at wavenumbers exceeding a time-dependent break wavenumber b(t). Below this wavenumber, E − = 0 and E + = ηk p z , where η and p are constants, and −1 < p < 1. At k z > b, the dominant term in the expansion of E + (E − ) scales like k −1 z (k −2 z ), and the ratios of E + (b + ) to ηb p + and E + (b + ) to E − (b + ) are fixed functions of p, where b + is a wavenumber infinitesimally larger than b. For p = −0.5, E + (b + )/ηb p + = 5/3 and E + (b + )/E − (b + ) = 10, in approximate agreement with the numerical results (see also the right panel of Figure 4). In order to understand the time evolution illustrated in Figure 2, it is instructive to first consider the case in which E ± = c ± k α ± z (6.10) within some interval (k z1 , k z2 ), where c ± and α ± are constants. Equation (2.11) implies that, within this interval, (6.11) If α ∓ > −1, then ln E ± grows at a rate that increases with k z , causing E ± to increase and "harden," in the sense that the best-fit value of α ± within the interval (k z1 , k z2 ) increases. If α ∓ = −1, then E ± does not change. If −2 < α ∓ < −1, then ln E ± decreases at a rate that increases with k z , which causes the best-fit value of α ± within the interval α ∓ -1 -2 0 e ± decreases and flattens e ± decreases and steepens e ± grows and flattens Figure 3. In this figure, it is assumed that the frequency spectra are initially power laws of the form e ± ∝ f α ± , and that α + and α − are both negative. According to Equation (6.11), parametric decay alters both the amplitude and slope of e ± in the manner shown. For example, if e − ∝ f −1.5 , then E + ∝ k −1.5 z , and Equation (6.11) implies that E + decreases at a rate that increases with kz. This in turn implies that e + decreases at a rate that increases with f , so that e + steepens.
(k z1 , k z2 ) to decrease. If α ∓ = −2, then ln E ± decreases at the same rate at all k z , and α ± remains unchanged. Finally, if α ∓ < 2, then E ± decreases at a rate that decreases with k z , which causes the best-fit value of α ± in the interval (k z1 , k z2 ) to increase. These rules are summarized in Figure 3 and apply to e ± ∝ f α ± as well as E ± ∝ k α ± z . Returning to Figure 2, in the early stages of the numerical calculation, e − grows most rapidly at those frequencies at which γ − in Equation (3.1) is largest -namely, the highf end of the f −0.5 range of e + . By the time e − reaches a sufficient amplitude that e + and e − evolve on the same timescale, e − develops a peaked frequency profile extending from some frequency f = f low to some larger frequency f = f high , as illustrated in the upper-right panel of Figure 2. Near f low , de − /df > 0 (i.e., α − > 0), which causes e + to grow. † Near f high , α − < −1, which causes e + to decrease. Thus, e + steepens across the interval (f low , f high ) until it attains a 1/f scaling, at which point e − stops growing between f low and f high . However, at frequencies just below f low , e + and e − both continue to grow, causing f low to decrease. At the same time, e + continues to decrease at larger f where α − < −1. Together, the growth of e + just below f low and the damping of e + at larger f cause the f −1 range of e + to broaden in both directions, i.e., towards both smaller and larger frequencies.
The unique scaling of e − consistent with an e + spectrum ∝ f −1 that is a decreasing function of time is e − ∝ f −2 . Moreover, the scalings e + ∼ f −1 and e − ∼ f −2 are, in a sense, stable, as can be inferred from Figure 3. For example, if α − increases from −2 to a slightly larger value, then e + decreases at a rate that increases with f , causing α + to decrease to a value slightly below −1. This causes e − to decrease at a rate that increases with f , thereby causing α − to decrease back towards −2. A similar "spectral restoring force" arises for any other small perturbation to the values α + = −1 and α − = −2.
It is worth emphasizing, in this context, that the analytic solution presented in Appendix B is approximate rather than exact. As the spectral break frequency decreases past some fixed frequency f 3 , the values of e + and e − at f 3 suddenly jump, but they do not jump to the precise values needed to extend the e + ∼ f −1 and e − ∼ f −2 scalings to smaller f . Instead, the spectra need further "correcting" after the break frequency has swept past in order to maintain the scalings e + ∼ f −1 and e − ∼ f −2 in an approximate way. Also, the decrease in e − that occurs after t = 4 hr is a consequence of the sub- † Although the caption of Figure 3 excludes positive α ∓ to allow use of the words "flattens" and "steepens" in the figure, Equation (6.11) applies for positive α ∓ .

Comparison with Helios Measurements
In the (average) plasma rest frame, the equations of incompressible MHD can be written in the form ∂z ± ∂t + (z ∓ ± v A ) · ∇z ± = −∇Π, (6.12) where z ± = δv ∓ δB/ √ 4πρ are the Elsasser variables, δv and δB are the velocity and magnetic-field fluctuations, ρ is the mass density, v A = B 0 / √ 4πρ is the Alfvén velocity, and Π is the total pressure divided by ρ (Elsasser 1950). Although the solar wind is compressible, Equation (6.12) provides a reasonable approximation for the noncompressive, AW-like component of solar-wind turbulence. As Equation (6.12) shows, the advection velocity of a z ± fluctuation is z ∓ ± v A . This implies, as shown by Maron & Goldreich (2001), that z ± fluctuations propagate along magnetic field lines perturbed by z ∓ . As a consequence, in the solar wind, when the rms magnetic-field fluctuation δB in associated with inward-propagating AWs (z − ) is much smaller than the background magnetic field B 0 , the outward-propagating AWs (z + ) propagate to a good approximation along the direction of B 0 . This is true even if the rms magneticfield fluctuation δB out associated with z + is comparable to B 0 . In the fast solar wind at r < 0.3 AU, the (fractional) cross helicity is high (i.e., E + ≫ E − ), and δB in is indeed small compared to B 0 (Bavassano et al. 2000;Cranmer & van Ballegooijen 2005). Moreover, the background magnetic field at r = 0.3 AU is nearly in the radial direction, because the Parker-spiral magnetic field begins to deviate appreciably from the radial direction only at larger r in the fast wind (Verscharen et al. 2015). Hence, in high-crosshelicity fast-wind streams at r = 0.3 AU, the function e + defined by Equations (6.3) and (6.4) corresponds to a good approximation to the frequency spectrum of outwardpropagating AWs observed by a spacecraft in the solar wind. It is not clear, however, how well e − corresponds to the observed spectrum of inward-propagating AWs, because the inward-propagating AWs follow field lines perturbed by the outward-propagating AWs, which can be inclined relative to the radial direction by a substantial angle. Figure 4 reproduces the t = 4 hr and t = 8 hr panels of Figure 2, but with the same axis ranges as those in Figure 2-2c of Tu & Marsch (1995) to facilitate comparison. Figure 4 also includes a third panel that shows the spectra at t = 32 hr. The e + spectrum in the t = 8 hr panel of Figure 2 shares a number of properties with the e + spectrum in Figure 2-2c of Tu & Marsch (1995), in addition to the f −0.5 scaling at small f that was built in to the numerical calculation as an initial condition. In particular, e + ≫ e − at all frequencies, e + ∼ f −1 at f 3 × 10 −4 Hz, and there is a bump in the e + spectrum at the transition between the f −0.5 and f −1 scaling ranges of e + .
Although this comparison is suggestive, it is not entirely clear how to map time in the numerical calculation to heliocentric distance in the solar wind, because the plasma parameters in the numerical calculation are independent of position and time, whereas they depend strongly upon heliocentric distance in the solar wind. For example, the turbulence is weaker (in the sense of smaller δv 0 /v A ) the closer one gets to the Sun. (See also the discussion following Equation (7.10).) Also, the choice of initial conditions in the numerical calculation artificially prolongs the linear stage of evolution, since in the solar wind there are sources of inward-propagating waves other than parametric instability, such as non-WKB reflection (Heinemann & Olbert 1980;Velli 1993). Nevertheless, as a baseline for comparison, the travel time of an outward-propagating AW from the photosphere to 0.3 AU in the fast-solar-wind model developed by Chandran & Hollweg (2009) is approximately 12 hr.

Discussion of Approximations and Relevance to the Solar Wind
This section critically assesses the assumptions underlying the results in Sections 2 through 6 and the degree to which these assumptions apply to the fast solar wind between r = 10R ⊙ (the approximate perihelion of the Parker Solar Probe) and r = 0.3 AU.

The Weak Turbulence Approximation
A central assumption of the analysis is the weak-turbulence criterion in Equation (2.1). Since E + and E − differ in the solar wind, Equation (2.1) is really two conditions, where ω + nl (ω − nl ) is the inverse of the timescale on which nonlinear interactions modify outward-propagating (inward-propagating) AWs. The contribution to ω ± nl from the parametric instability is The contribution to ω ± nl from one other type of nonlinear interaction is estimated in Section 7.3. The estimate of ∂E ± /∂t in Equation (7.2) follows from Equation (2.11) and setting E ∓ ∼ |k z | α ∓ with α ∓ not very close to −1. A rough upper limit on ω ± nl,PI results from replacing k z E ∓ in Equation (7.2) with (δv ∓ ) 2 , where (δv + ) 2 is the meansquare velocity fluctuation associated with outward-propagating AWs, and (δv − ) 2 is the mean-square velocity fluctuation associated with inward-propagating AWs. This leads to a rough upper limit on ω ± nl,PI because (δv ± ) 2 includes contributions from all wavenumbers and is much larger than the value of k z E ± at some k z . Equation (7.1), with ω nl ∼ ω nl,PI , is thus satisfied provided Bavassano et al. (2000) analyzed Helios measurements of fluctuations in the fast solar wind at r = 0.4 AU and found that (δv − ) 2 ≪ (δv + ) 2 ≃ (60 km/s) 2 . As mentioned above, the typical value of v A in the fast solar wind at r = 0.3 AU is ∼ 150 km/s (Marsch et al. 1982;Marsch & Tu 1990). Near r = 0.3 AU, B 0 ∼ 1/r 2 , ρ ∼ 1/r 2 , and v A ∼ 1/r, and so the typical value of v A in fast-solar-wind streams at r = 0.4 AU is ∼ 112.5 km/s. These measurements indicate that in the fast solar wind at r = 0.4 AU. Since δv ± /v A decreases as r decreases below 0.4 AU (Cranmer & van Ballegooijen 2005;Chandran & Hollweg 2009), the condition ω + nl,PI ≪ |k z |v A is well satisfied at r < 0.4 AU, and the condition ω − nl,PI ≪ |k z |v A is at least marginally satisfied at r < 0.4 AU.
It is worth noting that weak turbulence theory fails when applied to resonant interactions between three AWs, because such interactions occur only when one of the AWs has zero frequency, violating the weak-turbulence ordering (Schekochihin et al. 2012;Meyrand et al. 2015). In contrast, the AW/slow-wave interactions in parametric decay do not involve a zero-frequency mode. Weak turbulence theory is thus in principle a better approximation for the nonlinear evolution of the parametric instability than for incompressible MHD turbulence.

The Low-β Assumption
The assumption that β ≪ 1 is not satisfied at r 0.3 AU ≃ 65R ⊙ , where β is typically ∼ 1, but is reasonable at r 20R ⊙ (Chandran et al. 2011). It is possible that the β ≪ 1 theory presented here applies at least at a qualitative level provided β is simply 1, and indeed this possibility motivates the comparison of the present model with Helios observations. However, further work is needed to investigate how the results of this paper are modified as β increases to values ∼ 1.

Neglect of Other Types of Nonlinear Interactions
Another approximation in Sections 2 through 6 is the neglect of all nonlinear interactions besides parametric decay. One of the neglected interactions is the shearing of inward-propagating AWs by outward-propagating AWs, which makes a contribution to ω − nl that depends on the perpendicular length scale of the AWs. At the perpendicular outer scale L ⊥ (the overall correlation length of the AWs measured perpendicular to B 0 ), the contribution to ω − nl from shearing is approximately is the critical-balance parameter (Goldreich & Sridhar 1995;Ng & Bhattacharjee 1996;Lithwick et al. 2007). Equation (7.5) does not apply when χ is much larger than 1, but direct numerical simulations suggest that χ 1 at r 10R ⊙ for the bulk of the AW energy (J. Perez, private communication). Thus, at r 10R ⊙ , (7.7) As AWs propagate away from the Sun, they follow magnetic field lines, which leads to the approximate scaling L ⊥ ∝ B −1/2 0 . In the WKB limit, the AW frequency in the Sun's frame k z (U + v A ) is independent of r. The scaling k z ∼ 1/(U + v A ) thus serves as a rough approximation for outward-propagating AWs in the turbulent solar wind. At the coronal base (just above the transition region), where k z and L ⊥ have the values k zb and L ⊥b , the value of k zb L ⊥b for the energetically dominant AWs launched by the Sun can be estimated (in essence from the critical-balance condition) as δv and v Ab are the values of δv + and v A at the coronal base (Goldreich & Sridhar 1995;van Ballegooijen & Asgari-Targhi 2016). Together, these scalings lead to the estimate where B 0b and U b are the values of the background magnetic field and solar-wind outflow velocity at the coronal base. Between r = 10R ⊙ and r = 60R where f max is the super-radial expansion factor (Kopp & Holzer 1976). In the fast solar wind within this range of radii, U + v A ≃ 700 − 800 km/s, which is comparable to for the energetically dominant fluctuations launched by the Sun. If we set f max = 9, δv b = 30 km/s, and v Ab = 900 km/s, then Equation (7.9) becomes for the energetically dominant AWs launched by the Sun. Equations (7.7) and (7.10) suggest that it is reasonable to neglect the shearing of inward-propagating AWs by outward-propagating AWs at r 10R ⊙ . On the other hand, at smaller radii, shearing could suppress the growth of inward-propagating AWs that would otherwise result from the parametric instability. Also, the requirement that k z L ⊥ > 1 in order for ω − nl,PI to exceed ω − nl,⊥ could prevent the f −1 range from spreading to frequencies below some minimum (r-dependent) value.
The other nonlinearities in the weak-turbulence wave kinetic equations that are neglected in this paper include interactions involving fast magnetosonic waves, the turbulent mixing of slow waves by AWs, phase mixing of AWs by slow waves, and the shearing of outward-propagating AWs by inward-propagating AWs (Chandran 2008). In-situ measurements indicate that fast waves account for only a small fraction of the energy in compressive fluctuations at 1 AU (Yao et al. 2011;Howes et al. 2012;Klein et al. 2012). Also, fast waves propagating away from the Sun undergo almost complete reflection before they can escape into the corona (Hollweg 1978). These findings suggest that nonlinear interactions involving fast waves have little effect upon the conclusions of this paper. The turbulent mixing of slow waves by AWs acts as an additional slowwave damping mechanism and is thus unlikely to change the conclusions of this paper, which already assume strong slow-wave damping. Phase mixing of AWs by slow waves transports AW energy to larger k ⊥ at a rate that increases with |k z | (Chandran 2008). Although the fractional density fluctuations between r = 10R ⊙ and r = 0.3 AU are fairly small (see, e.g., Tu & Marsch 1995;Hollweg et al. 2010), phase mixing could affect the parallel AW power spectra, and further work is needed to investigate this possibility. The shearing of outward-propagating AWs by inward-propagating AWs is enhanced by non-WKB reflection, which makes this shearing more coherent in time (Velli et al. 1989). The resulting nonlinear timescale for outward-propagating AWs is roughly r/(U + v A ) (Chandran & Hollweg 2009), where U is the solar-wind outflow velocity. This timescale is comparable to the AW propagation time from the Sun to heliocentric distance r, and hence to the parametric-decay timescale at the small-f end of the 1/f range of e + . How this shearing modifies E + (k z ), however, is not clear. For example, shearing by inward-propagating AWs may transport outward-propagating-AW energy to larger k ⊥ = k 2 x + k 2 y at a rate that is independent of |k z |, in which case this shearing would reduce E + (k z ) by approximately the same factor at all k z , leaving the functional form of E + (k z ) unchanged.

Neglect of Spatial Inhomogeneity
In this paper, it is assumed that the background plasma is uniform and stationary. In the solar wind, however, as an AW propagates from the low corona to 0.3 AU, the properties of the ambient plasma seen by the AW change dramatically, with β increasing from ∼ 10 −2 to ∼ 1 and δv rms /v A increasing from ∼ 0.02 to ∼ 0.5 (Bavassano et al. 2000;Cranmer & van Ballegooijen 2005;Chandran et al. 2011). Further work is needed to determine how this spatial inhomogeneity affects the nonlinear evolution of the parametric instability.

Approximate Treatment of Slow-Wave Damping
A key assumption in Sections 2 through 6 is that slow waves are strongly damped, and this damping is implemented by neglecting terms in the wave kinetic equations that are proportional to the slow-wave power spectrum S ± k . There are two sources of error in this approach. First, damping could modify the polarization properties of slow waves, thereby altering the wave kinetic equations. Second, even if S ± k is much smaller than the AW power spectrum A ± k , the neglected parametric-decay terms in the wave kinetic equations for AWs that are proportional to S ± k could still be important, because they contain a factor of β −1 , which is absent in the terms that are retained. This factor arises from the fact that the fractional density fluctuation of a slow wave is ∼ β −1/2 times larger than the fractional magnetic-field fluctuation of an AW with equal energy. These neglected terms act to equalize the 3D AW power spectra A + and A − , and hence to equalize E + and E − . If these neglected terms were in fact important, they could invalidate the solutions presented in Section 6, in which E + ≫ E − . However, in situ observations indicate that E + ≫ E − in the fast solar wind at r = 0.3 AU (Marsch & Tu 1990;Tu & Marsch 1995), which suggests that the neglect of these terms is reasonable. Further work is needed to investigate these issues more carefully.

Conclusion
In this paper, weak turbulence theory is used to investigate the nonlinear evolution of the parametric instability in low-β plasmas. The analysis starts from the wave kinetic equations describing the interactions between AWs and slow waves in weak compressible MHD turbulence. To account for the strong damping of slow waves in collisionless plasmas, terms containing the slow-wave energy density are dropped. The equations allow for all wave-vector directions, but are integrated over the wave-vector components perpendicular to the background magnetic field B 0 (k x and k y ), which leads to equations for the 1D power spectra E + and E − that depend only on the parallel wavenumber k z and time. During parametric decay in a low-β plasma, an AW decays into a slow wave propagating in the same direction and a counter-propagating AW with a frequency slightly smaller than the frequency of the initial AW. The total number of AW quanta is conserved, and the reduction in AW frequencies leads to an inverse cascade of AW quanta towards smaller ω and k z . The energy of each AW quantum is ω, and the decrease in ω during each decay corresponds to a decrease in the AW energy, which is compensated for by an increase in the slow-wave energy. The subsequent damping and dissipation of slow-wave energy results in plasma heating.
The main results of this paper concern the parametric decay of a population of AWs propagating in one direction, say parallel to B 0 , when the counter-propagating AWs start out with much smaller amplitudes. If the initial frequency spectrum e + of the parallelpropagating AWs has a peak frequency f 0 (at which f e + is maximized) and an "infrared" scaling f p at smaller f with −1 < p < 1, then e + acquires a 1/f scaling throughout a range of frequencies that spreads out in both directions from f 0 . At the same time, the anti-parallel-propagating AWs acquire a 1/f 2 spectrum within this same frequency range. If the plasma parameters and infrared e + spectrum are chosen to match conditions in the fast solar wind at a heliocentric distance of 0.3 AU, and the AWs are allowed to evolve for a period of time that is roughly two-thirds of the AW travel time from the Sun to 0.3 AU, the resulting form of e + is similar to the form observed by the Helios spacecraft in the fast solar wind at 0.3 AU. Because the background plasma parameters are time-independent in the analysis of this paper but time-dependent in the plasma rest frame in the solar wind, it is not clear how to map the time variable in the present analysis to heliocentric distance. Nevertheless, the similarity between the spectra found in this paper and the spectra observed by Helios suggests that parametric decay plays an important role in shaping the AW spectra observed in the fast solar wind at 0.3 AU, at least for wave periods 1 hr.
The frequency f * that dominates the AW energy is the maximum of (f e + + f e − ). At the beginning of the numerical calculation presented in Section 6, f * is approximately f 0 = 0.01 Hz. At t = 8 hr in this numerical calculation, f * is the smallest frequency at which e + ∼ f −1 and e − ∼ f −2 , which is ∼ 3 × 10 −4 Hz. This decrease in f * is a consequence of the aforementioned inverse cascade, which transports AW quanta from the initial peak frequency to smaller frequencies. Inverse cascade offers a way to reconcile the observed dominance of AWs at hour-long timescales at 0.3 AU with arguments that the Sun launches most of its AW power at significantly shorter wave periods (Cranmer & van Ballegooijen 2005;van Ballegooijen & Asgari-Targhi 2016).
Further work is needed to relax some of the simplifying assumptions in this paper, including the low-β approximation, the assumption of spatial homogeneity, the simplistic treatment of slow-wave damping, and the neglect of nonlinear interactions other than parametric decay. Further work is also needed to evaluate the relative contributions of parametric decay and other mechanisms to the generation of 1/f spectra in the solar wind. For example, Matthaeus & Goldstein (1986) argued that the f −1 spectrum seen at r = 1 AU at 3 × 10 −6 Hz < f < 8 × 10 −5 Hz is a consequence of forcing at the solar surface, and Velli et al. (1989) argued that the shearing of outward-propagating AWs by the inward-propagating AWs produced by non-WKB reflection causes the outwardpropagating AWs to acquire an f −1 spectrum.
NASA's Parker Solar Probe (PSP) has a planned launch date in the summer of 2018 and will reach heliocentric distances less than 10R ⊙ . The FIELDS (Bale et al. 2016) and SWEAP (Kasper et al. 2015) instrument suites on PSP will provide the first-ever in-situ measurements of the magnetic-field, electric-field, velocity, and density fluctuations in the solar wind at r < 0.29 AU. Although the issues mentioned in the preceding paragraph are sources of uncertainty, the results of this paper lead to the following predictions that will be tested by PSP. First, the 1/f range of e + in fast-solar-wind streams at r < 0.3 AU and f 3 × 10 −4 Hz is produced in situ by parametric decay. As a consequence, the 1/f range of e + in fast-solar-wind streams will be much more narrow at small r than at r = 0.3 AU. As AWs propagate away from the Sun, the frequency range (f min , f max ) in which e + ∼ 1/f spreads out in both directions from the near-Sun peak frequency (the maximum of f e + ). Thus, f min will be larger closer to the Sun, and f max will be smaller. Finally, during epochs in which the local magnetic field is aligned with the relative velocity between the plasma and the spacecraft (see the discussion in Section 6.2), the spectrum of e − will scale like 1/f 2 in the frequency interval (f min , f max ). I thank Phil Isenberg for discussions about the work of Cohen & Dewar (1974) and the three anonymous reviewers for helpful comments that led to improvements in the manuscript. This work was supported in part by NASA grants NNX15AI80, NNX16AG81G, and NNX17AI18G, NASA grant NNN06AA01C to the Parker Solar Probe FIELDS Experiment, and NSF grant PHY-1500041.
where H(x) is the Heaviside function, a n , b, and c n are functions of time, and η and p are constants, with − 1 < p < 1.

(B 7)
Equations (B 5) through (B 7) imply that the break wavenumber b(t) decreases in time, and that the ratios E + (b + )/ηb p + and E + (b + )/E − (b + ) remain constant, where b + is a wavenumber infinitesimally larger than b. For example, if the infrared spectral index p is −1/2, then