Anisotropic Magnetohydrodynamic Turbulence Driven by Parametric Decay Instability: The Onset of Phase Mixing and Alfvén Wave Turbulence

Since the discovery of the high-temperature solar atmosphere and supersonic outflow from the Sun, coronal heating and solar wind acceleration have become one of the most important research subjects in astrophysics. The Parker Solar Probe (Fox et al. 2016) is going to be launched to resolve these long-standing problems in solar physics.

In the open field regions, the wave/turbulence model is promising as the heating mechanism (Cranmer 2012). Several photospheric observations have shown that transverse waves have sufficient energy flux at the photosphere (Fujimura & Tsuneta 2009; Chitta et al. 2012), and a certain portion of these waves propagate through the chromosphere into the corona to power the solar wind (De Pontieu et al. 2007; McIntosh et al. 2011). Nonthermal coronal line broadening (Banerjee et al. 2009; Hahn & Savin 2013) also indicates that there are sufficiently strong Alfvén waves in the corona.

While the global energetics of the corona and solar wind are clarified by wave observations, the energy cascading mechanism is still under investigation. We note that, because there is a large gap between the energy-containing scale and the dissipation scale of Alfvén waves in the corona and solar wind, the cascading rate gives the approximate plasma heating rate. There are three main promising cascading mechanisms: Alfvén wave turbulence (AWT; Dobrowolny et al. 1980; Perez & Chandran 2013), phase mixing (PM; Heyvaerts & Priest 1983; Magyar et al. 2017), and parametric decay instability (PDI; Suzuki & Inutsuka 2005; Tenerani & Velli 2013).

AWT is triggered by counter-propagating Alfvén waves (Kraichnan 1965; Goldreich & Sridhar 1995). In the corona and solar wind, partial reflection due to the field-aligned gradient of Alfvén speed (Velli 1993; Verdini & Velli 2007) can drive AWT. Some theoretical models reproduce the corona and solar wind self-consistently based on the AWT heating scenario (Cranmer et al. 2007; Verdini et al. 2010; Chandran et al. 2011). However, recent 3D simulations indicate that heating by AWT is insufficient (Perez & Chandran 2013; van Ballegooijen & Asgari-Targhi 2017).

PM works when the Alfvén speed is inhomogeneous in the direction perpendicular to the mean magnetic field (Heyvaerts & Priest 1983). As the horizontal density variation is observed in the corona (Raymond et al. 2014), PM plays a role in the coronal heating. In fact, a similar process called resonant absorption is observed to work in the solar atmosphere (Okamoto et al. 2015). A recent magnetohydrodynamic (MHD) simulation has shown that PM can generate a turbulent structure, and this process is called generalized phase mixing (Magyar et al. 2017).

PDI is an instability of the Alfvén wave (Sagdeev & Galeev 1969; Goldstein 1978), in which an Alfvén wave decays into a forward slow-mode wave and backward Alfvén wave in a low-beta regime. PDI can be a source of coronal heating because slow shocks are generated, heating up the corona and solar wind (Suzuki & Inutsuka 2005; Shoda et al. 2018a). Several processes in the solar wind are possibly attributable to PDI, such as the cross-helicity evolution (Malara et al. 2000; Shoda & Yokoyama 2016), density fluctuation (Bowen et al. 2018; Shoda et al. 2018b), and power spectrum formation (Chandran 2018). Both numerical simulation (Suzuki & Inutsuka 2005; Shoda et al. 2018a) and radio-wave observation (Miyamoto et al. 2014) show large density fluctuations in the wind acceleration region, possibly generated by PDI. These results indicate that the plasmas in that region are in a state of PDI-driven turbulence.

In this Letter we conduct a simulation of PDI-driven turbulence to clarify which of the processes mentioned above is the most responsible for energy cascading (and dissipation). Multi-dimensional simulations of PDI-driven turbulence have already been performed, focusing mainly on the spectral property (Ghosh & Goldstein 1994) or the dependence on dimensions (Del Zanna et al. 2001), although no simulations have ever been performed focusing on the cascading process. We perform a detailed analysis to reveal the dominant process in each phase of instability and turbulence.

