Nonlinear Evolution of Ion Kinetic Instabilities in the Solar Wind

TL;DR

This study uses a novel 3D hybrid model to investigate the nonlinear evolution of ion kinetic instabilities in the solar wind, revealing their role in heating and magnetic fluctuation dissipation near the Sun.

Contribution

First application of a 3D hybrid PIC-fluid model to simulate ion kinetic instabilities in the solar wind, capturing nonlinear evolution and dissipation processes.

Findings

Self-consistent non-Maxwellian ion VDFs produced

Temperature anisotropies observed in unstable ion populations

Wave spectra consistent with solar wind observations

Abstract

In-situ observations of the solar wind (SW) plasma from 0.29 to 1AU show that the protons and alpha particles are often non-Maxwellian, with evidence of kinetic instabilities, temperature anisotropies, differential ion streaming, and associated magnetic fluctuations spectra. The kinetic instabilities in the SW multi-ion plasma can lead to preferential heating of alpha particles and the dissipation of magnetic fluctuation energy, affecting the kinetic and global properties of the SW. Using for the first time a three-dimensional hybrid model, where ions are modeled as particle using the Particle-In-Cell (PIC) method and electrons are treated as fluid, we study the onset, nonlinear evolution and dissipation of ion kinetic instabilities. The Alfven/ion-cyclotron, and the ion drift instabilities are modeled in the region close to the Sun (~10R_s). Solar wind expansion is incorporated in the model. The model produces self-consistent non-Maxwellian velocity distribution functions (VDFs) of unstable ion populations, the associated temperature anisotropies, and wave spectra for several typical SW instability cases in the nonlinear growth and saturation stage of the instabilities. The 3D hybrid modeling of the multi-ion SW plasma could be used to study the SW acceleration region close to the Sun that will be explored by the Parker Solar Probe mission.

Tables

Table 1. Table 1. : The parameters of the initial states for the 3D hybrid modeling cases in the present study. In all cases the proton population is initially anisotropic, β ∥ = β ∥ , p = β ∥ , α = β e \beta_{\parallel}=\beta_{\parallel,p}=\beta_{\parallel,\alpha}=\beta_{e} , and n α = 0.05 ​ n e subscript 𝑛 𝛼 0.05 subscript 𝑛 𝑒 n_{\alpha}=0.05n_{e} .

Case	Instability type	${\beta }_{\parallel }$	$T_{⟂,\alpha }/T_{\parallel ,\alpha }$	$V_d\left[V_A\right]$	$ϵ$
1	Drift	0.04	1	2	0
2	Ion-Cyclotron	0.04	10	0	0
3	Ion-Cyclotron	1	10	0	0
4	Ion-Cyclotron	1	10	0	${10}^{-3}$

Equations

$$\frac{d{\text{mathbfit}x}_s}{dt}={\text{mathbfit}v}_s$$

$$m_s\frac{d{\text{mathbfit}v}_s}{dt}=Z_{s}e\left(\text{mathbfit}E+\frac{{\text{mathbfit}v}_s\times \text{mathbfit}B}{c}\right),$$

$$e{n}_e\left(\text{mathbfit}E+\frac{{\text{mathbfit}v}_e\times \text{mathbfit}B}{c}\right)+\text{mathbfit}\nabla p_e=0,$$

$$R(t)=R_0+U_{0}t,$$

$$a(t)=R(t)/R_0=1+\frac{U_0}{R_0}t.$$

$$x=x^{\prime }+R,y=a{y}^{\prime },z=a{z}^{\prime }.$$

$$\displaystyle\frac{d\mbox{\bf v}}{dt}=\frac{d\mbox{\bf v}^{\prime}}{dt^{\prime}}+\frac{\dot{a}}{a}\mbox{\bf P}\cdot\mbox{\bf v}^{\prime},\ \mbox{\bf P}=\left(\begin{array}[]{ccc}0&0&0\\ 0&1&0\\ 0&0&1\end{array}\right),$$

$$\frac{d{x}^{\prime }}{dt}=v_x^{\prime },\frac{d{y}^{\prime }}{dt}=\frac{1}{a}{v}_y^{\prime },\frac{d{y}^{\prime }}{dt}=\frac{1}{a}{v}_y^{\prime }.$$

