Energy Transfer by Nonlinear Alfvén Waves in the Solar Chromosphere and Its Effect on Spicule Dynamics, Coronal Heating, and Solar Wind Acceleration

We used 1D magnetohydrodynamic equations based on the axial symmetry assumption of the magnetic flux tube. The surface of the axisymmetric flux tube is defined by the poloidal and toroidal axes, which are denoted in this study by x and ϕ. The basic equations in cgs units are written as follows:

The mass conservation law is presented by

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{1}{A}\displaystyle \frac{\partial }{\partial x}(\rho {v}_{x}A)=0,\end{eqnarray} \tag{ 1 }$$

where ρ, v_x, and A are the mass density, poloidal component of velocity, and cross section of the flux tube, respectively.

The energy conservation law is presented by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{p}{\gamma -1}+\displaystyle \frac{1}{2}\rho {v}^{2}+\displaystyle \frac{{B}^{2}}{8\pi }\right)\\ \ \ +\displaystyle \frac{1}{A}\displaystyle \frac{\partial }{\partial x}\left[A\left\{\left(\displaystyle \frac{\gamma p}{\gamma -1}+\displaystyle \frac{\rho {v}^{2}}{2}+\displaystyle \frac{{B}_{\phi }^{2}}{4\pi }\right){v}_{x}-\displaystyle \frac{{B}_{x}}{4\pi }({B}_{\phi }{v}_{\phi })\right\}\right]\\ \ \ \ \ \ =\rho {v}_{x}\displaystyle \frac{\partial }{\partial x}\left(\displaystyle \frac{{{GM}}_{\odot }}{r}\right)-\displaystyle \frac{1}{A}\displaystyle \frac{\partial }{\partial x}({{AF}}_{c})-{Q}_{{\rm{rad}}},\end{array}\end{eqnarray} \tag{ 2 }$$

where p, B_x, B_ϕ, v_ϕ, and γ are the gas pressure, poloidal and toroidal components of the magnetic field, toroidal component of velocity, and the specific heat ratio, which is set to 5/3, respectively. ${v}^{2}={v}_{x}^{2}+{v}_{\phi }^{2}$ and ${B}^{2}={B}_{x}^{2}+{B}_{\phi }^{2}$. G and ${M}_{\odot }$ are the gravitational constant and the solar mass. r is the distance from the Sun center. F_c and ${Q}_{\mathrm{rad}}$ represent the heat conduction flux and radiative cooling term, respectively, as described in Section 2.3.

The poloidal component of the equation of motion is presented by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial (\rho {v}_{x})}{\partial t}+\displaystyle \frac{\partial p}{\partial x}+\displaystyle \frac{1}{A}\displaystyle \frac{\partial }{\partial x}\left\{\left(\rho {v}_{x}^{2}+\displaystyle \frac{{B}_{\phi }^{2}}{8\pi }\right)A\right\}\\ \qquad -\rho {v}_{\phi }^{2}\displaystyle \frac{\partial \mathrm{ln}\sqrt{A}}{\partial x}-\rho \displaystyle \frac{\partial }{\partial x}\left(\displaystyle \frac{{{GM}}_{\odot }}{r}\right)=0.\end{array}\end{eqnarray} \tag{ 3 }$$

The toroidal component of the equation of motion is presented by

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {v}_{\phi })}{\partial t}+\displaystyle \frac{1}{A\sqrt{A}}\displaystyle \frac{\partial }{\partial x}\left\{A\sqrt{A}\left(\rho {v}_{x}{v}_{\phi }-\displaystyle \frac{{B}_{x}{B}_{\phi }}{4\pi }\right)\right\}=0.\end{eqnarray} \tag{ 4 }$$

The toroidal component of the induction equation is presented by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{\phi }}{\partial t}+\displaystyle \frac{1}{\sqrt{A}}\displaystyle \frac{\partial }{\partial x}\left(\sqrt{A}({v}_{x}{B}_{\phi }-{v}_{\phi }{B}_{x})\right)=0.\end{eqnarray} \tag{ 5 }$$

The poloidal magnetic flux conservation is presented by

$$\begin{eqnarray}&&{B}_{x}A={\rm{const.}}\end{eqnarray} \tag{ 6 }$$

Finally, we note that the poloidal axis x is not always parallel to the radial axis r. They are related to each other as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{dx}}{{dr}}=\sqrt{1+{\left(\displaystyle \frac{d\sqrt{A}}{{dr}}\right)}^{2}}.\end{eqnarray} \tag{ 7 }$$

The assumption of the background magnetic field is described here in detail. The cross section of the flux tube A is related to r through the filling factor f as $A(r)=4\pi {r}^{2}f(r)$. f determines the geometry of the flux tube. We consider the axisymmetric magnetic flux tube from the photosphere to interplanetary space. The outer boundary of our simulation is set to 0.5 au. In the lower atmosphere, the magnetic flux tube expands exponentially such that the magnetic pressure inside the flux tube balances out with the ambient plasma gas pressure, which decreases with the scale height ${H}_{\mathrm{ph}}={R}_{g}{T}_{\mathrm{eff}}/({\mu }_{\mathrm{ph}}{g}_{\odot })$. Here, ${R}_{g}=8.31\times {10}^{7}$ erg K⁻¹ mol⁻¹ is the gas constant, ${T}_{\mathrm{eff}}=5770$ K, ${\mu }_{\mathrm{ph}}=1.3$ is the mean molecular weight on the photosphere, and ${g}_{\odot }$ is the gravitational acceleration on the solar surface. The filling factor f in this layer is expected to be ${f}_{\mathrm{atm}}(r)={f}_{\mathrm{ph}}\exp \{{r}_{\odot }/(2{H}_{\mathrm{ph}})(1-{r}_{\odot }/r)\}$, where ${f}_{\mathrm{ph}}$ is the coverage of the open magnetic flux tube on the photosphere. By using ${f}_{\mathrm{atm}}$, ${B}_{x}={r}_{\odot }^{2}{f}_{\mathrm{ph}}{B}_{\mathrm{ph}}/({r}^{2}{f}_{\mathrm{atm}})$ satisfies the condition that ${B}_{x}^{2}/(8\pi )={p}_{\mathrm{atm}}$, where ${p}_{\mathrm{atm}}$ is the solution of the hydrostatic equilibrium. In the lower atmosphere, where $r={r}_{\odot }+h$ ($h\ll {r}_{\odot }$), we obtain ${f}_{\mathrm{atm}}(h)={f}_{\mathrm{ph}}{e}^{h/(2{H}_{\mathrm{ph}})}$. This exponential expansion of flux tube is assumed to stop at some height where it merges with the neighboring flux tube. Above this height (i.e., the merging height H_m), the magnetic pressure dominates the gas pressure and the flux tube extends vertically. The poloidal magnetic field strength in this layer is assumed to be almost constant around $\overline{B}={B}_{\mathrm{ph}}{f}_{\mathrm{ph}}/{f}_{\mathrm{atm}}({H}_{m})={B}_{\mathrm{ph}}{e}^{-{H}_{m}/(2{H}_{\mathrm{ph}})}$ through the upper chromosphere and coronal base. Thus, $\overline{B}={B}_{\mathrm{ph}}{e}^{-{H}_{m}/(2{H}_{\mathrm{ph}})}$ roughly represents the area-averaged magnetic field strength in the coronal hole from which the solar wind emanates. It should be noted that various flux tube models in the lower atmosphere have been considered, for example, by Hasan et al. (2003, 2005) and Cranmer & van Ballegooijen (2005). The different magnetic field geometries would lead to different results. Their significance should be tested in future studies as long as we rely on the 1D simulation.

