The interaction between the solar wind and solar system bodies, such as planets, satellites, and asteroids, is one of the fundamental global-scale phenomena in space plasma physics. In the present study, the electromagnetic environment around a small dielectric body with a weak intrinsic magnetic field is studied by means of a first-principle kinetic plasma simulation, which is a challenging task in space plasma physics as well as high-performance computing. Due to several computational limitations, five-dimensional full electromagnetic Vlasov simulations with two configuration space and three velocity space coordinates are performed with two different spatial resolutions. The Debye-scale charge separation is not solved correctly in the simulation run with a low spatial resolution, while all the physical processes in collisionless plasma are included in the simulation run with a high spatial resolution. The direction comparison of electromagnetic fields between the two runs shows that there is small difference in the structure of magnetic field lines. On the other hand, small-scale fine structures of electrostatic fields are enhanced by the electric charge separation and the charge accumulation on the surface of the body in the high-resolution run, while these structures are absent in the low-resolution runs. These results are consistent with the conventional understanding of plasma physics that the structure and dynamics of global magnetic fields, which are generally described by the magneto-hydro-dynamics (MHD) equations, are not affected by electron-scale microphysics.

Background

The Earth is surrounded by the collisionless space plasma and the intrinsic magnetic field of the Earth (magnetosphere). The plasma is regarded as the ‘fourth state’ of the matter in which a certain portion of the particles are ionized to ions with positive charge and electrons with negative charge. The space plasma is a highly nonlinear and multi-scale medium where charged particles and electromagnetic fields are strongly coupled with each other without collision among plasma particles. The electromagnetic environment around a solar system body is formed by the interaction between a space plasma flow from the sun (the solar wind) and the magnetosphere of the body. The structure and dynamics of the electromagnetic environment are quite different depending on the strength of the intrinsic magnetic field and the condition of the surface of the bodies. The intrinsic magnetic field generally acts as an obstacle to the solar wind (e.g., for the case of Mercury, Earth, Jupiter, and Saturn). When the intrinsic magnetic field is weak, ionized gases from the upper atmosphere of a body act as a conducting obstacle (e.g., for the case of Venus and Mars). When a body has neither intrinsic magnetic field nor atmosphere, on the other hand, the solar wind plasma directly impacts the body. The interaction between the solar wind and an unmagnetized body is one of the fundamental issues in space plasma physics (Halekas et al. 2012). The Earth’s satellite Moon is one of the examples of a dielectric/insulating body without intrinsic magnetic field and atmosphere.

The computer simulation is an essential approach in studies of highly nonlinear phenomena involving space plasma. There are numerous types of self-consistent plasma simulations which treat plasma according to various approximations. The interaction between the solar wind and a solar system body is regarded as a ‘global-scale’ phenomenon in which magneto-hydro-dynamics (MHD) is dominant and is conventionally studied by means of three-dimensional (3D) MHD simulations where all of the plasma particles are treated by the single-fluid approximation. On the other hand, one of the major purposes of the present study is to extend the conventional single-fluid simulation to a first-principle kinetic simulation, which is a challenging task in space plasma physics as well as high-performance computing. Table 1 shows typical examples and characteristics for global simulations of solar system bodies with various simulation techniques. The MHD equations can handle solar system bodies of any size since there is no characteristic spatial scale in the MHD equations. However, the charge separation and the charge accumulation onto the surface of the bodies cannot be handled with the single-fluid approximation. The charge separation and the charge accumulation can be handled by treating both ion motion and electron motion individually. Since the difference between ion dynamics and electron dynamics results in the ion inertial effect (or well-known Hall effect), the typical spatial scale of solar system bodies in the ion-electron multi-fluids is of several tens of the ion inertial length (di). The fluid approximations, however, cannot handle single-particle effects such as gyro motion. In the ‘hybrid’ particle-in-cell (PIC) simulation where ions are treated as individual charged particles while electrons are treated as a charge-neutralizing massless fluid, the typical spatial scale of solar system bodies is several tens of the gyro radius of thermal ions (ri), which is a realistic scale for the Moon (e.g., Travnicek et al. 2005; Kallio 2005; Holmstrom et al. 2010, 2012; Kallio et al. 2012). However, the charge separation and the charge accumulation onto the surface of the bodies cannot be handled in the hybrid PIC simulation, which strongly affect the electromagnetic environment around the body (Kimura and Nakagawa 2008).