$$\frac{\partial \text{mbox}{\mathbf{B}}^{\prime }}{\partial t^{\prime }}+\text{mbox}{\mathbf{B}}^{\prime }({\nabla }^{\prime }\cdot \text{mbox}{\mathbf{U}}_i^{\prime })-(\text{mbox}{\mathbf{B}}^{\prime }\cdot {\nabla }^{\prime })\text{mbox}{\mathbf{U}}_i^{\prime }+(\text{mbox}{\mathbf{U}}_i^{\prime }\cdot {\nabla }^{\prime })\text{mbox}{\mathbf{B}}^{\prime }$$

$$+\frac{c}{4\pi e{n}_e}{\nabla }^{\prime }\times \left[({\nabla }^{\prime }\times \text{mbox}{\mathbf{B}}^{\prime })\times \text{mbox}{\mathbf{B}}^{\prime }\right]=-\frac{2U_0}{R_0}{B}_x^{\prime }\hat{x}-\frac{U_0}{R_0}{B}_y^{\prime }\hat{y}-\frac{U_0}{R_0}{B}_z^{\prime }\hat{z},$$

$$\text{mbox}{\mathbf{E}}^{\prime }=$$

Full text

\setlastpage\inarticletrue{opening}

1 Introduction

\ilabel

int:sec

In-situ observations of the solar wind (SW) plasma at distances $>$0.29AU by Helios, Ulysses, ACE, and Wind spacecraft show that the $\alpha$ particle population is usually hotter than that of the proton population, and flows faster by an Alfvén speed in the fast SW streams. The VDFs of protons and ions exhibit non-Maxwellian features, such as temperature anisotropy with respect to the background magnetic field, beams (with stronger departures from Maxwellian for the $\alpha$ particles), and differential ion streaming (Marsch et al., 1982). These effects are stronger in the fast SW streams (compared to slow SW), increasing in magnitude closer to the Sun. Kinetic instabilities, such as ion-cyclotron, mirror and firehose, play an important role in shaping the SW plasma properties, as evident from observations at 1AU with WIND (Bale et al., 2009; Maruca, Kasper, and Gary, 2012) and other spacecraft data.

Observed break points in the magnetic fluctuations power spectra indicate the inertial and kinetic dissipation ranges of magnetic field fluctuations (see the review, Bruno and Carbone, 2013), with the break point of the dissipation range aligned with the proton cyclotron frequencies at the various heliocentric distances (e.g., Bourouaine et al., 2012; Telloni, Bruno, and Trenchi, 2015). Direct evidence of electromagnetic ion cyclotron (EMIC) waves near the proton cyclotron frequency was found in the SW by STEREO at 1AU (e.g., Jian et al., 2009, 2014) and from Helios and MESSENGER at 0.3 AU (Jian et al., 2010). Preferential acceleration and heating of the $\alpha$ particle populations was associated with wave activity in SW observations at 1AU (Kasper, Lazarus, and Gary, 2008; Kasper et al., 2013; Bourouaine, Marsch, and Neubauer, 2011b, a; Bourouaine et al., 2013; Maruca, Kasper, and Bale, 2011; Maruca, Kasper, and Gary, 2012; Borovsky and Gary, 2014; Borovsky, 2016; Durovcová, Němeček, and Šafránková, 2019).

Recently launched NASA’s Parker Solar Probe (PSP) mission, is going to provide detailed in-situ observations as close as 10$R_{s}$, aimed at studying the SW acceleration and heating processes (Fox et al., 2016). In anticipation of the PSP SW ions observations (Kasper et al., 2016) we developed a 3D hybrid model of SW plasma with protons and $\alpha$ particles and study the onset, evolution, and dissipation of ion driven kinetic instabilities, using our 3D hybrid kinetic, expanding computational box model (see model details below).