The flux tube is assumed to expand superradially again in the extended corona such that the interplanetary space is filled with the open flux tube. We characterize this expansion with the coronal loop height H_l. The functional form of the filling factor in this layer ${f}_{\mathrm{wind}}(r)$ is suggested by Kopp & Holzer (1976). Based on these considerations, the profile of the filling factor f(r) is determined as follows:

$$\begin{eqnarray}&&{f}_{\mathrm{atm}}(r)={f}_{m}\tanh \left[\displaystyle \frac{{f}_{\mathrm{ph}}}{{f}_{m}}\exp \left\{\displaystyle \frac{{r}_{\odot }}{2{H}_{\mathrm{ph}}}\left(1-\displaystyle \frac{{r}_{\odot }}{r}\right)\right\}\right],\end{eqnarray} \tag{ 8 }$$

where ${f}_{m}={f}_{\mathrm{ph}}{e}^{{H}_{m}/(2{H}_{\mathrm{ph}})}$

$$\begin{eqnarray}&&{f}_{\mathrm{wind}}(r)=\displaystyle \frac{{e}^{(r-{r}_{\odot }-{H}_{l})/{\sigma }_{l}}+{f}_{m}-(1-{f}_{m}){e}^{-({H}_{l}/{\sigma }_{l})}}{{e}^{(r-{r}_{\odot }-{H}_{l})/{\sigma }_{l}}+1}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\begin{array}{rcl}\hat{f}(r) & = & {f}_{\mathrm{atm}}(r)+\displaystyle \frac{1}{2}\left(\max [{f}_{\mathrm{wind}}(r),{f}_{m}]-{f}_{\mathrm{atm}}(r)\right)\\ & & \times \,\left\{1+\tanh \left(\displaystyle \frac{r-{r}_{\odot }-{H}_{l}}{{H}_{l}}\right)\right\}\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&f(r)={f}_{\mathrm{ph}}+(1-{f}_{\mathrm{ph}})\displaystyle \frac{\hat{f}(r)-\hat{f}({r}_{\odot })}{1-\hat{f}({r}_{\odot })}.\end{eqnarray} \tag{ 11 }$$

The key parameters of f(r) are ${f}_{\mathrm{ph}}$, H_m, and H_l. ${\sigma }_{l}$ in Equation (9) is set to H_l. The manner by which the properties of solar and stellar wind depend on ${f}_{\mathrm{ph}}$ has already been well investigated in previous studies (Suzuki 2006; Suzuki et al. 2013). Thereafter, we use the fixed value of 1/1600 for ${f}_{\mathrm{ph}}$ by referring to Suzuki et al. (2013). Note that, when ${f}_{\mathrm{ph}}=1/1600$, the magnetic field strength at $r=1\,\mathrm{au}$ is 2.1nT, which is within the typical observed value (Wang et al. 2000). The configuration of the magnetic flux tube with ${f}_{\mathrm{ph}}$ = 1/1600 is depicted in Figure 1. As shown in this figure, the merging height H_m is the parameter defining the magnetic field strength $\overline{B}$ from the chromosphere up to the lower corona. The higher merging height corresponds to a weaker magnetic field $\overline{B}$, and, in particular, ${H}_{m}/{H}_{\mathrm{ph}}=8,12$ are used in this study. It should be noted that ${H}_{m}/{H}_{\mathrm{ph}}=8,10,12$ correspond to $\overline{B}=29,11,4$ G, respectively. These magnetic field strengths are comparable to the typical value for the area-averaged magnetic field strength in the coronal hole (3–36 G near the solar activity maximum and 1–7 G close to the minimum, according to Harvey et al. 1982; see also review by Wiegelmann & Solanki 2004). By adopting a higher coronal loop height H_l, the magnetic field strength in the upper corona can be larger (Figure 1), but H_l/r_⊙ is fixed at 0.1.

Zoom In Zoom Out Reset image size

The equation of state is $p=\rho {R}_{g}T/[{\mu }_{\mathrm{ph}}(1-\chi (T)/2)]$, where $\chi (T)$ is the ionization degree as a function of temperature, which is calculated by referring to Carlsson & Leenaarts (2012). The radiative cooling ${Q}_{\mathrm{rad}}$ is given by the empirical formulae, which is composed of three distinct terms, i.e., the photospheric radiation ${Q}_{\mathrm{ph}}$, chromospheric radiation ${Q}_{\mathrm{ch}}$, and coronal radiation ${Q}_{\mathrm{cr}}$:

$$\begin{eqnarray}&&{Q}_{\mathrm{rad}}=(1-{\xi }_{1})(1-{\xi }_{2}){Q}_{\mathrm{ph}}+{\xi }_{1}(1-{\xi }_{2}){Q}_{\mathrm{ch}}+{\xi }_{2}{Q}_{\mathrm{cr}},\end{eqnarray} \tag{ 12 }$$

where ${\xi }_{1}$ and ${\xi }_{2}$ are assumed to be as follows:

$$\begin{eqnarray}&&{\xi }_{1}=\displaystyle \frac{1}{2}\left[1+\tanh \left(\displaystyle \frac{r-{r}_{\odot }}{{H}_{\mathrm{ph}}}-3\right)\right],\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\xi }_{2}=\exp \left(-4\times {10}^{-20}{\int }_{\infty }^{r}{n}_{{\rm{H}}{\rm\small{I}}}{dr}^{\prime} \right),\end{eqnarray} \tag{ 14 }$$

${n}_{{\rm{H}}{\rm\small{I}}}=(1-\chi (T))\rho /{m}_{p}$ is the neutral hydrogen density. Each term in Equation (12) is defined as follows:

$$\begin{eqnarray}&&{Q}_{\mathrm{ph}}=4\rho {\kappa }_{{\rm{R}}}{\sigma }_{\mathrm{SB}}{T}^{4}\max \left(\displaystyle \frac{{T}^{4}}{{T}_{\mathrm{ref}}^{4}}-1,-{e}^{-{\left(r-{r}_{\odot }\right)}^{2}/{H}_{\mathrm{ph}}^{2}}\right),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\rm{where}}\ \ {T}_{\mathrm{ref}}={T}_{\mathrm{eff}}{\left(\displaystyle \frac{3}{4}\rho {\kappa }_{{\rm{R}}}{H}_{\mathrm{ph}}+\displaystyle \frac{1}{2}\right)}^{1/4},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{Q}_{{\rm{c}}{\rm{h}}}=4.9\times {10}^{9}[{\rm{e}}{\rm{r}}{\rm{g}}\,{{\rm{g}}}^{-1}\,{s}^{-1}]\rho ,\,{Q}_{{\rm{c}}{\rm{r}}}=\chi (T)(\rho /{m}_{p}{)}^{2}{\rm{\Lambda }}(T).\end{eqnarray} \tag{ 17 }$$

