Three-dimensional simulation of the fast solar wind driven by compressible magnetohydrodynamic turbulence

TL;DR

This study uses 3D compressible MHD simulations to demonstrate that the fast solar wind is heated and accelerated by turbulence driven by parametric decay instability and Alfven wave reflection, aligning with the wave/turbulence driven scenario.

Contribution

It provides a self-consistent 3D simulation showing how compressible MHD turbulence, driven by parametric decay instability, accelerates and heats the fast solar wind.

Findings

Turbulence is imbalanced and anisotropic.

Density fluctuations originate from parametric decay instability.

Alfven wave reflection triggers turbulence leading to wind acceleration.

Abstract

Using a three-dimensional compressible magnetohydrodynamic (MHD) simulation, we have reproduced the fast solar wind in a direct and self-consistent manner, based on the wave/turbulence driven scenario. As a natural consequence of Alfvenic perturbations at its base, highly compressional and turbulent fluctuations are generated, leading to heating and acceleration of the solar wind. The analysis of power spectra and structure functions reveals that the turbulence is characterized by its imbalanced (in the sense of outward Alfvenic fluctuations) and anisotropic nature. The density fluctuation originates from the parametric decay instability of outwardly propagating Alfven waves and plays a significant role in the Alfven wave reflection that triggers turbulence. Our conclusion is that the fast solar wind is heated and accelerated by compressible MHD turbulence driven by parametric decay instability and resultant Alfven wave reflection.

Equations

$$\frac{\partial }{\partial t}\rho$$

$$\displaystyle\frac{\partial}{\partial t}\left(\rho\mbox{\boldmath$v$}\right)$$

$$\displaystyle\frac{\partial}{\partial t}\mbox{\boldmath$B$}$$

$$\frac{\partial }{\partial t}e$$

$$\displaystyle e=\frac{p}{\gamma-1}+\frac{1}{2}\rho\mbox{\boldmath$v$}^{2}+\frac{\mbox{\boldmath$B$}^{2}}{8\pi},\ \ \ p_{T}=p+\frac{\mbox{\boldmath$B$}^{2}}{8\pi}.$$

$$\displaystyle\nabla\approx\mbox{\boldmath$e$}_{r}\frac{\partial}{\partial r}+\mbox{\boldmath$e$}_{\theta}\frac{1}{r}\frac{\partial}{\partial\theta}+\mbox{\boldmath$e$}_{\phi}\frac{1}{r}\frac{\partial}{\partial\phi},$$

$$\displaystyle\mbox{\boldmath$g$}=-\frac{GM_{\odot}}{r^{2}}\mbox{\boldmath$e$}_{r}.$$

$$\displaystyle\mbox{\boldmath$q$}_{\rm cnd}\cdot\mbox{\boldmath$e$}_{r}$$

$$\displaystyle\mbox{\boldmath$q$}_{\rm cnd}\cdot\mbox{\boldmath$e$}_{\theta}$$

$$\xi (r)=\min \left(1,\frac{r_{\mathrm{sw}}^2}{r^2}\right),$$

$$-1.44\times 10^{-2}\mathrm{rad}\le \theta ,\varphi \le 1.44\times 10^{-2}\mathrm{rad}.$$

$$\displaystyle\mbox{\boldmath$z$}_{\perp}^{\pm}=z_{\theta}^{\pm}\mbox{\boldmath$e$}_{\theta}+z_{\phi}^{\pm}\mbox{\boldmath$e$}_{\phi},\ \ \ \ z_{\theta,\phi}^{\pm}=v_{\theta,\phi}\mp\frac{B_{\theta,\phi}}{\sqrt{4\pi\rho}}.$$

$$\displaystyle\int E^{\pm}\left(k_{\perp}\right)dk_{\perp}=\frac{1}{L_{\perp}^{2}}\int r^{2}d\theta d\phi\ \mbox{\boldmath$z$}_{\perp}^{2}\left(\theta,\phi\right),$$