2.1. Basic Equations and Solver

We numerically solve the ideal MHD equations

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\rho +{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(\rho {\boldsymbol{v}})+{\rm{\nabla }}\cdot \left({p}_{T}\hat{{\boldsymbol{I}}}+\rho {\boldsymbol{vv}}-\displaystyle \frac{1}{4\pi }{\boldsymbol{BB}}\right)=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\boldsymbol{B}}+{\rm{\nabla }}\cdot ({\boldsymbol{vB}}-{\boldsymbol{Bv}})=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}e+{\rm{\nabla }}\cdot \left[(e+{p}_{T}){\boldsymbol{v}}-\displaystyle \frac{1}{4\pi }{\boldsymbol{B}}({\boldsymbol{v}}\cdot {\boldsymbol{B}})\right]=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&e=\displaystyle \frac{p}{{\rm{\Gamma }}-1}+\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi },\ \ \ {p}_{T}=p+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi }.\end{eqnarray} \tag{ 5 }$$

Γ = 1.0001 is used instead of the adiabatic value Γ = 5/3 to keep the system almost isothermal. Because we are interested in the cascading process, we exclude the diffusion terms; mechanical energy is thermalized by numerical diffusion at the grid scale.

The HLLD Riemann solver (Miyoshi & Kusano 2005) with WENOZ reconstruction (Borges et al. 2008) and the SSP Runge–Kutta method (Shu & Osher 1988) are used for numerical calculation. This is a high-resolution (fifth-order-in-space and third-order-in-time) shock-capturing numerical scheme, which is appropriate for the simulation of PDI. Numerical ${\rm{\nabla }}\cdot {\boldsymbol{B}}$ error is removed with the hyperbolic cleaning method (Dedner et al. 2002).

2.2. Numerical Setup

We solve the basic Equations (1)–(5) in a rectangular simulation box $[0,L]\times [0,L]\times [0,3L]$. We apply $288\,\times 288\times 864$ grid points to resolve the computational domain, so that the grid size is isotropic. The periodic boundary conditions are used for all of the boundaries. The initial condition is as follows:

$$\begin{eqnarray}&&\rho ={\rho }_{0}+\delta {\rho }_{0},p={c}_{s}^{2}{\rho }_{0},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{v}_{x}=\eta {v}_{A0}\sin ({k}_{0}z-{\omega }_{0}t),{B}_{x}=-{B}_{0}{v}_{x}/{v}_{A0}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{v}_{y}=\eta {v}_{A0}\cos ({k}_{0}z-{\omega }_{0}t),{B}_{y}=-{B}_{0}{v}_{y}/{v}_{A0}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{v}_{z}=0,{B}_{z}={B}_{0}.\end{eqnarray} \tag{ 9 }$$

Note that the box-averaged density, momentum, magnetic field, and total energy are time-independent because of the periodic boundary. ${k}_{0}=2\pi /L$ and ${\omega }_{0}={v}_{A0}{k}_{0}$ are the wavenumber and frequency of the initial wave, respectively. The parent wave period is therefore given as ${\tau }_{0}=2\pi /{\omega }_{0}=L/{v}_{A0}$. c_s and ${v}_{A0}={B}_{0}/\sqrt{4\pi {\rho }_{0}}$ denote the background sound and Alfvén velocities, respectively. η is the normalized amplitude of the initial wave. In this Letter, we fix $\eta =0.2$ and ${c}_{s}/{v}_{A0}=0.2$.

The density fluctuation is imposed to trigger PDI. We calculated two cases: a 3D setting and a 1D setting. In the 3D setting, the fluctuations are given as

$$\begin{eqnarray}&&\delta {\rho }_{0}/{\rho }_{0}=\varepsilon \sin ({k}_{0}x)\sin ({k}_{0}y)N(z),\end{eqnarray} \tag{ 10 }$$

where ε ∼ 10⁻³ is the initial amplitude of density fluctuation and N(z) denotes a random noise function of z. The 1D setting is defined as