Recently, ion kinetic instabilities in SW electron-proton-alpha ($e-p-\alpha$) plasma were studied extensively with 2.5D hybrid models, including effects such as SW expansion and background inhomogeneity, and turbulence (e.g., Ofman, Viñas, and Maneva, 2014; Ozak, Ofman, and Viñas, 2015; Maneva, Ofman, and Viñas, 2015; Ofman, Viñas, and Roberts, 2017; Maneva and Poedts, 2018). The advantage of the hybrid models over other kinetic modeling methods, such as full PIC and Vlasov solvers is that they can model the nonlinear evolution of the ion kinetic instabilities for parameter ranges, physical size, and duration not accessible by other methods, relaxing the approximations on mass ratio and speed of light used in PIC methods (the trade-off is the limitation of hybrid modeling to plasmas where electron kinetic evolution or instabilities could be neglected or are uncoupled from the ions), and the severe resolution limitations in velocity space of multi-dimensional Vlasov codes (e.g., Perrone et al., 2014). Three-dimensional hybrid codes were developed for the study of SW $e-p$ plasma turbulence with more general (than 2.5D) description of the wave-particle interactions (e.g., Vasquez, Markovskii, and Chandran, 2014; Vasquez, 2015; Franci et al., 2018). Here, for the first time, we use 3D hybrid model of $e-p-\alpha$ SW plasma to demonstrate the anisotropic ion heating as a result of ion drift and temperature anisotropy instabilities, and the associated high frequency Alfvén/ion-cyclotron wave spectra in the SW plasma using parameters relevant to the heliocentric distances that will be sampled for the first time by the recently launched PSP mission. The paper is organized as follows: in Section \irefhyb:sec we describe the details of the hybrid model, in Section \irefnum:sec we present the numerical results, and the summary and conclusions are in Section \irefcon:sec.

2 Hybrid Model

\ilabel

hyb:sec We use the 3D hybrid model where the protons and $\alpha$ particles are described kinetically as particles using the Particle-In-Cell (PIC) method, while electrons are treated as neutralizing massless background fluid. The parallelized 3D hybrid code is an extension of the 1.5D hybrid code initially developed by Winske and Omidi (1993), and parallelized 2.5D hybrid models developed to study multi-ion SW plasma (e.g., Ofman and Viñas, 2007; Ofman, Viñas, and Maneva, 2014; Ofman, Viñas, and Roberts, 2017). The equations of motion (with the usual notations for the variables) are solved for each particle of all ion species ($s$) subject to the Lorentz force for 3D position x and velocity v vectors:

[TABLE]

where $m_{s}$ is the particles mass, $Z_{s}$ is the charge number, $e$ is the electron charge, $c$ is the speed of light, ${\mathbfit E}$ is the electric field, and ${\mathbfit B}$ is the magnetic field. The electron momentum equation is solved in 3D by using negligible electron inertia leading to the generalized Ohm’s law

[TABLE]

where the scalar electron pressure equation of state $p_{e}=k_{B}n_{e}T_{e}$ is used for closure, where $T_{e}$ is the electron temperature, and $n_{e}$ is the electron number density. For stability a resistive term $en_{e}\eta{\mathbfit J}$ is added to the r.h.s. of Equation \irefohm:eq, where $\eta$ is an empirical small resistivity coefficient, and ${\mathbfit J}$ is the local current density. The above equations are supplemented by the quasi-neutrality condition $n_{e}=n_{p}+2n_{\alpha}$, where $n_{p}$ is the proton population and $n_{\alpha}$ is the $\alpha$ particle population number density, respectively. Maxwell’s equations $\nabla\times{\mathbfit B}=\frac{4\pi}{c}{\mathbfit J},$ and $\nabla\times{\mathbfit E}=-\frac{1}{c}\frac{\partial{\mathbfit B}}{\partial t}$ are solved on the 3D spatial grid, with grid sizes of up to $128^{3}$ and 64 to 128 particles per cell (ppc) in the present study. The proton and $\alpha$ particle equations of motions are advanced in time as the particle motions respond to the updated Lorentz force at each time step. The particle and field equations are integrated in time using the Rational Runge-Kutta (RRK) method (Wambecq, 1978) whereas the spatial derivatives are calculated by pseudospectral FFT method in 3D, and periodic boundary conditions are applied. The equations of motions of the ions are solved using Cartesian coordinates with the background uniform magnetic field $B_{0}$ direction defined along ${\rm\bf x}$. The spatial coordinates are normalized using the proton inertia length $\delta=c/\omega_{pp}$, where $\omega_{pp}=(4\pi n_{p}e^{2}/m_{p})^{1/2}$ is the proton plasma frequency, the time is in units of the inverse proton cyclotron frequency, $\Omega_{p}^{-1}$, where $\Omega_{p}=\frac{eB_{0}}{m_{p}c}$, and the velocities are in units of the Alfvén speed, $V_{A}$. In the present study the spatial resolution is up to $0.75\delta$, and the time step is on the order of $0.01\Omega_{p}^{-1}$. The hybrid modeling method has been tested and used successfully in many previous studies, and details of the normalization were published (see, e.g., Winske and Omidi, 1993; Ofman and Viñas, 2007; Ofman, 2010; Ofman, Viñas, and Maneva, 2014; Ofman, Viñas, and Roberts, 2017).