Table 1

Typical example for global simulations of solar system bodies with various simulation techniques

The first-principle full kinetic (PIC and Vlasov) simulations are unique ways to handle all the important physics of the interaction between the solar wind and a dielectric/insulating body. The full PIC simulation treats both ions and electrons as individual charged particles. The typical spatial scale of solar system bodies in the full PIC simulation is of several tens of the electron Debye length (λDe) which is the minimum spatial scale in space plasma (e.g., Guio and Pecseli 2004; Kimura and Nakagawa 2008; Nakagawa and Kimura 2011; Nakagawa 2013). It is rather difficult to handle solar system bodies with larger spatial scales in the full PIC simulation. This is because the number density of plasma particles is highly inhomogeneous in the global magnetosphere and it is not easy to balance computational load among computational domains (or computing processes) in massively parallel PIC simulation. The full Vlasov simulation used in the present study is an alternative to the full PIC simulation. The Vlasov simulation follows the temporal development of distribution functions defined in the position-velocity phase space based on the Vlasov (collisionless Boltzmann) equation. The full Vlasov method is free from numerical noise and allows us to relax the constraint on the grid spacing to be shorter than the Debye length. The Vlasov simulation can use a grid spacing longer than the Debye length as long as the E × B drift motion of thermal plasma is resolved appropriately (Umeda et al. 2009, 2010). Thus, the Vlasov simulation can handle solar system bodies with the spatial scale of several tens of the gyro radius of thermal electrons (re) (Umeda et al. 2011; Umeda 2012a; Umeda and Ito 2014). One critical drawback of the Vlasov simulation is that huge computational resources are required. Hence, we use massive-parallel computers with the help of high-performance computing technology.

Recently, a global Vlasov simulation of a small dielectric body with a small magnetosphere has been carried out (Umeda and Ito 2014). In this previous study, the spatial grid size was set to be ∆x = ∆y ≫ λDe due to the limitation of computational resources (hereafter, this previous simulation run is referred to as the ‘low-resolution run’). In the low-resolution run, the Debye shielding, which is an important physical process on the minimum spatial scale, was not included. In this paper, we present a preliminary result from a ‘high-resolution run’ with ∆x = ∆y ~ λDe which has been obtained by using the K computer. The purpose of this paper is to make a direct comparison between the previous low-resolution run and the present high-resolution run for discussing the effect of electron-scale kinetics to MHD-scale fluid dynamics.

where E, B, J, ρ, μ0, ∈0, and c represent electric field, magnetic field, current density, charge density, magnetic permeability, dielectric constant, and light speed, respectively. The Vlasov equation (2) describes the development of the distribution functions by the electromagnetic (Lorentz) force, with the collision term in the right-hand side set to zero. The distribution function f(t, r, v) is defined as a function of time t, position r, and velocity v with the subscript s being the species of singly charged particles (e.g., s = i and e for ions and electrons, respectively). The quantities qs and ms are the charge and the mass of the particle species s. The motion of charged particles results in the electric current. By integrating the Vlasov equation (2) over the velocity v and getting the sum for all the particle species s, the continuity equation for charge is obtained, which describes the necessary condition for the electric current density,

$$ \frac{\partial \rho }{\partial t}+\nabla \cdot \mathbf{J}=0 $$

(3)

These equations are regarded as the ‘first principle’ of the collisionless plasma. Note that the fluid equations are derived by taking zeroth moment (mass conservation), the first moment (momentum conservation), and the second moment (energy conservation) of the Vlasov equation (2).

Our Vlasov code solves two spatial dimensions (x, y) and three velocity dimensions (vx, vy, vz), which is so-called two-and-half dimension (2.5D). The numerical schemes are based on a conservative semi-Lagrangian scheme with several recent improvements (Schmitz and Grauer 2006; Umeda 2008; Umeda et al. 2009, 2012a) and hybrid parallelization with MPI library and OpenMP (Umeda et al. 2012b). The detail descriptions on the numerical schemes are found in the literature (Umeda 2012b) and are not repeated here.