$$\begin{eqnarray}&&\delta {\rho }_{0}/{\rho }_{0}=\displaystyle \frac{1}{2}\varepsilon N(z).\end{eqnarray} \tag{ 11 }$$

The factor 1/2 is to make the rms value of $\delta {\rho }_{0}$ the same as that in the 3D setting.

First, we discuss the time evolution of volume-averaged characteristic quantities: the normalized rms density fluctuation $\delta {\rho }_{\mathrm{rms}}/{\rho }_{0}$, the rms angle between current and magnetic field ϕ_j,B, Elsässer energies ${E}^{\pm }$ and normalized cross-helicity σ_c. $\delta {\rho }_{\mathrm{rms}}$ is defined as

$$\begin{eqnarray}&&\delta {\rho }_{\mathrm{rms}}=\sqrt{\langle {(\rho -{\rho }_{0})}^{2}\rangle }={(\rho -{\rho }_{0})}_{\mathrm{rms}},\end{eqnarray} \tag{ 12 }$$

where $\langle X\rangle$ and ${X}_{\mathrm{rms}}$ denote the volume-average and rms operators, respectively:

$$\begin{eqnarray}&&\langle X\rangle =\displaystyle \frac{1}{V}\int d{\boldsymbol{x}}X,\ \ {X}_{\mathrm{rms}}=\sqrt{\displaystyle \frac{1}{V}\int d{\boldsymbol{x}}{X}^{2}},\end{eqnarray} \tag{ 13 }$$

where $V=3{L}^{3}$ denotes the volume of the simulation box. ϕ_j,B is given by

$$\begin{eqnarray}&&{\phi }_{j,B}={\left[{\arccos }^{-1}\left(\displaystyle \frac{{\boldsymbol{j}}\cdot {\boldsymbol{B}}}{| {\boldsymbol{j}}| | {\boldsymbol{B}}| }\right)\right]}_{\mathrm{rms}},\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{j}}=\displaystyle \frac{c}{4\pi }{\rm{\nabla }}\times {\boldsymbol{B}}.\end{eqnarray} \tag{ 15 }$$

Normalized Elsässer energies are defined as:

$$\begin{eqnarray}&&{E}^{\pm }=\left\langle \displaystyle \frac{1}{4}\rho {{{\boldsymbol{\zeta }}}^{\pm }}^{2}\right\rangle /{E}_{0},\end{eqnarray} \tag{ 16 }$$

where ${E}_{0}={\rho }_{0}{v}_{A0}^{2}$ and ${{\boldsymbol{\zeta }}}^{\pm }={{\boldsymbol{v}}}_{\perp }\mp {{\boldsymbol{B}}}_{\perp }/\sqrt{4\pi \rho }$ are Elsässer variables (Elsässer 1950). The subscript ⊥ denotes the components perpendicular to the mean field ${B}_{0}{{\boldsymbol{e}}}_{z}$. The normalized cross-helicity σ_c is calculated from ${E}^{\pm }$ as (Del Zanna et al. 2001)

$$\begin{eqnarray}&&{\sigma }_{c}=\displaystyle \frac{{E}^{+}-{E}^{-}}{{E}^{+}+{E}^{-}}.\end{eqnarray} \tag{ 17 }$$

In Figure 1, we show the time evolution of the volume-averaged quantities defined above by solid lines. The corresponding values calculated with the 1D setting are also shown by dashed lines in each panel. The blue lines in Figure 1(a) are fitted lines that give the growth rates γ in the linear phase.

Zoom In Zoom Out Reset image size

Figure 1(a) shows that the 3D growth rate ($\gamma /{\omega }_{0}=0.106$) is 26% smaller than the 1D growth rate ($\gamma /{\omega }_{0}=0.144$). This indicates that the 3D structure of density fluctuation works to reduce the growth rate of PDI. Note that both values are smaller than the analytical value ($\gamma /{\omega }_{0}=0.157$) by Goldstein (1978). $\delta {\rho }_{\mathrm{rms}}/{\rho }_{0}$ in the 1D calculation shows an oscillation that is anti-correlated with E⁺ (Figures 1(a), (c)). This is a resonant energy exchange between forward Alfvén and sound waves (Shoda & Yokoyama 2016).