${\kappa }_{{\rm{R}}}=0.2$ cm² g⁻¹ pertains to the Rosseland opacity on the photosphere. ${\sigma }_{\mathrm{SB}}$ is the Stefan–Boltzmann constant. ${Q}_{\mathrm{ch}}$ and ${Q}_{\mathrm{cr}}$ are the same functions used in Hori et al. (1997), which are always positive. ${\rm{\Lambda }}(T)$ in Equation (17) is the radiative loss function for an optically thin plasma. ${Q}_{\mathrm{ph}}$ in Equation (15) is allowed to be negative where ${e}^{-{(r-{r}_{\odot })}^{2}/{H}_{\mathrm{ph}}^{2}}\sim 1$, which represents the radiative heating.

The heat conductive flux is presented by

$$\begin{eqnarray}&&{F}_{{\rm{c}}}=-\kappa (T)\displaystyle \frac{\partial T}{\partial x},\end{eqnarray} \tag{ 18 }$$

where $\kappa (T)$ is the heat conductivity as a function of the temperature. That is composed of the collisional and collisionless terms:

$$\begin{eqnarray}&&\kappa (T)=q{\kappa }_{\mathrm{coll}}+(1-q){\kappa }_{\mathrm{sat}},\end{eqnarray} \tag{ 19 }$$

where $q=\max (0,\min (1,1-0.5{\kappa }_{\mathrm{coll}}/{\kappa }_{\mathrm{sat}}))$. ${\kappa }_{\mathrm{coll}}(T)$ is adopted from Nagai (1980), which agrees with the Spitzer–Härm heat conductivity (Spitzer & Härm 1953) ${\kappa }_{0}{T}^{5/2}$ (${\kappa }_{0}={10}^{-6}$ in CGS unit) when $T\gt {10}^{6}$ K. ${\kappa }_{\mathrm{sat}}$ is given by

$$\begin{eqnarray}&&{\kappa }_{\mathrm{sat}}=\displaystyle \frac{3}{2}{{pv}}_{e,\mathrm{thr}}\displaystyle \frac{r}{T},\end{eqnarray} \tag{ 20 }$$

where ${v}_{e,\mathrm{thr}}$ is the thermal speed of the electron. ${\kappa }_{\mathrm{sat}}$ represents the saturation of heat flux caused by the collisionless effect (Parker 1964; Bale et al. 2013). The above expression of ${\kappa }_{\mathrm{sat}}$ means that the transition of heat conductivity from ${\kappa }_{\mathrm{coll}}$ to ${\kappa }_{\mathrm{sat}}$ occurs around $r\sim {\lambda }_{e,\mathrm{mfp}}$ (${\lambda }_{e,\mathrm{mfp}}$ is the electron mean free path) and that the heat flux is limited to $\displaystyle \frac{3}{2}\alpha {{pv}}_{e,\mathrm{thr}}$ in the distance where $T\sim {r}^{-\alpha }$ ($\alpha =0.2-0.4$ for winds faster than 500 km s⁻¹; Marsch et al. 1989). Based on the foregoing heat conductivity, heat conduction is solved by the super-time-stepping method (Meyer et al. 2012, 2014).

We set the static atmosphere with a temperature of 10⁴ K as the initial state. The temperature on the bottom boundary is promptly cooled down to ${T}_{\mathrm{eff}}=5770$ K after the initiation of the simulation. The mass density and poloidal magnetic field strength on the photosphere are ${\rho }_{\mathrm{ph}}=2.5\times {10}^{-7}$ g cm⁻³ and ${B}_{\mathrm{ph}}=1560$ G, respectively. To excite the outwardly propagating Alfvén wave, the toroidal velocity v_ϕ on the bottom boundary is oscillated artificially, which represents convective motion on the solar photosphere. We consider it as a frequency-dependent fluctuation with the following power spectrum:

$$\begin{eqnarray}&&{v}_{\mathrm{conv}}^{2}\propto {\int }_{{\nu }_{\min }}^{{\nu }_{\max }}{\nu }^{-1}d\nu ,\end{eqnarray} \tag{ 21 }$$

where v_conv is the free parameter corresponding to the amplitude of the convective velocity. ${\nu }_{\min }^{-1}$ and ${\nu }_{\max }^{-1}$ are 30 minutes and 20 s, respectively. The phase offsets of fluctuation are randomly assigned. The amplitude of fluctuation v_conv is the subject of survey in this study, e.g., ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}=$ 0.07, 0.14, 0.21, 0.42, 0.85 (${c}_{s\mathrm{ph}}=\sqrt{\gamma {R}_{g}{T}_{\mathrm{eff}}/{\mu }_{\mathrm{ph}}}=7.8$ km s⁻¹ is the adiabatic sound speed on the photosphere). This parameter range includes the typical velocity of horizontal convective motion, 1.1 km s⁻¹ (Matsumoto & Kitai 2010).

To excite the purely outward Alfvén waves on the bottom boundary, the toroidal magnetic field B_ϕ is determined by ${B}_{\phi }=-\sqrt{4\pi \rho }{v}_{\phi }$. This means that the Elsässer variables (i.e., ${z}_{\mathrm{out}}={v}_{\phi }-{B}_{\phi }/\sqrt{4\pi \rho }$, and ${z}_{\mathrm{in}}={v}_{\phi }+{B}_{\phi }/\sqrt{4\pi \rho }$) on the bottom boundary satisfy the conditions, i.e., ${z}_{\mathrm{out}}=2{v}_{\phi }$ and ${z}_{\mathrm{in}}=0$. The longitudinal velocity component v_x on the bottom boundary is also given as the fluctuation with the amplitude v_conv, the power spectrum similar to that of the foregoing, and the randomly assigned phase offsets. We performed a few simulations with v_x = 0 on the bottom boundary. We were able to confirm that ${v}_{x}\ne 0$ on the photosphere does not have any influence on the solar wind structure, but the spicule height can depend on it.

The upper boundary is treated as the free boundary. A total of 19,200 grids are placed nonuniformly in between. The numerical scheme is based on the HLLD Riemann solver (Miyoshi & Kusano 2005) with the second-order MUSCL interpolation and the third-order TVD Runge–Kutta method (Shu & Osher 1988).

