Many important problems in astrophysics, space physics, and geophysics involve flows of (possibly ionized) gases in the vicinity of a spherical object, such as a star or planet. The geometry of such a system naturally favors numerical schemes based on a spherical mesh. Despite its orthogonality property, the polar (latitude-longitude) mesh is ill suited for computation because of the singularity on the polar axis, leading to a highly non-uniform distribution of zone sizes. The consequences are (a) loss of accuracy due to large variations in zone aspect ratios, and (b) poor computational efficiency from a severe limitations on the time stepping. Geodesic meshes, based on a central projection using a Platonic solid as a template, solve the anisotropy problem, but increase the complexity of the resulting computer code. We describe a new finite volume implementation of Euler and MHD systems of equations on a triangular geodesic mesh (TGM) that is accurate up to fourth order in space and time and conserves the divergence of magnetic field to machine precision. The paper discusses in detail the generation of a TGM, the domain decomposition techniques, three-dimensional conservative reconstruction, and time stepping.

There is a great need in several areas of astrophysics and space-physics to carry out high order of accuracy, divergence-free MHD simulations on spherical meshes. This requires us to pay careful attention to the interplay between mesh quality and numerical algorithms. Methods have been designed that fundamentally integrate high order isoparametric mappings with the other high accuracy algorithms that are needed for divergence-free MHD simulations on geodesic meshes. The goal of this paper is to document such algorithms that are implemented in the geodesic mesh version of the RIEMANN code. The fluid variables are reconstructed using a special kind of WENO-AO algorithm that integrates the mesh geometry into the reconstruction process from the ground-up. A novel divergence-free reconstruction strategy for the magnetic field that performs efficiently at all orders, even on isoparametrically mapped meshes, is then presented. The MHD equations are evolved in space and time using a novel ADER predictor algorithm that is efficiently adapted to the isoparametrically mapped geometry. The application of one-dimensional and multidimensional Riemann solvers at suitable locations on the mesh then provides the corrector step. The corrector step for the magnetic field uses a Yee-type staggering of magnetic fields. This results in a scheme with divergence-free update for the magnetic field. The use of ADER enables a one-step update which only requires one messaging operation per complete timestep. This is very beneficial for parallel processing. Several accuracy tests are presented as are stringent test problems. PetaScale performance is also demonstrated on the largest available supercomputers.

Interstellar Pickup ions (PUIs) play a significant part in mediating the solar wind (SW) interaction with the interstellar medium. In this paper, we examine the details of spatial variation of the PUI velocity distribution function (VDF) in the SW by solving the PUI transport equation. We assume the PUI distribution is isotropic resulting from strong pitch-angle scattering by wave-particle interaction. A three-dimensional model combining the MHD treatment of the background SW and neutrals with a kinetic treatment of PUIs throughout the heliosphere and the surrounding local interstellar medium (LISM) has been developed. The model generates PUI power law tails via second-order Fermi process. We analyze how PUIs transform across the heliospheric termination shock (TS) and obtain the PUI phase space distribution in the inner heliosheath including continuing velocity diffusion. Our simulated PUI spectra are compared with observations made by New Horizons, Ulysses, Voyager 1, 2 and Cassini, and a satisfactory agreement is demonstrated. Some specific features in the observations, for example, a cutoff of PUI VDF at $v = V_{SW}$ and a $f \propto v^{-5}$ tail in the reference frame of the SW, are well represented by the model.

This report describes a new magnetohydrodynamic numerical model based on a hexagonal spherical geodesic grid. The model is designed to simulate astrophysical flows of partially ionized plasmas around a central compact object, such as a star or a planet with a magnetic field. The geodesic grid, produced by a recursive subdivision of a base platonic solid (an icosahedron), is free from control volume singularities inherent in spherical polar grids. Multiple populations of plasma and neutral particles, coupled via charge-exchange interactions, can be simulated simultaneously with this model. Our numerical scheme uses piecewise linear reconstruction on a surface of a sphere in a local two-dimensional "Cartesian" frame. The code employs HLL-type approximate Riemann solvers and includes facilities to control the divergence of magnetic field and maintain pressure positivity. Several test solutions are discussed, including a problem of an interaction between the solar wind and the local interstellar medium, and a simulation of Earth's magnetosphere.