Figure 1(b) shows that the alignment between ${\boldsymbol{j}}$ and ${\boldsymbol{B}}$ occurs because of turbulence. This is sometimes called selective decay (Biskamp 2003); the magnetic field approaches force-free field turbulence. ϕ_j,B cannot be 0 because the system has a finite cross-helicity (Stribling & Matthaeus 1991).

From Figure 1(c), we can tell the energy dissipation rate of Alfvén waves. In the later phase of the 1D calculation, E⁻ is almost constant because there exists no significant physical mechanism that dissipates backward Alfvén waves. Meanwhile, in the nonlinear phase of the 3D calculation, both E⁺ and E⁻ decreases with time, and the larger component has the smaller decay rate. This is the behavior of dynamic alignment (Dobrowolny et al. 1980; Biskamp 2003); the minor Elsässer variable decays faster.

Next, we discuss the structure of turbulence on 2D planes. Figure 2 shows the spatial distributions of the normalized Alfvén speed ${v}_{A}/{v}_{A0}$ (where ${v}_{A}={B}_{z}/\sqrt{4\pi \rho }$) on the xy plane (upper panels) and the normalized amplitude B_y/B_z on the xz plane (lower panels) of different t: from left to right t = 0 (initial condition), $t=12{\tau }_{0}$ (saturation phase), t = 18τ₀ (early nonlinear phase), and t = 24τ₀ (fully nonlinear phase). The inhomogeneity in the upper panel is the trigger of PM, while the filament-like structure is observed in the lower panels because of PM.

As time proceeds, the structure of B_y/B_z becomes aligned with the mean magnetic field ${B}_{0}{{\boldsymbol{e}}}_{z};$ the phase structure becomes progressively finer in the direction perpendicular to the mean field. Both PM and AWT can generate these structures; therefore, we cannot distinguish them here. PM certainly works in the nonlinear phase because there exists a perpendicular gradient of v_A (upper panels in Figure 2); AWT also works because there are bidirectional Alfvén waves and the dynamic alignment is observed (Figures 1(c) (d)). In Section 5, we discuss the dominances of PM and AWT.

Anisotropy also appears in the energy spectrum. To observe the time evolution of anisotropy, we calculate the normalized magnetic energy spectrum E_m with different times. Specifically, we concentrate on the 2D spectrum in the xy directions and xz directions averaged over the other direction, defined as