After several tens of hours, the solar wind in the simulation box reaches quasi-steady state with numerous wave signatures (Figure 2). Figure 3 shows the simulation results, including the snapshots of solar wind velocity, mass density, temperature profiles, and temporally averaged profiles of the Alfvén wave amplitude and Alfvén speed in the solar wind. The black and red lines in each figure correspond to the results in the cases of $\overline{B}=29$ and 4 G, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The top panel of Figure 3 shows that the solar wind in the $\overline{B}=4$ G case is found to be faster than that in $\overline{B}=29$ G case. The Alfvén speed at the coronal base is much higher in the $\overline{B}=29$ G case than in $\overline{B}=4$ G. In the outer space above the coronal loop height, where the magnetic field strengths in both cases are the same, the Alfvén speed in the $\overline{B}=4$ G case is larger than in $\overline{B}=29$ G, clearly indicating denser wind in the $\overline{B}=29$ G case. With regard to the higher Alfvén speed at the coronal base in $\overline{B}=29$ G, the Alfvén speed steeply declines above the coronal loop height due to the largely expanding magnetic flux tube. This induces the strong interference between the outward and inward Alfvén waves, resulting in the humps in the Alfvén wave amplitude profiles below $0.1{r}_{\odot }$.

The most significant discrepancy between the solar winds in different merging heights is found in the wind's mass-loss rate. Figure 4 shows the mass-loss rates as a function of the energy input from the photosphere (${F}_{{\rm{A}}0}={\rho }_{\mathrm{ph}}{v}_{\mathrm{conv}}^{2}{V}_{\mathrm{Aph}}$). The filled and open circles show the results for $\overline{B}=29$ and 4 G, respectively. While the wind's mass-loss rate monotonically increases with a larger energy input from the photosphere in the case of $\overline{B}=29$ G, that in $\overline{B}=4$ G is almost independent of the energy input. The mass-loss rate in $\overline{B}=4$ G is limited to ∼10⁻¹⁵ ${M}_{\odot }$ yr ⁻¹ even in the largest energy input case of ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}=0.42$, which is two orders of magnitude smaller than that in the $\overline{B}=29$ G case.

Zoom In Zoom Out Reset image size

Figure 5 shows the time-slice diagrams of the mass density in the lower atmosphere. The top of the chromosphere ($\rho /{\rho }_{\mathrm{ph}}\sim {10}^{-7}$) shows the upward and downward motions representing the spicule dynamics. Figures 5(a) and (b) are the results in the cases of $\overline{B}=29$ G for ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}=0.21$ and 0.42, respectively.

Zoom In Zoom Out Reset image size

The height of the spicule becomes taller with a larger ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}$. On the other hand, the height of the spicule in the $\overline{B}=4$ G case is less dependent on v_conv, as shown in Figures 5(c) and (d). The average spicule height, as a function of v_conv, is summarized in Figure 6. The line styles and symbols are the same as those used in Figure 4. The spicule height is measured by tracking the isothermal contour of 4 × 10⁴ K, the typical temperature of the transition layer (Heggland et al. 2011; Iijima & Yokoyama 2015). By fitting the oscillatory pattern of the isothermal contour with the trajectories of the Lagrange particles, the individual spicules are identified, which enables us to do statistical analysis.

Zoom In Zoom Out Reset image size

A common feature can be confirmed in the behaviors of the wind's mass-loss rate (Figure 4) and the average spicule height (Figure 6). The spicule becomes monotonically taller with a larger v_conv in the $\overline{B}=29$ G case, while in $\overline{B}=4$ G case, it is almost independent of v_conv.

The lesser dependence of the simulated solar wind on v_conv implies significant wave damping below the transition layer, i.e., in the chromosphere. The difference in the spicule dynamics between $\overline{B}=4$ and 29 G also suggests that the propagation of a chromospheric shock wave is qualitatively affected by the parameter $\overline{B}$. These possibilities are further investigated in the following section.

To investigate the energy transfer by Alfvén waves, the time-averaged Poynting flux of the magnetic tension force (${F}_{{\rm{A}}}=-{B}_{\phi }{v}_{\phi }{B}_{x}/(4\pi )$) is plotted as a function of height in Figure 7. The black and red lines correspond to the results in the cases of $\overline{B}=29$ and 4 G, respectively. Although the velocity amplitude on the photosphere is fixed at ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}=0.21$, the F_A below 1 Mm in the $\overline{B}=4$ G case is slightly larger than that in the $\overline{B}=29$ G case. This is caused by the reflection of the Alfvén waves at the merging height. The energy flux of the reflected (inward) Alfvén waves is plotted in Figure 7(b) using dotted lines, where ${F}_{{\rm{A}}}^{\mathrm{out},\mathrm{in}}=\displaystyle \frac{1}{4}\rho {z}_{\mathrm{out},\mathrm{in}}^{2}{V}_{{\rm{A}}x}$, ${F}_{{\rm{A}}}={F}_{{\rm{A}}}^{\mathrm{out}}-{F}_{{\rm{A}}}^{\mathrm{in}}$. As seen in this plot, the inward Alfvén waves below 1 Mm comes mainly from the merging height, above which the Alfvén speed exponentially increases (Figure 7(c)). The energy flux of the inward Alfvén waves below 1 Mm is, therefore, related to the outward energy flux at the merging height. This leads to the smaller net energy flux when the merging height is lower.

Zoom In Zoom Out Reset image size

The most remarkable feature in Figure 7(a) is the significant decrease in the energy flux around the transition layer in the case of $\overline{B}=4$ G (red line). Figure 8(a) shows the dependence of the F_A height profile on v_conv in the case of $\overline{B}=4$ G. Although a larger v_conv produces a larger F_A on the bottom boundary (${F}_{{\rm{A}}0}=6\times {10}^{8}\mbox{--}2\times {10}^{10}$ erg cm⁻² s⁻¹ for ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}=0.07\mbox{--}0.42$), the transmitted energy fluxes into the corona do not show the significant increase from ${F}_{{\rm{A}}}\sim {10}^{5}$ erg cm⁻² s⁻¹ (${F}_{{\rm{A}}}(A/{A}_{0})\,\sim$ a few ×10⁷ erg cm⁻² s⁻¹). In other words, the additional energy input associated with a larger v_conv is completely lost below the transition layer. This cannot be seen in the case of $\overline{B}=29$ G. Figure 8(b) shows that a larger energy input from the photosphere always leads to larger transmitted energy flux when $\overline{B}=29$ G.

Zoom In Zoom Out Reset image size

In the previous subsection, it was determined that the transmission of energy flux into the corona is limited to ∼10⁵ erg cm⁻² s⁻¹ when the merging height is higher ($\overline{B}=4$ G). This suggests that Alfvén waves cannot be responsible for a Poynting flux across the chromosphere that is larger than a certain upper limit in the case of $\overline{B}=4$ G. Therefore, how the oscillations of toroidal velocity and magnetic field depend on the poloidal magnetic field configuration in the chromosphere was investigated.