2.1 Expanding Box Model

A solar wind plasma parcel expands naturally as it travels away from the Sun into the heliosphere at the SW speed, $U_{0}$. In order to model the expansion of the SW, we use the Expanding Box model (EBM), developed initially by Grappin and Velli (1996) for an MHD fluid code and later adapted by Liewer, Velli, and Goldstein (2001) in their 1.5D hybrid code, (see, also Hellinger et al., 2003, 2005). We have implemented the EBM in our 2.5D hybrid code (Ofman, Viñas, and Moya, 2011), and recently applied for the first time to our 3D hybrid code. Here, we summarize the equations given in Ofman, Viñas, and Moya (2011), extended to the 3D hybrid model.

The heliocentric radial position, $R(t)$, of the SW plasma parcel as a function of time is approximated by

[TABLE]

where a constant SW velocity $U_{0}$ is assumed (i.e., applicable beyond $\sim 10R_{s}$), and $R_{0}$ is the reference radial distance. The dimensionless expansion factor $a(t)$ is

[TABLE]

The Galilean transformation $v^{\prime}=v-U_{0}$ is used to transform to the SW plasma rest frame. Since the expansion rate of the SW is ’slow’ the expansion parameter is small, $\epsilon=\frac{U_{0}}{R_{0}}t\ll 1$, where time is in units of $\Omega_{p}^{-1}$. Thus, the coordinates are transformed to the moving frame as follows:

[TABLE]

While the $x$ coordinate (along the radial direction) undergoes Galilean transformation, the $y$ and $z$ coordinates undergo stretching, or inflation. Thus, the ion equation of motions in the 3D hybrid model transforms as

[TABLE]

where we note that for $\epsilon\ll 1$, we can use the approximation $\frac{\dot{a}}{a}\sim\frac{U_{0}}{R_{0}}$, and where v′** components are given by

[TABLE]

For the magnetic field we have

[TABLE]

where the transformed ion bulk velocity is $\mbox{\bf U}_{i}^{\prime}=\mbox{\bf U}_{i}-\mbox{\bf U}_{0}$, and the electric field in the moving frame are

[TABLE]

The above transformations are implemented in the 3D hybrid code, including the derivatives with respect to $x^{\prime}$, $y^{\prime}$ and $z^{\prime}$, with the expansion parameter $\epsilon$. The typical expansion parameter, $\epsilon$ at $10R_{s}$ is $10^{-4}-10^{-5}$, while setting $\epsilon=0$, recovers the non-expanding model. In order to evaluate the effects of expansion, we present preliminary results with a higher expansion rate, due to computational limitations of the 3D hybrid model, while more realistic (smaller) values of $\epsilon$ requiring longer runs left for a future study.

3 Numerical Results

\ilabel

num:sec In Figures \iref3Dvs2D:fig-\irefk_spect:fig we present the results of the 3D hybrid modeling for the cases with parameters summarized in Table 1. The parameters of the instabilities in this study exceed the linear stability threshold in all cases as determined form solutions of Vlasov’s linear dispersion relation (see, e.g., Davidson and Ogden, 1975; Gary, 1993; Gary et al., 2001, 2003; Xie, Ofman, and Viñas, 2004; Ofman and Viñas, 2007). We have used 5% $\alpha$ particles number density (in terms of $n_{e}$). In Figure \iref3Dvs2D:fig we compare the evolution of the ion drift instability (Case 1), modeled with both, 2.5D and 3D hybrid codes with identical parameters. It is evident that the relaxation of the drift instability is rapid in terms of proton gyro-periods, and the drift velocity relaxes from 2 to $\sim$1.6 in 200$\Omega_{p}^{-1}$. At the same time, the ion populations are heated in the perpendicular direction and becomes anisotropic with maximal $T_{\perp,p}/T_{\parallel,p}\approx 2.1$, and $T_{\perp,\alpha}/T_{\parallel,\alpha}\approx 5.6$. The 2.5D and 3D hybrid models produce similar results, with the initial growth somewhat slower an higher final proton anisotropy in the 3D hybrid model. The final states of $\alpha$ population temperature anisotropy and drift velocity are close in both models. The similarity in the evolution is expected, since most of the power in the modes are in the parallel (to the magnetic field) direction. These results provide additional validation of our newly develop 3D hybrid model.