Simulation setup

Even with the fastest supercomputer in the world in 2014, it is impossible to use realistic physical parameters for the Moon in the first-principle full kinetic simulation. Hence, we need to ‘reduce’ the physical parameters to enable a numerical simulation with a limited computer resource. The detailed simulation parameters and setting are described in the literature (Umeda and Ito 2014) and are briefly reviewed here.

There exists a uniform plasma flow, the solar wind, in the simulation box at the initial state. The plasma flow carries a background magnetic field as the interplanetary magnetic field (IMF), which is the intrinsic magnetic field of the Sun. The physical parameters of the solar wind are as follows: the ion bulk velocity Us= 400 km/s, the ion thermal velocity Vti = 80 km/s, and the (ion and electron) number density N0 = 5 cm−2, which are typical values for the solar wind ions. It is also assumed that the ions and electrons have the same temperature (Ti/Te = 1). In order to handle the ion gyro motion in the full kinetic simulation, on the other hand, the magnitude of the background magnetic field is higher than the realistic value and is set to B0 = 614 nT. Thus, the ion plasma and gyro frequencies in the solar wind are fpi ~ 469 Hz and fci ~ 936 Hz, respectively. The (ion and electron) Debye length and the ion thermal gyro radius are λDi = λDe ~ 27.2 m and ri ~ 1.36 km, respectively. We also need to reduce the ion-to-electron mass ratio to mi/me = 25 for computational cost. That is, electrons are assumed to be about 70 times heavier than the real ones. Then, the electron thermal gyro radius is re ~ 272 m with Ti/Te = 1.

There is an inner boundary at \( \left|\mathbf{r}\right|=\sqrt{x^2+{y}^2}={R}_L \) treated as a dielectric circular (cylindrical) body, where RL denotes the radius of the body. We solve the Vlasov equation (2) only for |r| > RL while we solve the Maxwell equations (1) in the entire region. With these conditions, charged particles are accumulated on the surface of the body while electromagnetic waves propagate freely inside the body. Note that the emission of photo-electrons on the dayside surface of the body is not included in the present study. The radius of the body is set to be RL ~ 2.72 km, which corresponds to a typical size of an asteroid (approximately 1/640 of the radius of the Moon). The ratio of the radius of the body to the gyro radius of thermal ions is RL/ri = 2 (RL/re = 10 to the gyro radius of thermal electrons, RL /di = 1 to the ion inertial length, RL/de = 5 to the electron inertial length, and RL/λDi,e = 100 to the Debye length, respectively). Note that the realistic ratio is RL/ri = 20 for the Moon.

The intrinsic magnetic field (magnetosphere) of the body is modeled by a 2D dipole,

where M denotes the dipole moment and d = r − d0 denotes the distance from the center of the dipole (d0), with \( \mathbf{M}=\left(0,-6.4\uppi {B}_0{R}_L^2\right) \) and d0 = (−0.6RL, −0.3RL). The background magnetic field is imposed in the x − z plane with the direction of 45° with respect to the x-axis (i.e., \( {B}_{x0}={B}_{z0}={B}_0/\sqrt{2} \)), which corresponds to the equatorial spiral angle of the IMF at 1 AU due to the Parker spiral (Parker 1958). This orientation of IMF is chosen in order to include both ion gyro motion and magnetic convection in the simulation (x − y) plane. For the convection of the background magnetic field, an ambient electric field is also imposed in the y direction with the magnitude of Ey0 = USBz0. The initial configuration of the magnetic field lines is seen in Figure 1.

Figure 1

Temporal development of the ion density together with the configuration of magnetic field lines. Spatial profiles of ion density normalized to the initial density N0 at different times are shown by the hot-scale colormap. The magnetic field lines and the location of the inner boundary are indicated by the cyan lines and the magenta circle. The results from the low-resolution run (Umeda and Ito 2014) and the high-resolution run are shown in the top and bottom panels, respectively.