Modeling of the global heliosphere seeks to investigate the interaction of the solar wind with the partially ionized local interstellar medium. Models that treat neutral hydrogen self-consistently and in great detail, together with the plasma, but that neglect magnetic fields, constitute a sub-category within global heliospheric models. There are several different modeling strategies used for this sub-category in the literature. Differences and commonalities in the modeling results from different strategies are pointed out. Plasma-only models and fully self-consistent models from four research groups, for which the neutral species is modeled with either one, three, or four fluids, or else kinetically, are run with the same boundary parameters and equations. They are compared to each other with respect to the locations of key heliospheric boundary locations and with respect to the neutral hydrogen content throughout the heliosphere. In many respects, the models' predictions are similar. In particular, the locations of the termination shock agree to within 7% in the nose direction and to within 14% in the downwind direction. The nose locations of the heliopause agree to within 5%. The filtration of neutral hydrogen from the interstellar medium into the inner heliosphere, however, is model dependent, as are other neutral results including the hydrogen wall. These differences are closely linked to the strength of the interstellar bow shock. The comparison also underlines that it is critical to include neutral hydrogen into global heliospheric models.

We investigate heliosheath energetic neutral atom (ENA) fluxes at keV energies, by assuming that the heliosheath proton distribution can be approximated by a kappa-distribution. The choice of the kappa parameter derives from observational data of the solar wind (SW). This has direct applications to the upcoming IBEX mission. We will look at all-sky ENA maps within the IBEX energy range (10 eV to 6 keV), as well as ENA energy spectra in several directions. We find that the use of kappa, as opposed to a Maxwellian, gives rise to greatly increased ENA fluxes above 1 keV, while medium energy fluxes are somewhat reduced. We show how IBEX data can be used to estimate the spectral slope in the heliosheath, and that the use of kappa reduces the differences between ENA maps at different energies. We also investigate the effect introducing a kappa-distribution has on the global interaction between the SW and the local interstellar medium (LISM), and find that there is generally an increase in energy transport from the heliosphere into the LISM, due to the modified profile of ENA's energies. This results in a termination shock that moves out by 4 AU, a heliopause that moves in by 9 AU and a bow shock 25 AU farther out, in the nose direction.

At present, the heliosphere is embedded in a warm low density interstellar cloud that belongs to a cloud system flowing through the local standard of rest with a velocity near ~18 km/s. The velocity structure of the nearest interstellar material (ISM), combined with theoretical models of the local interstellar cloud (LIC), suggest that the Sun passes through cloudlets on timescales of < 10^3 - 10^4 yr, so the heliosphere has been, and will be, exposed to different interstellar environments over time. By means of a multi-fluid model that treats plasma and neutral hydrogen self-consistently, the interaction of the solar wind with a variety of partially ionized ISM is investigated, with the focus on low density cloudlets such as are currently near the Sun. Under the assumption that the basic solar wind parameters remain/were as they are today, a range of ISM parameters (from cold neutral to hot ionized, with various densities and velocities) is considered. In response to different interstellar boundary conditions, the heliospheric size and structure change, as does the abundance of interstellar and secondary neutrals in the inner heliosphere, and the cosmic ray level in the vicinity of Earth. Some empirical relations between interstellar parameters and heliospheric boundary locations, as well as neutral densities, are extracted from the models.

We present initial results from our new two-dimensional (radius and latitude), self-consistent model of galactic cosmic rays in the heliosphere. We focus on the latitudinal variations in the solar wind flow caused by the energetic particles. Among other things our results show that the cosmic rays significantly modify the latitudinal structure of the solar wind flow downstream of the termination shock. Specifically, for A>0 (corresponding to the present solar minimum) the wind beyond the shock is driven towards the equator, resulting in a faster wind flow near the current sheet, while for A<0 the effect is reversed and the wind turns towards the pole, with a faster flow at high latitudes. We attribute this effect to the latitudinal gradients in the cosmic ray pressure, caused by drifts, that squeeze the flow towards the ecliptic plane or the pole, respectively.