$$\displaystyle\frac{\partial}{\partial t}\mbox{\boldmath$z$}^{\pm}+\left(\mbox{\boldmath$z$}^{\mp}\cdot\nabla\right)\mbox{\boldmath$z$}^{\pm}=0,\ \ \mbox{\boldmath$z$}^{\pm}=\mbox{\boldmath$v$}\mp\mbox{\boldmath$B$}/\sqrt{4\pi\rho}.$$

$$E^{\pm }\propto k_{\perp }^{-2\pm \alpha },$$

$$E^{\pm }\propto k_{\perp }^{-3/2\pm \tilde{\alpha }},$$

$$S{F}_2^{\pm }\left(l_{\parallel },l_{\perp }\right)$$

$${\gamma }_{\mathrm{PDI}}={\tilde{\gamma }}_{\mathrm{GD}}2\pi f_0-{\gamma }_{\mathrm{acc}}-{\gamma }_{\exp },$$

$${\omega }_{\mathrm{ref}}=\left(v_r+v_{A,r}\right)\frac{\partial }{\partial r}\ln \left(r{\rho }^{1/4}\right),$$

Full text

magnetohydrodynamics (MHD) – methods: numerical – Sun: solar wind – turbulence

1 Introduction

One of the most important unsolved problems in astrophysics is the driving mechanism of the solar wind. In addition to the close relation to the coronal heating problem (Parker, 1958), understanding solar wind acceleration is required for stellar rotational evolution (e.g. Brun & Browning, 2017) and for space weather forecasting at Earth and at exoplanets (e.g. Garraffo et al., 2016). It is now widely accepted that the ultimate energy source of the solar wind comes from the surface convection. Only $0.1-1\%$ of the photospheric energy flux is sufficient to drive the solar wind (Withbroe & Noyes, 1977). Indeed, some observations confirm the sufficient upward energy transport (De Pontieu et al., 2007; McIntosh et al., 2011), whereas we should note that these observations are still controversial (see e.g. Thurgood et al., 2014). An unsolved issue regarding solar wind acceleration is the thermalization process, specifically the nature of solar wind turbulence where the solar wind is accelerated. In-situ observations near $1{\rm\ au}$ indicate that the turbulent dissipation (cascading) accounts for ongoing heating of the solar wind plasma against adiabatic expansion (Carbone et al., 2009). However, since the plasma condition is very different between the Earth’s orbit ($r/R_{\odot}\approx 200$ where $R_{\odot}$ denotes the solar radius) and the wind acceleration region ($r/R_{\odot}\approx 10$), it is risky to simply assume the same situation. In fact, several observations show that the density fluctuation is large near the Sun ($r/R_{\odot}\approx 1-10$, see Miyamoto et al., 2014; Hahn et al., 2018).

In this study, based on the wave/turbulence-driven (WTD) scenario , we perform three-dimensional, compressible MHD simulation of the fast solar wind. The simulation is conducted in a self-consistent manner; direct calculation of MHD equations enables us to consider the evolutions of the mean field and fluctuation simultaneously. The compressibility is critical for two reasons. First, due to the compression of plasma, formation of shock waves is allowed, which can contribute to the heating of the solar wind . Second, the parametric decay instability (PDI) is incorporated. PDI is an instability of Alfvén wave and can grow in the wind acceleration region (Suzuki & Inutsuka, 2006; Tenerani & Velli, 2013; Shoda et al., 2018b; Réville et al., 2018; Chandran, 2018) and activate the turbulence by introducing various energy cascading channels (e.g. Shoda & Yokoyama, 2018). In fact, reduced MHD model cannot account for the solar wind heating without density fluctuation that is likely to be generated by PDI. In addition to the compressibility, three-dimensionality is crucial for solving turbulence. In general, lower-dimensional (1D or 2D) simulations show different behavior compared with 3D ones (e.g. Shoda & Yokoyama, 2018). Therefore, both compressibility and three-dimensionality appear to be crucial for the study of solar wind turbulence.

2 Method

We simulate the fast solar wind from the polar region in the solar minimum. The basic equations are ideal MHD equations with gravity and thermal conduction in the spherical coordinates:

[TABLE]

where

[TABLE]

\hat{\mbox{\boldmathI}} stands for the unit tensor, $g$ and \mbox{\boldmathq}_{\rm cnd} are the gravitational acceleration and thermal conductive flux, respectively. We employ the adiabatic specific heat ratio of monatomic gas $\gamma=5/3$. In solving these equations, the spherical coordinate system is used on the $\theta=\pi/2$ plane that symmetrizes $\theta$ and $\phi$ with respect to $\nabla$ operator as

[TABLE]

where \mbox{\boldmathe}_{r,\theta,\phi} denotes the unit vector in each direction. Due to the small horizontal ($\theta$ and $\phi$) extension of our numerical domain, this approximation yields at most $0.1\%$ error compared with the usual spherical coordinate. Note that we replace $r\sin\theta$ with $r$ in Eq. (6) and the deviation between the two is in the order of $(\theta-\pi/2)^{2}$ near $\theta=\pi/2$, which yields $10^{-3}$ in our setting. Since we are simulating polar wind, $g$ is given as

[TABLE]

The radiative cooling is not considered because we do not solve the atmosphere below the transition region where radiation plays a role. The thermal conduction instead dominates the energy balance in the corona and solar wind. We employ a conductive flux that mimics the Spitzer-Härm type one (Spitzer & Härm, 1953) as

[TABLE]

where $\kappa_{0}=10^{-6}$ in cgs unit, $\xi(r)$ is a quenching due to the large mean free path of electron (Hollweg, 1974) and \hat{b}_{r,\theta,\phi}=B_{r,\theta,\phi}/\left|\mbox{\boldmathB}\right|. $\phi$ component of \mbox{\boldmathq}_{\rm cnd} is given in a similar way as $\theta$ component. The quenching term $\xi(r)$ is given as a function of $r$:

[TABLE]

where $r_{\rm sw}=5R_{\odot}$ in this study.

An additional quenching $f_{q}=0.1$ is used to avoid the severe restriction of time step. Numerical results do not depend on the value of $f_{q}$ because the conductions in $\theta$ and $\phi$ directions are sufficiently fast to homogenize the temperature on the horizontal plane. Although the conductive flux in our model is not the same as Spitzer-Härm type one, the most important effect of thermal conduction, that is cooling by the radial heat transport, is appropriately solved. Thus, our simplified conductive flux is appropriate for the solar wind simulation.

The numerical domain extends from the coronal base ($r=1.02R_{\odot}$) to $r=20R_{\odot}$ with horizontal size $20{\rm\ Mm}$ at the bottom, which yields the range of $\theta$ and $\phi$ as

[TABLE]

An additional numerical domain with coarser grids is prepared beyond the top boundary up to $r\approx 10^{3}R_{\odot}$, far enough to ensure that no fluctuations reach within the time of simulation. This method is validated because no physical quantities can propagate back, against the super-Alfvénic solar wind, to the numerical domain in the quasi-steady state. As an initial condition, we impose the isothermal Parker wind with temperature $T=1.1\times 10^{6}{\rm\ K}$ with radially extending magnetic field embedded. The bottom boundary is as follows. The density, temperature and radial magnetic field are fixed to $\rho=8.5\times 10^{-16}{\rm\ g\ cm^{-3}}$, $T=6\times 10^{5}{\rm\ K}$ and $B_{r}=2{\rm\ G}$, respectively, without any variation in $\theta$ and $\phi$ directions.

If we consider the expansion factor of magnetic field, $B_{r}=2{\rm\ G}$ is consistent with the source of the fast solar wind (Fujiki et al., 2015),

The radial derivative of $v_{r}$ is fixed to zero, which allows the supply of mass into the numerical domain. The inward Elsässer variables $z_{\theta,\phi}^{-}$ are set to be transmissive at the bottom. Here the upward ($z_{\perp}^{+}$) and downward ($z_{\perp}^{-}$) Elsässer variables are defined as

[TABLE]