In this paper, we compare two simulation runs. One is the low-resolution run (Umeda and Ito 2014), and the other is the high-resolution run performed on the K computer. In the low-resolution run, the spatial grid size is set to be Δx = Δy = 10λDi and the spatial simulation box is taken to be −8RL ≤ x ≤ 8RL and −12RL ≤ y ≤ 12RL. Thus, a total of Nx × Ny = 160 × 240 grid points are used for the two spatial dimensions. In the high-resolution run, the spatial grid size is set to be Δx = Δy = λDi and the simulation box is taken to be −9.6 ≤ x ≤ 9.6RL and −12.8RL ≤ y ≤ 12.8RL. Thus, a total of Nx × Ny = 1,920 × 2,560 grid points are used for the two spatial dimensions. The velocity space is taken to be −15Vti ≤ vx,y,z ≤ 15Vti for ions and −10Vte ≤ vx,y,z ≤ 10Vte for electrons, with a total of Nvx × Nvy × Nvz = 60 × 60 × 60 grid points for the three velocity dimensions in both runs.

Simulation on the K computer

The high-performance computing infrastructure (HPCI) system research project was started in Japan in 2012 to promote the efficient use of the national flagship K computer. A limited number of research proposals have been accepted as HPCI system research projects using the K computer. The K computer has a total of 82,944 computational nodes, and each node has a single SPARC64 VIIIfx processor (which corresponds to the total of 663,552 CPU cores) and 16-GB memory. The theoretical peak performance of the K computer is about 10 Peta floating operations per second (FLOPS), and this is why the computer system is named as the ‘K’ computer (here, the Japanese unit ‘K’ correspond to 10,000,000,000,000,000 = 10 Peta), The present study used 6,144 computational nodes (i.e., 49,152 CPU cores) which corresponds to about 7% of the entire computational resource.