$$\begin{eqnarray}&&{E}_{m}({k}_{x},{k}_{y})=\displaystyle \frac{{\rm{\Delta }}{k}_{x}{\rm{\Delta }}{k}_{y}}{3L}\int {dz}\displaystyle \frac{{\boldsymbol{B}}{({k}_{x},{k}_{y},z)}^{2}}{{B}_{0}^{2}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{E}_{m}({k}_{x},{k}_{z})=\displaystyle \frac{{\rm{\Delta }}{k}_{x}{\rm{\Delta }}{k}_{z}}{L}\int {dy}\displaystyle \frac{{\boldsymbol{B}}{({k}_{x},y,{k}_{z})}^{2}}{{B}_{0}^{2}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{k}_{x}={\rm{\Delta }}{k}_{y}=\displaystyle \frac{2\pi }{L},\ \ \ {\rm{\Delta }}{k}_{z}=\displaystyle \frac{2\pi }{3L}\end{eqnarray} \tag{ 20 }$$

where, for example, ${\boldsymbol{B}}({k}_{x},y,z)$ denotes the Fourier transformation of ${\boldsymbol{B}}(x,y,z)$ with respect to x.

Figure 3 shows E_m (k_x, k_y) (upper panels) and ${E}_{m}({k}_{x},{k}_{z})$ (lower panels) at t = 12τ₀ (left), t = 15τ₀ (middle), and t = 18τ₀ (right), respectively. The upper panels show isotropic behavior, indicating that the turbulence is axisymmetric with respect to the mean magnetic field. In contrast, the lower panels show anisotropic distributions; the contour is elongated along k_x in the nonlinear phase. This elongation shows that the cascading proceeds faster in the perpendicular (x) direction than in the parallel (z) direction, which is consistent with the classical theory of AWT (Goldreich & Sridhar 1995). In fact, the typical aspect ratio (${k}_{z}/{k}_{x}$) of the spectrum is approximately the same as the nonlinearity ($\delta v/{v}_{A}$), and this indicates the critical balance of AWT (Goldreich & Sridhar 1995).

Zoom In Zoom Out Reset image size

Both PM (phase mixing) and AWT (Alfvén wave turbulence) yield anisotropy, and this similarity makes it difficult to distinguish the physical processes. Here, by calculating the timescale of each process, we aim to clarify which process is dominant in a certain phase.

To estimate the timescale of PM and AWT, we derive the analytical expression of the timescale of PM and AWT. Because both processes are caused by Alfvén waves, reduced MHD (RMHD) approximation is convenient. In the absence of parallel flow and parallel gradients of density and magnetic field, the Alfvén wave propagation is described as follows (Verdini & Velli 2007):

$$\begin{eqnarray}&&[\partial /\partial t\pm {v}_{{\rm{A}}\parallel }{{\rm{\nabla }}}_{\parallel }]{{\boldsymbol{\zeta }}}^{\pm }=-({{\boldsymbol{\zeta }}}^{\mp }\cdot {{\rm{\nabla }}}_{\perp }){{\boldsymbol{\zeta }}}^{\pm },\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\parallel }={{\boldsymbol{e}}}_{z}\displaystyle \frac{\partial }{\partial z},\ \ \ {{\rm{\nabla }}}_{\perp }={{\boldsymbol{e}}}_{x}\displaystyle \frac{\partial }{\partial x}+{{\boldsymbol{e}}}_{y}\displaystyle \frac{\partial }{\partial y}.\end{eqnarray} \tag{ 22 }$$

The second term on the left-hand side is the term for propagation and PM, while the first term on the right-hand side corresponds to AWT. Even though the RMHD approximation does not hold in our simulation, Equation (21) approximately describes the Alfvén wave propagation and provides a general understanding of the physical process that is not related to compressibility.

We begin by deriving the timescale of PM. Nonlinear terms are ignored here because PM is a linear process triggered by the perpendicular variance of the Alfvén speed. Taking the perpendicular rotation (${{\rm{\nabla }}}_{\perp }\times$) of Equation (21), we obtain

$$\begin{eqnarray}&&[\partial /\partial t\pm {v}_{{\rm{A}}\parallel }{{\rm{\nabla }}}_{\parallel }]({{\rm{\nabla }}}_{\perp }\times {{\boldsymbol{\zeta }}}^{\pm })=\mp ({{\rm{\nabla }}}_{\perp }{v}_{{\rm{A}}\parallel })\times \displaystyle \frac{\partial }{\partial z}{{\boldsymbol{\zeta }}}^{\pm }.\end{eqnarray} \tag{ 23 }$$

The right-hand side term represents the phase mixing. From this equation, the growth rate of PM (γ_PM) is approximately given as follows:

$$\begin{eqnarray}&&{\gamma }_{\mathrm{PM}}\approx | {{\rm{\nabla }}}_{\perp }{v}_{{\rm{A}}\parallel }| .\end{eqnarray} \tag{ 24 }$$

Next, we calculate the timescale of AWT. Ignoring the linear term in Equation (21), the decay rate of Alfvén wave turbulence is given as

$$\begin{eqnarray}&&{\gamma }_{\mathrm{AWT}}^{\pm }\approx \displaystyle \frac{| ({\tilde{{\boldsymbol{\zeta }}}}^{\mp }\cdot {{\rm{\nabla }}}_{\perp }){{\boldsymbol{\zeta }}}^{\pm }| }{{\langle {{\boldsymbol{\zeta }}}^{\pm }\rangle }_{\mathrm{rms}}},\end{eqnarray} \tag{ 25 }$$

where ${\tilde{{\boldsymbol{\zeta }}}}^{\mp }$ denote the fluctuating parts of Elsässer variables:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{\zeta }}}}^{\mp }={{\boldsymbol{\zeta }}}^{\mp }-\displaystyle \frac{1}{{L}^{2}}\int {dxdy}{{\boldsymbol{\zeta }}}^{\mp }.\end{eqnarray} \tag{ 26 }$$