The amplitude of the upward Elsässer variable at the bottom boundary is $64{\rm\ km\ s^{-1}}$ corresponding to the observed non-thermal velocity of $32{\rm\ km\ s^{-1}}$ (e.g. Banerjee et al., 2009; Landi & Cranmer, 2009).

The typical horizontal length scale of the upward Elsässer variable is fixed to the horizontal size of the simulation domain. The basic equations are numerically integrated by the combination of 3rd-order SSP Runge–Kutta method (Shu & Osher, 1988) and HLLD Riemann solver (Miyoshi & Kusano, 2005) with spatial reconstruction a combination of 2nd-order MUSCL (van Leer, 1979) and 5th-order MP5 (Suresh & Huynh, 1997) methods. The number of grid points is $(6600,192,192)$ in $(r,\theta,\phi)$ directions, respectively. The super-time-stepping method is used to solve the thermal conduction (Meyer et al., 2014). To remove the numerically generated finite \nabla\cdot\mbox{\boldmathB}, we employ the hyperbolic cleaning method (Dedner et al., 2002).

3 Results & Discussion

Figure 1 shows the snapshots of temperature (Panel a), radial velocity (Panel b), anti-Sunward Elsässer variable (Panel c) and Sunward Elsässer variable (Panel d), respectively, on the meridional plane after the quasi-steady state is achieved. The maximum temperature exceeds $10^{6}{\rm\ K}$ and the termination radial velocity approximates $600{\rm\ km\ s^{-1}}$, ensuring the successful reproduction of the fast solar wind. As a natural consequence of fast thermal conduction, no fine structuring is observed in the temperature map. The panel of $v_{r}$ shows that, in addition to the gradual acceleration of the solar wind, ubiquitous local (or discontinuous) enhancements are observed. According to the previous 1D simulations (Suzuki & Inutsuka, 2005; Shoda et al., 2018a), these fluctuations are large-amplitude slow mode waves that can at least partially contribute to the heating of the solar wind.

Panels c and d in Figure 1 show the evolution of waves and turbulence. Note that $z_{\phi}^{+}$ and $z_{\phi}^{-}$ correspond to anti-Sunward and Sunward Alfvén wave characteristics in the linear regime. $z_{\phi}^{+}$ maintains the coherent structure in the entire simulation domain while $z_{\phi}^{-}$ shows an evidence of strong turbulent distortion.

It is evident, especially in $4\lesssim r/R_{\odot}\lesssim 10$, that $z_{\phi}^{-}$ has much finer transverse structure than $z_{\phi}^{+}$. This feature is more quantitatively observed in the Elsässer power spectra with respect to perpendicular wave number defined as

[TABLE]

where $k_{\perp}$ is the wave number perpendicular to the mean field direction ($r$ axis) and $L_{\perp}$ denotes the horizontal extension of the simulation domain at $r$. Note that $E^{\pm}\left(k_{\perp}\right)$ reflects the spatial structures of anti-Sunward and Sunward Elsässer variables perpendicular to the mean field. Solid lines in Figure 2 show $E^{\pm}\left(k_{\perp}\right)$ calculated at $r/R_{\odot}=6.46$ (Panel a) and $r/R_{\odot}=17.3$ (Panel b). Also shown by dashed lines correspond to the power-law fitting in the inertial range ($2\leq k_{\perp}L_{\perp}/(2\pi)\leq 20$). We observe different inertial-range power indices between $E^{+}$ and $E^{-}$ for both locations. Specifically, the anti-Sunward component has flatter (harder) power spectrum than that of Sunward one. This spectral behavior is consistent with structure difference between $z_{\phi}^{+}$ and $z_{\phi}^{-}$ observed in Figure 1, because the flatter power spectrum is associated with finer structures.

Another interesting feature is that, as $r$ gets larger, the power indices of both $E^{+}$ and $E^{-}$ approach the Kolmogorov’s index $-5/3$, which is observed in the magnetic power spectrum in the solar wind (e.g. Bruno & Carbone, 2013).