In Figure \irefVDF_drift:fig we show the velocity distribution functions (VDFs) of the protons and $\alpha$ particle populations for Case 1 at $t=90\Omega_{p}^{-1}$. The anisotropy and the non-Maxwellian features are evident in the $V_{x}-V_{z}$ velocity space of both ion species VDFs. The larger (than proton) anisotropy is seen in the $\alpha$ population VDF. In the lower panel of Figure \irefVDF_drift:fig we show cuts of the VDFs at $V_{z}=0$, obtained by integrating the number of particles in velocity bins (i.e., $\Delta V_{x}\ll 1$) for each species, and the best-fit Maxwellians. It is evident that protons are nearly Maxwellian in the parallel direction, while $\alpha$ particle population has a Maxwellian core with a non-thermal tail, likely produced by the nonlinearity of the drift instability, i.e., wave-particle interactions and pitch-angle scattering populating the VDF tails. The 3D hybrid modeling results are similar to the previous 2.5D hybrid modeling results of this instability (see, e.g., Ofman, Viñas, and Maneva, 2014).

The temporal evolution of the perpendicular magnetic fluctuations power, $|B_{\perp}|^{2}$, are shown in Figure \irefbperp_drift:fig. The perpendicular fluctuations are primarily due to the wave power as demonstrated in previous studies (for example, see the dispersion relation in Figure 12 in Ofman, Viñas, and Roberts (2017) obtained from 2.5D hybrid model). The initial rapid growth of the magnetic power accompanies the initial growth of the drift instability evident in Figure \iref3Dvs2D:fig. The peak perpendicular magnetic power is reached at $t\approx 55\Omega_{p}^{-1}$, followed by nonlinear saturation and gradual dissipation. This is consistent with the perpendicular heating of the ions, as evident in the evolution of their anisotropies in Figure \iref3Dvs2D:fig, suggesting the Alfvén/cyclotron wave absorption as the heating process.

The evolution of the temperature anisotropies, and parallel and perpendicular kinetic energies of protons and $\alpha$ particle populations for Case 2 are shown in Figure \irefAniso10:fig. The ion-cyclotron instability is driven by the initial temperature anisotropy of the $\alpha$ particles, $T_{\perp,\alpha}/T_{\parallel,\alpha}$=10, while protons are initially isotropic. It is evident that the temperature anisotropy of the $\alpha$ population is gradually decreasing with time, while protons remain nearly isotropic (note the small $y$-axis range of the proton anisotropy plot).

The decrease of $\alpha$ temperature anisotropy is due to the decrease of their perpendicular kinetic energy, accompanied by an increase of the parallel energy due to velocity phase space diffusion. At the same time, the proton parallel and perpendicular kinetic energies remain nearly constant. At $t=600\Omega_{p}^{-1}$ the $\alpha$ particle population temperature anisotropy has decreased by a factor of two to $\sim 5$, approaching gradually the equilibrium state. Due to computational limitation, the run was not continued for a longer duration.

The temporal evolution of the temperature anisotropy of proton and $\alpha$ particle populations for Cases 3 and 4 (Table 1) are shown in Figure \irefexpansion:fig. The initial $T_{\perp,\alpha}/T_{\parallel,\alpha}=10$ and initial $\beta_{\parallel,\alpha}=\beta_{\parallel,p}=\beta_{e}=1$ (where we use $n_{e}$ in the definition of $\beta$ of all species). As expected from linear theory, the relaxation of the ion-cyclotron instability driven by the $\alpha$ temperature anisotropy is faster in this case compared to the low-$\beta_{\parallel}$ case shown in Figure \irefAniso10:fig. The relaxation of the $\alpha$ temperature anisotropy occurs in $\sim 100\Omega_{p}^{-1}$ and an asymptotic state with $T_{\perp,\alpha}/T_{\parallel,\alpha}\approx 2$ is reached and remains to the end of the run. The protons are initially isotropic and undergo perpendicular heating due to the absorptions of the wave-spectrum produced by the relaxation of the $\alpha$ population. The maximal proton anisotropy is $\sim 1.17$ for the non-expanding case, relaxing towards isotropic state. When gradual expansion is introduced in the model ($\epsilon=10^{-3}$), it is evident that the relaxation of the temperature anisotropy is more rapid than in the non-expanding case, with further perpendicular cooling of the protons. This is expected due to the ‘stretching’ of the perpendicular coordinates as the parcel of plasma expands in the heliosphere, and consistent with previous 2.5D hybrid EBM modeling results (e.g., Ofman, Viñas, and Moya, 2011; Ofman, Viñas, and Maneva, 2014).