Figure 9 shows the twisting motion of the magnetic flux tube in the chromosphere. Figures 9(a1) and (b1) show the time-slice diagram of density in the lower atmosphere when $\overline{B}=29$ and 4 G, respectively. Figures 9(a2) and (b2) show the nonlinearity of Alfvén wave amplitude. Because of the weaker $\overline{B}$, ${v}_{\phi }/{V}_{{\rm{A}}x}$ are higher in the $\overline{B}=4$ G case. In addition, the toroidal velocities above and below the merging height (horizontal dashed lines) often have opposite signs when $\overline{B}=4$ G. Such an antiphase oscillation is rarely seen when $\overline{B}=29$ G. This difference is more clearly seen in Figures 9(a3) and (b3). These panels show the comparison of the low-frequency component of the v_ϕ oscillation ($\nu \lt 1$ mHz). The antiphase oscillation mentioned above appears in Figure 9(b3).

Zoom In Zoom Out Reset image size

Figures 10 and 11 depict the typical time sequence of the magnetic field lines in the cases of $\overline{B}=4$ and 29 G, respectively. When the merging height is low and $\overline{B}$ is large, the upper part of the flux tube above the merging height is twisted as its lower part rotates (Figure 11). On the other hand, Figure 10 shows that the upper part of the flux tube is counterrotating against the lower part, thereby causing the formation of the break in the magnetic field line. The close-up view around such a break in the magnetic field line is shown in Figure 12, which corresponds to the rectangle area in Figure 9. The break in the magnetic field line is represented by the dashed line in this figure, which agrees with the characteristics at ${v}_{x}+{B}_{x}/\sqrt{4\pi \rho }$. Figure 9 shows that this signature appears transiently and is associated with a compression that is strong enough to significantly enhance the plasma β in the downstream. The break in the magnetic field line is, therefore, identified as the intermediate shock.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The previous subsections revealed that energy transfer by Alfvén waves is restricted in the case of a weak magnetic field ($\overline{B}=4$ G). Aside from such a nearly incompressible wave, the propagation of magnetoacoustic shocks, including slow and fast shocks, is possibly dependent on the magnetic field configuration in the chromosphere. In fact, Figure 6 shows the dependence of the average spicule height on v_conv changes in accordance with the magnetic field strength $\overline{B}$. For a comprehensive discussion, we investigated the propagation of slow and fast shocks.

The relatively strong compressible wave can be distinguished as propagating spiky signatures with $-{\partial }_{x}{v}_{x}\gt 0$. After tracing these signatures, the Alfvén Mach number of the shock wave (M_A) is calculated using the following formula (the derivation is described in Appendix A):

$$\begin{eqnarray}&&{M}_{{\rm{A}}}=-\displaystyle \frac{1}{{V}_{{\rm{A}}}}{\left(\displaystyle \frac{\partial {v}_{x}}{\partial x}\right)}^{-1}\left\{\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial {p}_{\mathrm{tot}}}{\partial x}-\displaystyle \frac{\partial }{\partial x}\left(\displaystyle \frac{{{GM}}_{\odot }}{r}\right)\right\},\end{eqnarray} \tag{ 22 }$$

where ${p}_{\mathrm{tot}}=p+{B}_{\phi }^{2}/(8\pi )$. By expressing the fast- and slow-mode Mach numbers with ${M}_{f}={M}_{{\rm{A}}}{V}_{{\rm{A}}}/{V}_{\mathrm{fast}}$ and ${M}_{s}\,={M}_{{\rm{A}}}{V}_{{\rm{A}}}/{V}_{\mathrm{slow}}$, where ${V}_{\mathrm{fast}}$ and ${V}_{\mathrm{slow}}$ are the fast- and slow-mode speeds, the detected shock is specified as the fast shock when $| {M}_{f}-1| \lt | {M}_{s}-1|$, or, otherwise, the slow shock. This classification is justified when both fast and slow shocks are relatively weak, i.e., ${M}_{f}\sim 1$ and ${M}_{s}\sim 1$. By counting the fast (slow) shocks with M_f (M_s) propagating around the mass density ρ in the stratified atmosphere, the distribution function of M_f or M_s, with respect to ρ, is defined as follows:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dN}}{d{\mathrm{log}}_{10}\rho {dM}}(\rho ,M)\\ \quad =\displaystyle \frac{{dN}(\rho \in [{\rho }_{i},{\rho }_{i+1}],M\in [{M}_{j},{M}_{j+1}])}{({\mathrm{log}}_{10}{\rho }_{i+1}-{\mathrm{log}}_{10}{\rho }_{i})({M}_{j+1}-{M}_{j})},\end{array}\end{eqnarray} \tag{ 23 }$$

where ${dN}(\rho ,M)$ is the expected number of shocks characterized with ($\rho ,M$) in one snapshot and the subscriptions i and j represent the discretization.

Figure 13 shows the distribution functions calculated from the simulation results in the cases of $\overline{B}=29$ G (upper panels) and 4 G (lower panels). ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}$ is fixed at 0.21. The vertical dotted line in each panel corresponds to the mean mass density at the transition layer. The distribution around the transition layer is artificially sparse in all panels. This is because the shock crossing the transition layer is hardly detected in this analysis (we used the time series data over 50,000 s with an interval of 4 s. This interval is much longer than that of the shock-crossing timescale across the transition layer). The cross symbol represents the most frequently appearing Mach number in each bin of ρ. Thus, the gradual rise of cross symbols seen around $\rho \gtrsim {10}^{-12}$ g cm⁻³ in Figure 13(a2) shows the growth of the chromospheric slow shock. Note that the coronal slow shocks in Figures 13(a2) and (b2) are concentrated on ${M}_{s}\sim 1/\sqrt{\gamma }$ rather than ${M}_{s}\sim 1$. This is because we calculated M_s by assuming γ = 5/3 without any considerations of nonadiabatic effects. The phase speed of slow shock in the corona tends to be the isothermal sound speed $\sim \sqrt{p/\rho }$ due to the strong heat conduction. This leads to the underestimation of M_s by a factor of $\sim 1/\sqrt{\gamma }$ in the corona.

Zoom In Zoom Out Reset image size

The most remarkable feature in this figure is that the slow shock vanishes around $\rho \sim {10}^{-11}$ g cm⁻³ in Figure 13(b2). $\rho ={10}^{-11}$ g cm⁻³ is three orders of magnitude higher than the mass density at the transition layer and roughly corresponds to the mass density around the merging height. Therefore, this disappearance of slow shock is not related to the above-mentioned artificial sparse distribution around the transition layer. Instead, it is implied that the slow shock can be evanescent in the chromosphere when the magnetic field is weak.

We discuss here the causal relationship between the high nonlinearity of Alfvén waves in the chromosphere and the limit on the energy transmission into the corona (Section 4.1).