A brief description of the spectral difference is given as follows. In the regime of reduced MHD, neglecting the inhomogeneity of the background, the evolution of Alfvén waves is described as follows (e.g. Priest, 2014):

[TABLE]

Thus, the nonlinear wave-wave interaction is invoked by the collision of counter-propagating waves. Note, however, that, in the presence of background inhomogeneity, this is not the case because Elsässer variables are no longer pure characteristics of Alfvén waves (anomalous components, see Velli et al., 1989; Velli, 1993; Perez & Chandran, 2013). Eq. (14) shows that the energy cascading timescale of $z^{\pm}$ is determined by $z^{\mp}$. Specifically when $z^{+}>z^{-}$, the cascading of $z^{-}$ proceeds faster than $z^{+}$, leading to structure difference . In terms of relaxation process, this process is called dynamical alignment (e.g. Biskamp, 2003), in which the minor component of $z^{\pm}$ decays faster than the major one.

More quantitative explanation of the spectral imbalance is also given both numerically and analytically. The theory and simulation of the incompressible MHD turbulence show (Boldyrev & Perez, 2009)

[TABLE]

while the strong turbulence (EDQNM) theory predicts (Grappin et al., 1983)

[TABLE]

where $\alpha$ and $\tilde{\alpha}$ depend on the degree of imbalance $E^{+}/E^{-}$. Compressible MHD turbulence also possibly exhibits the similar spectral difference (Perez et al., 2012). Though not perfectly, our results are at least qualitatively consistent with these predictions. The summation of power indices shifts from $-3.5$ ($r/R_{\odot}=6.46$) to $-3.2$ ($r/R_{\odot}=17.3$), suggesting the weak-to-strong transition of turbulence.

The anisotropy is another factor that characterizes the structure of turbulence. In the presence of mean magnetic field, the structure of turbulence is expected to be anisotropic (e.g. Goldreich & Sridhar, 1995). To see the degree of anisotropy, we often use the 2nd-order structure function (Cho & Lazarian, 2003; Verdini et al., 2015) defined as

[TABLE]

where the bracket denotes the averaging operator over $\theta$, $\phi$ and $t$. Here the structure function is defined based on two assumptions; the turbulence is anisotropic with respect to $r$ axis and isotropic in the $\theta\phi$ plane. The first assumption is justified since the mean magnetic field is perfectly aligned with $r$ axis. In some simulations of the solar wind, the second assumption is violated because the wind expansion introduces another form of anisotropy (Dong et al., 2014). The flow direction of the solar wind is therefore an additional anisotropy axis that forms 3D anisotropy of the wind turbulence (Verdini et al., 2018). In our calculation, because we simulate the wind from the polar region, the flow direction is aligned with the mean field. The definition of the structure function is justified.

Figure 3 shows $SF_{2}^{\pm}$ measured at $r/R_{\odot}=10$. Both $SF_{2}^{+}$ and $SF_{2}^{-}$ rapidly increase in $l_{\perp}$ direction, showing that field-aligned structures are generated preferentially in the solar wind. This is consistent with previous works (e.g. Cho & Lazarian, 2003; Shoda & Yokoyama, 2018). A difference of anisotropy is also observed; the minor component ($SF_{2}^{-}$) shows larger degree of anisotropy than the major component ($SF_{2}^{+}$). This is consistent with the result of phenomenological study of Alfvén wave turbulence (Beresnyak & Lazarian, 2008). We also note that the fluctuations in the solar wind show a similar behavior (Wicks et al., 2011).

The results and discussion above show that the structure of solar wind turbulence is consistent with that of incompressible MHD (Alfvén-wave) turbulence. This does not mean that the compressibility is ignorable in the solar wind. In fact, the density fluctuation appears to be required to account for the heating rate in the solar wind (van Ballegooijen & Asgari-Targhi, 2016). Specifically, the density fluctuations act as reflectors of anti-Sunward Alfvén waves, playing an indirect but critical role in the onset of turbulence. Figure 4 shows the distribution, origin and role of density fluctuation. Panel (a) shows the turbulent and discontinuous structure of fractional density fluctuation $\delta\rho/\rho$ on the meridional plane. Note that the magnitude of density fluctuation is, at least locally, as large as the mean density. Panel (b) displays the magnitude of density fluctuation, and the growth rate of parametric decay instability (PDI) of Alfvén waves, $\gamma_{\rm PDI}$, in unit of ${\rm\ mHz}$. In the accelerating, expanding solar wind, $\gamma_{\rm PDI}$ is given as (Tenerani & Velli, 2013; Shoda et al., 2018b)