The low-resolution run requires a total of 400-GB memory for the entire computation, and the run has been performed on cluster supercomputer systems at university supercomputer centers. The high-resolution run requires a total of 50-TB memory for the entire computation. There are a limited number of supercomputer systems which have a total of more than 50-TB memory in Japan (see the TOP 500 supercomputer sites; http://www.top500.org/), although it is impossible to occupy almost the entire resource of these supercomputer systems for more than 3 months. Therefore, the K computer is a unique system for performing the high-resolution run. The number of nodes 6,144 is chosen due to both computational resource and performance. The simulation box of the high-resolution run is decomposed into 96 × 64 nodes with 20 × 40 × 60 × 60 × 60 grid points (which correspond to approximately 8 GB) per node. Our parallel Vlasov code has been well-tuned on the K computer and has achieved the computational performance of 138.634 Tera FLOPS with 6,144 nodes, which corresponds to a performance efficiency of 17.4% to the theoretical peak performance of 786.432 Tera FLOPS = 128 Giga FLOPS × 6,144 nodes (Umeda and Fukazawa 2014).

Figure 1 shows the temporal development of the ion density together with the configuration of magnetic field lines as a result of the interaction between the solar wind and a small dielectric body with a weak intrinsic magnetic field. The magnetic field lines are drawn as contour lines of the vector potential Az. The location of the inner boundary is indicated by the magenta circle at the center of the simulation box.

For the initial condition (ωcit = 0), the simulation box is filled with solar wind plasma propagating from the left to the right. The magnetic field lines are open on the nightside (right-hand side of the body) while they are closed on the dayside. As the time elapses, a very low-density region is formed on the nightside by absorption of the solar wind plasma on the surface of the body, which is known as the wake tail. The configuration of magnetic field lines is also modified to follow the shape of the wake tail. On the dayside, on the other hand, a high-density region is formed which indicates the formation of a bow shock-like structure by the compression of the solar wind plasma by the intrinsic magnetic field of the body. The profile of ion density shows asymmetry with respect to the x-axis for two reasons. One is that the location of the dipole intrinsic magnetic field is slightly shifted to − y direction, and the other is the gyro motion of ions which are reflected at the dayside magnetopause (Umeda and Ito 2014). The system reaches almost the quasi-steady state at ωcit = 6, although this timescale is not enough to follow the trajectories of the reflected ions.

One can notice from the square blocks around the body that the result in the top panels is obtained by the low-resolution run. It may be surprising, even for those who are familiar with plasma physics, that the difference between the two runs is very small although there is minor difference in the plasma environments on the dayside and the magnetic fields on the nightside. However, the present result seems to be consistent with the conventional understanding of plasma physics that the structure and dynamics of global magnetic fields, which are generally described by the MHD equations, are not affected by electron-scale microphysics.

In order to see the difference between the low-resolution and high-resolution runs, we analyze the electric field data. Figure 2 shows the electric field Ex, Ey, and Ez components together with the configuration of magnetic field lines at ωcit = 6. The magnitude is normalized by the ambient motional electric field Ey0.

Figure 2

Spatial profile of the electric fields together with magnetic field lines at ωcit = 6. The magnitudes of the electric field Ex, Ey, and Ez components are shown by the jet-scale colormap. The magnitude is normalized by the ambient motional electric field Eyo. The results from the low-resolution run (Umeda and Ito 2014) and the high-resolution run are shown in the top and bottom panels, respectively.

Although the pattern of the electric fields in both runs looks similar to each other, the difference of the electric fields is much clearer than that of the ion density and the magnetic fields. First, there are fine structures of strong electric field on the surface of the body in the high-resolution run, while they are not so clear in the low-resolution run. The magnitude of the electric field on the surface exceeds 100Ey0 in the high-resolution run, while it is approximately 10Ey0 in the low-resolution run. The electric field on the surface is especially strong at the boundary between open and closed magnetic field lines, because of the difference in the access of solar wind plasma to the surface of the body. It is also seen that the strong electric field develops along magnetic field lines. Second, there are numerical oscillations especially in the electric field Ey component. The first and second differences are due to the spatial resolution (∆x = ∆y). The Debye shielding is not effective, and the Debye sheath does not form in the low-resolution run. Although the magnitude of the electric field Ez component is small (approximately 0.1Ey0), the accumulation of the difference in the electric field Ez component during the propagation of the solar wind results in the small difference in the magnetic field structure on the nightside. It is also seen that there is wave activity along the magnetic field line in the electric field Ez component in the dayside magnetosphere, whose wavelength is estimated as λ ~ 0.4RL. This spatial scale is of the order of the local electron inertial length (i.e., λ ~ 2de,local with Ne,local ~ N0).

It is worth analyzing the source of these electric fields. Figure 3 shows the charge density ρ = e(Ni − Ne), the current density Jx, Jy, and Jz components together with the configuration of magnetic field lines at ωcit = 6. The magnitude of the charge density is normalized by the initial charge density of solar wind ions eN0, and the magnitude of the current density is normalized by eN0Vs.

Figure 3

Spatial profile of the charge and current densities together with magnetic field lines at ωcit = 6. The magnitude of the charge density ρ = e(Ni− Ne) and the magnitudes of current density Jx, Jy, and Jz components are shown by the jet-scale colormap. The magnitude of the charge density is normalized by the initial charge density of solar wind ions eN0, and the magnitude of the current density is normalized by eN0Vs. The results from the low-resolution run (Umeda and Ito 2014) and the high-resolution run are shown in the top and bottom panels, respectively.

The charge density in the vicinity of the surface of the body in the high-resolution run is positive, indicating the formation of the Debye sheath. The charge density in the vicinity of the surface of the body in the low-resolution run is also almost positive, suggesting that the accumulation of the electric charge on the surface of the body occurs but spurious (numerical) oscillations across the IMF are generated because of the absence of the Debye shielding.

There is a fine structure along the magnetic field in the charge density in the dayside magnetosphere (on the closed magnetic field lines) with the wavelength of λ ~ 0.2RL. This spatial scale is of the order of the local gyro radius of ions (i.e., λ ~ 0.5ri,local with |Blocal |approximately 1.3B0). There is also a strong negative charge density on the open magnetic field lines in the vicinity of the magnetosphere (i.e., cusp), since electrons move freely along magnetic field lines and can easily impact the body. The spatial difference of the charge density results in the strong electric field Ex and Ey components on the dayside. These fine structures of the charge density on the dayside are not resolved in the low-resolution run.

Spatial profiles of the current density Jx, Jy, and Jz components in the low- and high-resolution runs look similar to each other, although the magnitudes are larger in the high-resolution run than in the low-resolution run. There is a fine structure in the current density Jx and Jy components at (x, y) ~ (−2RL, −RL) with the wavelength of λ ~ 2RL, which is on the spatial scale of MHD. Note that the fine structure of the electric current is not so obvious in the electric field Ex and Ey components but is indicated in the magnetic field Bz component (not shown here), suggesting that the electric current corresponds to a diamagnetic current due to the compression and decompression of the magnetic field Bz component at the bow shock-like structure.

A preliminary result of the high-resolution (Δx = Δy = λDi,e) global Vlasov simulation on the interaction between the solar wind and a small and weakly magnetized body is presented. The direct comparison of electromagnetic fields and current and charge densities between the high-resolution run and the low-resolution (Δx = Δy = 10λDi,e) run (Umeda and Ito 2014) shows the influence of Debye-scale microphysics on the global structure of a small magnetosphere. The Vlasov simulation also can give us high-resolution distribution functions (Umeda and Ito 2014). Because of the difficulty in the transfer of distribution function data of approximately 30 TB from the K computer to the STE Lab, however, analysis of the distribution function is left as a future study.

The present study suggests that electron-scale kinetics hardly affects the global structure of magnetic field lines. On the other hand, electric fields are strongly enhanced by the electric charge separation and the charge accumulation on the surface of the body. The results are consistent with the conventional understanding of plasma physics that the structure and dynamics of global magnetic fields are generally governed by the MHD equations and are not affected by electrostatic fields on the electron scale, at least in the present two spatial dimensions. Note that the effect of the emission of photo-electrons on the dayside surface of the body should be included in a future study. For a full understanding of the structure and dynamics of magnetosphere, however, a six-dimensional Vlasov simulation (with three spatial dimensions) is essential, because the topology of magnetic fields and convection pattern is modified in three spatial dimensions from those in two spatial dimensions. As the in situ observation shows, there is no bow shock above the Lunar magnetic anomaly (Futaana et al. 2013), while the formation of the bow shock is confirmed in the present global simulation with two spatial dimensions. The full six-dimensional global Vlasov simulation is left as a far future study, since it would require massive computing resources beyond Exascale.

Acknowledgements

The authors are grateful to Prof. Yasuo Ogawa, the EPS editor-in-chief, for inviting us to the EPS frontier letter. One of authors (TU) thanks Yosuke Ito for his assistance in the data analysis. This study is supported by the MEXT/JSPS Grant-In-Aid (KAKENHI) for Challenging Exploratory Research No. 25610144. The high-resolution simulation run has been performed on the K computer at the RIKEN Advanced Institute for Computational Science as HPCI System Research Projects (ID: hp120092, hp140064, hp140081 and hp150069). The low-resolution simulation run has been performed on DELL PowerEdge R815 at Solar-Terrestrial Environment Laboratory (STEL), Fujitsu FX1, and HX600 at Information Technology Center, Nagoya University and Fujitsu RX200S6 at Research Institute for Information Technology, Kyushu University as a computational joint research program at STEL, a collaborative research project of computer science with HPC at Nagoya University, and JHPCN programs at Joint Usage/Research Center for Interdisciplinary Large-Scale Information Infrastructures in Japan.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

TU developed the simulation codes, analyzed the simulation data, and drafted the manuscript. KF helped to perform the computer simulations on the K computer. Both authors contributed to the interpretations and writing of the paper. Both authors read and approved the final manuscript.

Nakagawa T, Kimura S (2011) Role of the solar wind magnetic field in the interaction of a non-magnetized body with the solar wind: an electromagnetic 2-D particle-in-cell simulation. Earth Planets Space 63:477–486View ArticleGoogle Scholar

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.