Figure 14 shows the temporally averaged profiles of nonlinearity regarding B_ϕ, v_ϕ, and v_x. Note that we abbreviated $\langle {({B}_{\phi }/{B}_{x})}^{2}{\rangle }^{1/2}$, $\langle {({v}_{\phi }/{V}_{{\rm{A}}x})}^{2}{\rangle }^{1/2}$, and $\langle {({v}_{x}/{V}_{{\rm{A}}x})}^{2}{\rangle }^{1/2}$ into ${B}_{\phi }/{B}_{x}$, ${v}_{\phi }/{V}_{{\rm{A}}x}$, and ${v}_{x}/{V}_{{\rm{A}}x}$. ${B}_{\phi }/{B}_{x}$ and ${v}_{\phi }/{V}_{{\rm{A}}x}$ represent the nonlinearity of Alfvén waves. The maximum level of nonlinearity is always found around the merging height. The higher merging height is responsible for the higher maximum nonlinearity of torsional (B_ϕ and v_ϕ) and longitudinal (v_x) oscillation. In particular, the high level of nonlinearity of v_x corresponds to the large inertia of the magnetic flux tube.

Zoom In Zoom Out Reset image size

By using mass conservation (Equation (1)) and poloidal magnetic flux conservation (Equation (6)), the toroidal component of the equation of motion (Equation (4)) can be expressed as below:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{\phi }}{\partial t}+\displaystyle \frac{{v}_{x}}{\sqrt{A}}\displaystyle \frac{\partial (\sqrt{A}{v}_{\phi })}{\partial x}-\displaystyle \frac{{B}_{x}}{4\pi \rho \sqrt{A}}\displaystyle \frac{\partial (\sqrt{A}{B}_{\phi })}{\partial x}=0.\end{eqnarray} \tag{ 24 }$$

The variables in the above equation are the same as those used in Section 2.1. The second term represents the inertia term. It can be competitive against the Lorentz force (the third term) when the longitudinal oscillation is highly nonlinear. The last panel of Figure 14 shows the ratio of the temporally averaged absolute value of the inertia term in relation to that of the restoring term in Equation (24). Here, we define the following:

$$\begin{eqnarray}&&{f}_{\mathrm{inertia}}=\left\langle \left|\displaystyle \frac{{v}_{x}}{\sqrt{A}}\displaystyle \frac{\partial (\sqrt{A}{v}_{\phi })}{\partial x}\right|\right\rangle \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{f}_{\mathrm{restoring}}=\left\langle \left|\displaystyle \frac{{B}_{x}}{4\pi \rho \sqrt{A}}\displaystyle \frac{\partial (\sqrt{A}{B}_{\phi })}{\partial x}\right|\right\rangle .\end{eqnarray} \tag{ 26 }$$

The low ratio in the corona means that Alfvén waves can propagate without a significant nonlinear effect while the ratio around unity implies that wave propagation is strongly affected by the inertia term. In the case of the higher merging height ($\overline{B}=4$ G), the ratio reaches unity around the merging height, as shown with the red thick line in Figure 15. Due to this large inertia of the magnetic flux tube, the rotation of the upper part of the flux tube cannot be restored so easily by the twisting motion injected from the photosphere. That would result in the antiphase oscillation between the upper and lower parts of the flux tube (Figure 9). The highly sheared torsional flow and nonlinear longitudinal oscillation can cause the "fracture" of the flux tube, i.e., the formation of the intermediate shock (Figure 10). Once the intermediate shock is formed in the chromosphere, the Poynting flux associated with it hardly transmits into the corona. That is because the intermediate shock easily interacts with the slow and fast shocks or contact discontinuity, including the transition layer itself. Among these interactions, the collision of intermediate shock with the transition layer results in the transmitted waves composed of fast rarefaction wave and slow shock. As both of them have negative Poynting fluxes, the magnetic energy transferred by the chromospheric intermediate shock is, in this sense, confined below the transition layer until its dissipation.

Zoom In Zoom Out Reset image size

Figures 15 and 16 show examples of the formation of the chromospheric intermediate shock. In Figure 15, the intermediate shock deviates from the fast shock at around t = 120 s and collides with the downward slow shock at around t = 340 s. Although the upward fast shock generated by this collision has a positive Poynting flux, the other resultant waves, including the upward intermediate rarefaction wave, transport the magnetic energy downward. The formation of the intermediate shock in Figure 16 is a result of the head-on collision of upward and downward slow shocks, which is associated with the encounter with large shear flow. The upward intermediate shock finally becomes bidirectional fast shocks after the interaction with the other waves. The dissipation of the intermediate shock is clearly exemplified in Figure 17. In this scene, the interaction between the sequence of intermediate shocks and the downward slow shock results in the bidirectional slow shocks. As a result, the highly sheared magnetic field line is rapidly relaxed and the super-Alfvénic torsional flow is generated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 14 shows that wave nonlinearity such as ${B}_{\phi }/{B}_{x}$, ${v}_{\phi }/{V}_{{\rm{A}}x}$, and ${v}_{x}/{V}_{{\rm{A}}x}$ is highest around the merging height. It demonstrates that higher merging height (or weaker $\overline{B}$) and larger ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}$ are always associated with higher wave nonlinearity in the chromosphere. By focusing on the maximum values of the profiles plotted in Figure 14, the scaling relations between ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}$ and wave nonlinearity were summarized in Figure 18.

Zoom In Zoom Out Reset image size

Figure 18(a) shows that the ratio of the inertia force to the restoring force is clearly correlated to ${({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}})(\overline{B}/{B}_{\mathrm{ph}})}^{-1/2}$, i.e.,

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\mathrm{inertia}}/{f}_{\mathrm{restoring}}=0.28\left[({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}})\right.\\ \quad {\left.\times {\left(\overline{B}/{B}_{\mathrm{ph}}\right)}^{-1/2}\right]}^{-0.89}.\end{array}\end{eqnarray} \tag{ 27 }$$

This scaling is composed of the relation between ${f}_{\mathrm{inertia}}/{f}_{\mathrm{restoring}}$ with ${v}_{x}/{V}_{{\rm{A}}x}$ (Figure 18(b)) and that between ${v}_{x}/{V}_{{\rm{A}}x}$ with ${({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}})(\overline{B}/{B}_{\mathrm{ph}})}^{-1/2}$ (Figure 18(c)). In fact, Equations (25) and (26) indicate that ${f}_{\mathrm{inertia}}\sim {v}_{x}{v}_{\phi }/{\lambda }_{{\rm{A}}}$ and ${f}_{\mathrm{restoring}}\,\sim {B}_{x}{B}_{\phi }/(4\pi \rho {\lambda }_{{\rm{A}}})$, where ${\lambda }_{{\rm{A}}}$ is the wavelength of the Alfvén waves. Therefore, ${f}_{\mathrm{inertia}}/{f}_{\mathrm{restoring}}$ tends to be ${v}_{x}/{V}_{{\rm{A}}x}$ when ${v}_{\phi }\sim {B}_{\phi }/\sqrt{4\pi \rho }$. The amplitude of v_x basically follows the energy flux conservation for the longitudinal wave in the isothermal atmosphere. That means $\rho {v}_{x}^{2}/{B}_{x}\sim$ constant and ${v}_{x}/{V}_{{\rm{A}}x}\propto {B}_{x}^{-1/2}$. By using these scaling relations, it is inferred that there is a critical ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}$ or $\overline{B}/{B}_{\mathrm{ph}}$ across which the chromosphere is too highly nonlinear such that the Lorentz force associated with Alfvén wave propagation (${f}_{\mathrm{restoring}}$) can no longer twist the flux tube against the large inertia force (${f}_{\mathrm{inertia}}$). From Equation (27), we replace this critical condition of ${f}_{\mathrm{inertia}}\gtrsim {f}_{\mathrm{restoring}}$ with $0.28{[({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}})/\sqrt{\overline{B}/{B}_{\mathrm{ph}}}]}^{-0.89}\gtrsim 1$ or