[TABLE]

where $\tilde{\gamma}_{\rm GD}$ is the normalized growth rate (Goldstein, 1978; Derby, 1978) in the homogeneous system, $f_{0}=10^{-3}{\rm\ Hz}$ is the lowest frequency of the injected Alfvén waves, and $\gamma_{\rm acc}$ and $\gamma_{\rm exp}$ are the suppression of density fluctuation by the wind acceleration and expansion (see Shoda et al., 2018b). Note that we show $\gamma_{\rm PDI}/2$ rather than $\gamma_{\rm PDI}$ for better visualization. A clear spatial correlation between $\gamma_{\rm PDI}$ and $\delta\rho/\rho$ shows that the density fluctuation originates in PDI. Thus the density fluctuation comes from the PDI. Panel (c) shows how the density fluctuation affects the Alfvén wave propagation by displaying the reflection rate of anti-Sunward Alfvén waves given as (Heinemann & Olbert, 1980)

[TABLE]

where $v_{A,r}=B_{r}/\sqrt{4\pi\rho}$. To see the role of density fluctuation, we compare $\omega_{\rm ref}$ calculated in two ways. First, we calculate $\omega_{\rm ref}$ at each given time and space and average it over time and $\theta\phi$ plane. Second, we first average $v_{r}$, $v_{A,r}$ and $\rho$ in time and horizontal space, and calculate $\omega_{\rm ref}$ using Eq. (19). The former and latter correspond to the reflection rates with and without density fluctuations, because any fluctuations are smoothed out by averaging the background values ($v_{r}$, $v_{A,r}$, $\rho$). The comparison of these two values are shown in Panel (c). The density fluctuation enhances the reflection rate by a factor of 10 or larger, and therefore the density fluctuation plays a dominant role in the Alfvén wave reflection that triggers turbulence. To summarize, we have shown that the density fluctuation, excited by the PDI, plays a crucial role in the turbulence trigger, and therefore, the compressibility is far from negligible in the solar wind acceleration.

4 Conclusion

Our three-dimensional MHD simulation reproduces the fast solar wind as a natural consequence of Alfvén-wave injection from the coronal base, thus supporting the wave/turbulence-driven model of the fast solar wind. The turbulence is characterized by imbalance (Figure 2), anisotropy (Figure 3), and compressibility (Figure 4). The structure of turbulence is well described by the incompressible or reduced MHD turbulence in which the compressibility is ignored. To discuss the turbulent dissipation and heating, however, compressibility plays a crucial role because the wave reflection, the source of Alfvén wave turbulence, is driven dominantly by density fluctuations excited by the parametric decay instability of Alfvén waves.

We acknowledge the anonymous referee for a number of helpful and constructive comments. We would also like to thank Marco Velli and Takuma Matsumoto for many useful discussions and comments. Numerical computations were carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. M. Shoda is supported by Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows and by the NINS program for cross-disciplinary study (Grant Nos. 01321802 and 01311904) on Turbulence, Transport, and Heating Dynamics in Laboratory and Solar/Astrophysical Plasmas: ”SoLaBo-X”. T.K. Suzuki is supported in part by Grants-in-Aid for Scientific Research from the MEXT of Japan, 17H01105. M. Asgari-Targhi was supported under contract NNM07AB07C from NASA to the Smithsonian Astrophysical Observatory (SAO) and contract SP02H1701R from Lockheed Martin Space and Astrophysics Laboratory (LMSAL) to SAO. T. Yokoyama is supported by JSPS KAKENHI Grant Number 15H03640.