The VDFs of the proton and $\alpha$ particle populations for the ion-cyclotron instability (Case 3) are shown in Figure \irefVDF_IC:fig at the end of the evolution at $t=200\Omega_{p}^{-1}$ near the asymptotic quasi-steady state. The anisotropy and the non-Maxwellian features are evident in both proton and $\alpha$ particle VDFs in the $V_{x}-V_{z}$ phase space, with the larger anisotropy of $\alpha$ particles. The perpendicular velocity phase-space plane is nearly Maxwellian in both species, as evident in the circular shapes of the VDFs. This property is due to the dominance of the parallel propagating modes, with little effect on the perpendicular direction (see, Figure \irefk_spect:fig, below). The results are consistent with previous 2.5D modeling studies of this instability (e.g., Ofman, Viñas, and Moya, 2011).

The 2D power spectra of the magnetic fluctuations for the instabilities (Case 1 and Case 3) obtained from the 3D hybrid model are shown in Figure \irefk_spect:fig. It is evident that for both cases the peak power of the magnetic fluctuations is located near the $k_{x}$ axis, while the power at $|k_{\perp}|>0$ is significantly diminished. There is evidence for somewhat larger power in the oblique modes for the drift instability (Case 1), compared to the ion-cyclotron instability (Case 3). This result is consistent with the dominant growth of the parallel propagating modes expected from solutions of Vlasov’s linear dispersion relation (e.g., Ofman and Viñas, 2007). We have found that as the instabilities evolve and dissipate, the power at higher $k$ is dissipated, and the peak power in $k$-space is moving towards the origin in $k$-space (i.e., longer wavelengths).

4 Summary and Conclusions

\ilabel

con:sec

Motivated by past SW ion observations at 1AU and the inner heliosphere and by anticipated data from the PSP mission close to the Sun, we investigate kinetic instabilities in $e-p-\alpha$ SW plasma. The instabilities are initialized with temperature anisotropy ($T_{\perp}>T_{\parallel}$) of $\alpha$ particle populations, and an initial relative proton-$\alpha$ super-Alfvénic drift modeled for the first time with full 3D hybrid simulations and EBM implementation. We extend previous 2.5D hybrid modeling studies to more realistic 3D model and find general agreement with previous 2.5D result, with 3D effects affecting properties, such as the growth or relaxation rate of the instabilities in the nonlinear stage.

The instability leads to rapid relaxation of the super-Alfvénic drift, producing EMIC waves and associated perpendicular magnetic fluctuations power. The perpendicular velocity distributions become non-Maxwellian and $\alpha$ particle population is heated with strong perpendicular temperature anisotropy similar to previous 2.5D studies and spacecraft observations. We find that the SW protons are heated significantly in the perpendicular direction by the ion drift instability and the associated power of kinetic waves. We calculate the $k$-space power spectra of the waves and find that they are dominated by parallel propagating waves, with a small oblique component.

We demonstrate the various forms of self-consistent non-Maxwellian VDFs of protons and $\alpha$ particle populations produced by the modeled instabilities and wave-particle interactions with the associated wave spectra in the nonlinear and relaxed stages of the evolution. The results of the study are relevant to the understanding of the kinetic wave-particle processes that likely take place in the acceleration region of the SW multi-ion plasma close to the Sun ($\sim 10R_{s}$), the target investigation region of the recently launched NASA’s PSP mission, where the study of heating and acceleration of proton and $\alpha$ populations is one of the major goals.

Acknowledgments

The author acknowledges support by NASA cooperative agreement NNG11PL10A 670.154 to Catholic University of America. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

Disclosure of Potential Conflicts of Interest The author declares he has no conflicts of interest.