$$\begin{eqnarray}&&{v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}\gtrsim 4.2\sqrt{\overline{B}/{B}_{\mathrm{ph}}}.\end{eqnarray} \tag{ 28 }$$

This implies the following: first, for a given flux tube with $\overline{B}/{B}_{\mathrm{ph}}$, the energy input from the photosphere larger than ${F}_{{\rm{A}},\mathrm{cr}}={\rho }_{\mathrm{ph}}{v}_{\mathrm{conv},\mathrm{cr}}^{2}{V}_{\mathrm{Aph}}$ does not contribute to the coronal heating. As such,

$$\begin{eqnarray}&&{F}_{{\rm{A}},\mathrm{cr}}=2.4\times {10}^{12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}(\overline{B}/{B}_{\mathrm{ph}}).\end{eqnarray} \tag{ 29 }$$

When $\overline{B}=4$ G and ${B}_{\mathrm{ph}}=1560$ G, we find ${({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}})}_{\mathrm{cr}}\sim 0.21$ and ${F}_{{\rm{A}},\mathrm{cr}}=6.3\times {10}^{9}$ erg cm⁻² s⁻¹. Second, for a given convection velocity of ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}$, the magnetic flux tube with $\overline{B}/{B}_{\mathrm{ph}}\lt {(\overline{B}/{B}_{\mathrm{ph}})}_{\mathrm{cr}}=0.057{({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}})}^{2}$ is unable to guide the magnetic energy from the lower atmosphere to the corona. When ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}=0.21$ and ${B}_{\mathrm{ph}}=1560$ G, we find ${\overline{B}}_{\mathrm{cr}}=4$ G.

Finally, it is notable that the wave nonlinearity of ${B}_{\phi }/{B}_{x}$ follows ${B}_{\phi }/{B}_{x}\propto {({\overline{B}}^{-1/2})}^{1.15}$ (Figure 18(d)). This is accounted for by the energy flux conservation for the Alfvén waves propagating along the magnetic flux tube that expands like ${B}_{x}^{2}\propto \rho$. Because $\rho {v}_{\phi }^{2}{V}_{{\rm{A}}x}A=$ const., we find ${v}_{\phi }\propto {\rho }^{-1/4}$ and ${B}_{\phi }/{B}_{x}\,\sim {v}_{\phi }/{V}_{{\rm{A}}x}\propto {\rho }^{1/4}/{B}_{x}$. Thus, when ${B}_{x}^{2}\propto \rho$, it is obtained that ${B}_{\phi }/{B}_{x}\propto {B}_{x}^{-1/2}$. Hollweg (1971) and Shibata & Uchida (1985) discussed that large-amplitude Alfvén waves can be responsible for the longitudinal motion and derived the relationship of ${v}_{x}/{V}_{{\rm{A}}x}\propto {({B}_{\phi }/{B}_{x})}^{2}$. On the other hand, Figures 18(c) and (d) show ${v}_{x}/{V}_{{\rm{A}}x}=0.29{[({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}})/\sqrt{(\overline{B}/{B}_{\mathrm{ph}})}]}^{0.86}$ and ${B}_{\phi }/{B}_{x}\,=0.52{\left[\left({v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}\right)/\sqrt{\left(\overline{B}/{B}_{\mathrm{ph}}\right)}\right]}^{1.15}$, leading to the following scaling law:

$$\begin{eqnarray}&&{v}_{x}/{V}_{{\rm{A}}x}=0.51{\left({B}_{\phi }/{B}_{x}\right)}^{0.75}.\end{eqnarray} \tag{ 30 }$$

Figure 16 exhibits the formation of the intermediate shock as well as the disappearance of the upward slow shock. This reminds us of Figure 13, which shows that the slow shock is absent in the upper chromosphere when the magnetic field is weak. In addition to this head-on collision of the counterpropagating slow shocks, the head-on and rear-end collisions between the slow and intermediate shocks can disturb the upward propagation of the slow shock. These interactions would be encouraged in the highly nonlinear chromosphere, especially when the magnetic field is weak. This is because the crossing timescale of slow shock at a speed of $\sim {c}_{s}{B}_{x}/B$ becomes longer as nonlinearity increases.

This evanescence of the slow shock in the chromosphere could result in the following two consequences about the spicule dynamics. First, the ejection speed of the spicule would become smaller and less dependent on v_conv. Second, less frequent slow shocks in the upper chromosphere could reduce the chromospheric temperature, leading to a shorter density scale height in the chromosphere (Appendix B). As a result of the smaller ejection speed and shorter density scale height, the average spicule height in the weak magnetic field tends to be lower than that in the strong magnetic field.

5.4.1. Solar Wind

The typical fast solar wind proton flux observed around 1 au is ∼2 × 10⁸ cm⁻² s⁻¹ (Withbroe 1989; Wang 2010), comparable to the simulated value in the stronger magnetic field case ($\sim 2.1\times {10}^{8}$ cm⁻² s⁻¹ for $\overline{B}=29$ G), but inconsistent with that in the weaker magnetic field case ($\sim 0.20\times {10}^{8}$ cm⁻² s⁻¹ for $\overline{B}=4$ G). As discussed in Section 5.1, Poynting flux into the corona is limited to 10⁵ erg cm⁻² s⁻¹ when $\overline{B}=4$ G, which causes a significantly low-mass flux of solar wind. The inconsistency between the observed and simulated mass-loss rates in the $\overline{B}=4$ G case is, however, easily solved by considering the polarization of the Alfvén waves. In the present study, we used the axisymmetric coordinate system with linearly polarized Alfvén waves. The nonlinear propagation of circularly polarized Alfvén waves in nonsteady solar wind was simulated by Suzuki & Inutsuka (2006) and Shoda et al. (2018) using the local spherically symmetric coordinate system (Shoda & Yokoyama 2018). The differences between the axisymmetric and local spherically symmetric coordinate systems are summarized in Appendix C. We also conducted a similar parameter survey on the solar atmosphere and wind structure in the local spherically symmetric coordinate system with circularly polarized Alfvén waves. Consequently, it is confirmed that the results qualitatively agree with those in the axisymmetric coordinate system. That means, even in the local spherically symmetric coordinate system, there is an upper limit on the transmitted Poynting flux into the corona when the merging height is higher ($\overline{B}=4$ G). The wind's mass-loss rate and average spicule height become independent of the velocity amplitude on the photosphere (v_conv) when ${v}_{\mathrm{conv}}/{c}_{s\mathrm{ph}}\gtrsim 4.2\sqrt{\overline{B}/{B}_{\mathrm{ph}}}$ (Section 5.2). Figure 19 shows the mass-loss rates of solar wind as a function of the energy input from the photosphere. The filled and open circles are the same as those in Figure 4, and the square symbols are overplotted as the results in the local spherically symmetric coordinate system. The wind's mass-loss rate in the local spherically symmetric coordinate system with $\overline{B}=4$ G (open squares) appear to be constant for ${F}_{{\rm{A}}0}\gtrsim {10}^{10}$ erg cm⁻² s⁻¹. The upper limit of the mass-loss rate simulated in the local spherically symmetric coordinate system is, however, much higher than that in the axisymmetric coordinate system. This is partly because the circularly polarized Alfvén waves transfer twice as much magnetic energy as the linearly polarized Alfvén waves when their amplitudes are the same. In other words, the critical Poynting flux ${F}_{{\rm{A}},\mathrm{cr}}$ (Equation (29)) in the local spherically symmetric coordinate system is calculated as ${F}_{{\rm{A}},\mathrm{cr}}=2\,\times {\rho }_{\mathrm{ph}}{v}_{\mathrm{conv},\mathrm{cr}}^{2}{V}_{\mathrm{Aph}}$,