We evaluate the nonlinear operator as ${\tilde{{\boldsymbol{\zeta }}}}^{\mp }\cdot {{\rm{\nabla }}}_{\perp }$ instead of ${{\boldsymbol{\zeta }}}^{\mp }\cdot {{\rm{\nabla }}}_{\perp }$ because the perpendicularly uniform mode does not contribute to the wave distortion (energy cascading).

In Figure 4, we show the spatial distribution of ${j}_{z}/{j}_{\mathrm{rms}}$, γ_PM and ${\gamma }_{\mathrm{AWT}}^{\pm }$ on the xy plane in an early nonlinear phase (t = 14.4τ₀, upper panels) and a late nonlinear phase (t = 28.8τ₀, lower panels). Note that ${j}_{z}/{j}_{\mathrm{rms}}$ indicates the degree of development of the perpendicular cascading. In an early phase, γ_PM is larger than ${\gamma }_{\mathrm{AWT}}^{\pm }$ and it is spatially correlated with ${j}_{z}/{j}_{\mathrm{rms}}$ and ${\gamma }_{\mathrm{AWT}}^{\pm }$. This shows that the early nonlinear phase is dominated by PM-driven turbulence (Magyar et al. 2017). In the later phase, the magnitude relation is reversed; ${\gamma }_{\mathrm{AWT}}^{+}$ becomes the largest and this indicates that the later phase is characterized by Alfvén wave turbulence. Specifically, AWT in the later phase is imbalanced (${E}^{+}\ll {E}^{-}$, Figure 1(c)), and thus the dynamic alignment proceeds (Figure 1(d)).

Zoom In Zoom Out Reset image size

Figure 5 shows the rms values of γ_PM (red solid line), ${\gamma }_{\mathrm{AWT}}^{+}$ (blue solid line), and ${\gamma }_{\mathrm{AWT}}^{-}$ (blue dashed line) versus time. In addition to these, we also calculate the normalized growth rate of PDI (γ_PDI) from the dispersion relation given by Goldstein (1978) and show its time evolution with a black line.

Zoom In Zoom Out Reset image size

Figure 5 gives some important indications. First, the system is dominated by different processes depending on the phase. The initial phase ($0\leqslant t/{\tau }_{0}\lesssim 8$) is dominated by the growth of PDI because γ_PDI is the largest. Before the saturation to early nonlinear phase ($8\lesssim t/{\tau }_{0}\lesssim 16$), PM becomes active because of the large density fluctuation (Figure 1(a)). Finally, in the fully nonlinear phase ($16\lesssim t/{\tau }_{0}$), the system is characterized by imbalanced (${\gamma }_{\mathrm{AWT}}^{+}\gg {\gamma }_{\mathrm{AWT}}^{-}$) AWT.

Second, because the timescale of γ_PM is clearly correlated with $\delta {\rho }_{\mathrm{rms}}/{\rho }_{0}$, PM should be of importance in the large-density-fluctuation regions. In the corona and solar wind, such a region is either the coronal bottom (Raymond et al. 2014) or the wind acceleration region (Miyamoto et al. 2014). Specifically, the fast saturation of nonthermal line broadening in the corona (Hahn & Savin 2013) possibly comes from PM because of the presence of large density fluctuations near the coronal bottom (Raymond et al. 2014).

M.S. is supported by the Leading Graduate Course for Frontiers of Mathematical Sciences and Physics (FMSP) and Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows. T.Y. is supported by JSPS KAKENHI grant No. 15H03640. Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.