$$\begin{eqnarray}&&{F}_{{\rm{A}},\mathrm{cr}}=4.8\times {10}^{12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}(\overline{B}/{B}_{\mathrm{ph}}).\end{eqnarray} \tag{ 31 }$$

As a result, the upper limit of the transmitted Poynting flux into the corona is much larger than that in the axisymmetric coordinate system. As such, the resultant mass-loss rate can reach the observed level. Therefore, the above-mentioned inconsistency between the observed and simulated mass-loss rates in the axisymmetric coordinate system is merely an intrinsic problem of our 1D approximation.

Zoom In Zoom Out Reset image size

Suzuki et al. (2013) reported their simulation results for solar and stellar winds, which showed that the wind's mass-loss rate saturates due to the enhanced radiative loss in the corona. From the time-steady energy equation, they paid attention to the following energy conservation law (in the local spherically symmetric coordinate system):

$$\begin{eqnarray}&&\begin{array}{l}{\left.\left[\rho {v}_{r}A\left(\displaystyle \frac{{v}_{\perp }^{2}}{2}+\displaystyle \frac{{B}_{\perp }^{2}}{4\pi \rho }\right)-{B}_{r}A\displaystyle \frac{{v}_{\perp }{B}_{\perp }}{4\pi }\right]\right|}_{{r}_{\mathrm{tc}}}\\ \ \ \ \approx {\left.\left[\rho {v}_{r}A\displaystyle \frac{{v}_{r}^{2}}{2}\right]\right|}_{\infty }+{\displaystyle \int }_{{r}_{\mathrm{tc}}}^{\infty }{{AQ}}_{\mathrm{rad}}{dr}+{\left.\left[\rho {v}_{r}A\displaystyle \frac{{{GM}}_{\odot }}{r}\right]\right|}_{{r}_{\mathrm{tc}}},\end{array}\end{eqnarray} \tag{ 32 }$$

where ${v}_{\perp }$ and ${B}_{\perp }$ are the transverse components of the velocity and magnetic field. ${r}_{\mathrm{tc}}$ represents the top of the chromosphere, the position with the temperature $T=2\times {10}^{4}$ K, according to the definition by Suzuki et al. (2013). The foregoing expression means that the Poynting flux at $r={r}_{\mathrm{tc}}$ (left-hand side) is converted to the kinetic energy of wind (the first term in the right-hand side) as well as the radiative loss and gravitational potential energy (the second and third terms on the right-hand side, respectively). While the kinetic energy of wind is positively correlated to the Alfvén waves energy at the top of the chromosphere, the large energy transmission into the corona can make radiative energy loss dominant over the kinetic energy term. That leads to the saturation of the wind's mass-loss rate. This kind of saturation is also seen in our simulation. Figure 20 is analogous to Figure 8 in Suzuki et al. (2013). It presents the comparison between the left-hand side of Equation (32) (${({L}_{{\rm{A}}}f)}_{\mathrm{tc}}$) with the first and second terms on the right-hand side (${L}_{{\rm{K}},\mathrm{out}}$ and ${({L}_{{\rm{R}}}f)}_{\mathrm{tc}}$). As seen in this figure, ${({L}_{{\rm{A}}}f)}_{\mathrm{tc}}$ larger than ∼4 × 10²⁸ erg s⁻¹ leads to the saturation of ${L}_{{\rm{K}},\mathrm{out}}$, which is associated with the enhanced ${({L}_{{\rm{R}}}f)}_{\mathrm{tc}}$. The saturation level of ${L}_{{\rm{K}},\mathrm{out}}$ is almost consistent with that suggested by Suzuki et al. (2013) as shown below, indicated by the blue horizontal line in Figure 20,

$$\begin{eqnarray}&&{L}_{{\rm{K}},\mathrm{out},\mathrm{sat}}=2.05\times {10}^{28}\ {{\rm{erg}}{\rm{s}}}^{-1}\ {\left({B}_{\mathrm{ph}}{f}_{\mathrm{ph}}\right)}^{1.84}.\end{eqnarray} \tag{ 33 }$$

This saturation is, however, not expected in the case of the higher merging height ($\overline{B}=4$ G). That is because the transmission of the Alfvén wave energy itself is limited due to its high nonlinearity in the chromosphere, as discussed in the previous subsection. This is why the open circles and squares are absent above a certain level of ${({L}_{{\rm{A}}}f)}_{\mathrm{tc}}$, indicated by the vertical thick lines in Figure 20.

Zoom In Zoom Out Reset image size

5.4.2. Spicule

As for the chromospheric intermediate shock, Snow & Hillier (2019) found that the decoupling of the neutral fluid against plasma can cause the intermediate shock when reconnection occurs in the partially ionized plasma. Our study and their study suggest that the chromospheric intermediate shock would be observed ubiquitously over a wide range of spatial scales in near future. The effect of partially ionized plasma can appear especially for the propagation of high-frequency Alfvén waves (Soler et al. 2019) and should be considered in future studies. Several limitations should be imposed on the application of the results of our study with regard to the real solar atmosphere and wind. The present study is based on a 1D approximation (symmetry assumption), flux tube model, and simplified radiation. Our study ignores solar rotation, collisionless effects, and various wave dissipation mechanisms, including phase mixing and turbulent dissipation (Cranmer et al. 2007; Shoda et al. 2018). Nevertheless, it is worthwhile to emphasize that highly nonlinear Alfvén waves in the chromosphere could restrict the energy transfer from the photosphere to the corona. As such, our findings highlight the importance of the magnetic field configuration in the chromosphere in terms of the diversity of both solar and stellar atmosphere and wind structures.