# Electromagnetic Particle-in-Cell Simulations of the Solar Wind Interaction with Lunar Magnetic Anomalies

Jan Deca, A. Divin, G. Lapenta, Bertrand Lembège, S. Markidis, M. Horányi

► To cite this version:

Jan Deca, A. Divin, G. Lapenta, Bertrand Lembège, S. Markidis, et al. Electromagnetic Particle-in-Cell Simulations of the Solar Wind Interaction with Lunar Magnetic Anomalies. Physical Review Letters, 2014, 112 (15), pp.151102. 10.1103/PHYSREVLETT.112.151102. hal-00984865

HAL Id: hal-00984865

https://hal.science/hal-00984865

Submitted on 12 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Electromagnetic Particle-in-Cell Simulations of the Solar Wind Interaction with Lunar Magnetic Anomalies

J. Deca,¹,²,* A. Divin,³ G. Lapenta,¹ B. Lembège,² S. Markidis,⁴ and M. Horányi⁵

¹Centre for mathematical Plasma Astrophysics (CmPA), Department of Mathematics, KU Leuven, Celestijnenlaan 200B, bus 2400 B-3001 Leuven, Belgium

²Laboratoire Atmosphères, Milieux, Observations Spatiales (LATMOS), Université de Versailles à Saint Quentin, 11 Boulevard D'Alembert, 78280 Guyancourt, France

³Swedish Institute of Space Physics (IRF), SE-751 21 Uppsala, Sweden

⁴High Performance Computing and Visualization (HPCViz), KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden

⁵Laboratory for Atmospheric and Space Physics (LASP), University of Colorado, 1234 Innovation Drive, Boulder, Colorado 80303-7814, USA

(Received 27 November 2013; published 17 April 2014)

We present the first three-dimensional fully kinetic and electromagnetic simulations of the solar wind interaction with lunar crustal magnetic anomalies (LMAs). Using the implicit particle-in-cell code IPIC3D, we confirm that LMAs may indeed be strong enough to stand off the solar wind from directly impacting the lunar surface forming a mini-magnetosphere, as suggested by spacecraft observations and theory. In contrast to earlier magnetohydrodynamics and hybrid simulations, the fully kinetic nature of IPIC3D allows us to investigate the space charge effects and in particular the electron dynamics dominating the near-surface lunar plasma environment. We describe for the first time the interaction of a dipole model centered just below the lunar surface under plasma conditions such that only the electron population is magnetized. The fully kinetic treatment identifies electromagnetic modes that alter the magnetic field at scales determined by the electron physics. Driven by strong pressure anisotropies, the mini-magnetosphere is unstable over time, leading to only temporal shielding of the surface underneath. Future human exploration as well as lunar science in general therefore hinges on a better understanding of LMAs.

DOI: 10.1103/PhysRevLett.112.151102

PACS numbers: 96.50.Ek, 52.25.Dg, 52.65.Rr, 96.20.Jz

Unlike Earth or Mercury, our Moon has no intrinsic magnetic field and is therefore not shielded from the impinging solar wind by a magnetosphere. On the other hand, it does possess small-scale crustal magnetic fields regions, or lunar magnetic anomalies (LMAs), which can range up to a few hundred kilometers in size with surface magnetic field strengths up to hundreds of nanoTeslas [1–6]. Recent work has shown that LMAs can deflect solar wind protons enough to form a so-called mini-magnetosphere, a density cavity shielding the lunar surface from the impinging solar wind plasma [7–15]. Understanding the detailed physics of the solar wind interaction with LMAs, including magnetic shielding and surface charging, is thus vital to evaluate its implications for lunar exploration.

Because of the small spatial scales of LMAs relative to the solar wind ion gyroradius, the dominating physical processes in the interaction are highly nonadiabatic [16–20]. Magnetohydrodynamics (MHD) [8–10] or hybrid [16,17] simulations, lacking the ability to investigate the effects of charge separation, are therefore insufficient for detailed modeling of the near-surface lunar plasma environment. A kinetic model is an absolute must.

Using the fully kinetic and electromagnetic implicit particle-in-cell (PIC) code IPIC3D [21], we provide the answer. By self-consistently simulating the solar wind plasma interaction with LMAs, we identify the electron-dominated

physics near a small-scale magnetic dipole, resembling an idealized LMA plasma environment in which electrons are magnetized in contrast to the ion population. The implicit moment method [22–24] implemented in IPIC3D is designed especially to overcome numerical constraints conventional explicit PIC codes suffer from, hence providing the ideal framework for multiscale LMA simulations under various plasma conditions.

To accomplish this, the original implicit algorithm is modified to accommodate an external dipole magnetic field component **B'**, which is superimposed on the self-consistent (internal) magnetic field **B**:

$$ \mathbf{B}'(\mathbf{r}) = \frac{\mu_0}{4\pi} \left( \frac{3(\mathbf{m} \cdot \mathbf{r})\mathbf{r}}{r^5} - \frac{\mathbf{m}}{r^3} \right), $$

with the source located at the origin and **m** the dipole moment (in A m²). The plane YZ is parallel to the lunar surface, and the X direction is parallel to the unperturbed solar wind flow. The external magnetic field **B**' is introduced in both the particle mover and the field solver (a more specific discussion on the implicit dielectric susceptibility is provided by Ricci, Lapenta, and Brackbill [25]).

Secondly, open boundary conditions (BCs) are modeled by injecting new particles having a thermal spread $v_{th,s}$,
with s the species, and a bulk (solar wind) flow velocity $v_{sw}$, thus creating a uniform drifting Maxwellian plasma. Particles leaving the computational boundaries are excluded from the simulation. One of the boundaries ($x=0$) acts as a perfect absorber and represents the lunar surface. The BCs for the fields are set as follows: $\mathbf{B}_{\text{boundary}} = \mathbf{B}_{\text{IMF}}$ and $\mathbf{E}_{\text{boundary}} = -v_{sw} \times \mathbf{B}_{\text{IMF}}$, where $\mathbf{B}_{\text{IMF}}$ is the solar wind magnetic field.

In the current work, we adopt an electron-ion isothermal solar wind plasma at 1 AU, as a reference, in a temperature and velocity regime close to typical quiet solar wind conditions. The following physical parameters are utilized [26]: the plasma density is $n_0 = 3 \text{ cm}^{-3}$, corresponding to an ion inertial length $d_i \sim 130$ km; the ion and electron temperatures are $T_i = T_e = 15 \text{ eV}$, $v_{sw} = (-600, 0, 0) \text{ km s}^{-1}$, and $\mathbf{B}_{\text{IMF}} = (0, 3, 0) \text{ nT}$. Since the solar wind parameters can fluctuate significantly, we choose to adopt this rather high solar wind temperature to improve numerical stability. The dipole moment $\mathbf{m} = (0, M_d, 0)$, where $M_d = 11.2 \times 10^{12} \text{ A m}^2$, resembles the strongest component of the two-dipole model for the Reiner Gamma magnetic anomaly region by Kurata et al. [11]. The source is placed 13 km below the absorbing lunar surface.

Three-dimensional large-scale simulations with a realistic speed of light and electron mass are rather expensive, and therefore reduced values are used: the ion-to-electron mass ratio is $m_i/m_e = 256$ [27], in which $m_i$ is the reference value of mass for further calculations; and the speed of light value is $c/v_{sw} \approx 59$, creating a mesosonic solar wind flow with $v_{\text{th},i} = 3.8 \times 10^4 \text{ m/s} \ll v_{sw}$ and $v_{\text{th},e} = 6.1 \times 10^4 \text{ m/s} \sim v_{sw}$ (by using a realistic mass ratio, the latter would be $v_{\text{th},e} = 1.6 \times 10^6 \text{ m/s}$). The size of the computational domain measures ($L_x, L_y, L_z$) = (1, 1, 1)$d_i$ [28], with the absorbing surface located at $x=0$ and the dipole placed at the point $(-0.1, L_y/2, L_z/2)$. The grid size is $N_x \times N_y \times N_z = 256^3$ with 64 particles per cell per specie initially. Note that the electron scales are well resolved in our simulation: the electron skin depth $d_e = 0.0625 d_i = 16\Delta x$, with $\Delta x$ the grid spacing. The time step is set relative to the ion plasma frequency: $\Delta t = 0.0375\omega_{pi}^{-1}$; that is, thermal electrons pass only $0.2\Delta x$ in $16\Delta t$. Our numerical scheme does not require us to resolve the Debye length [22–24]. Throughout this Letter, the length and time scales are referenced to the initial plasma conditions and not to the local values. Note that the inertial length scales $d_i$ and $d_e$ are considered as a convenient reference frame only. The plasma dispersion scales $r_G$ and $r_{DM}$, however, are at all times resolved by our numerical model [29,30].

Figure 1 shows the 2D electron (top) and ion (bottom) charge density profiles along the dipole axis (XY plane), after the simulation has reached a quasisteady state. Superimposed are magnetic field lines. As the solar wind impinges on the dipolar structure, perpendicular to the lunar surface, both ions and electrons are deflected and

<figure><img src="image_2.png" /><figcaption>FIG. 1 (color online). 2D electron (top) and ion (bottom) charge density profiles, scaled to the initial density $n_0$ and along the dipole axis (Y direction) at $z = L_z/2$, after the simulation has reached a quasisteady state. The solar wind is flowing perpendicular (in the $-X$ direction) to the lunar surface at $x=0$. Superimposed in white are magnetic field lines. Note that only the lower half of the simulation domain is shown in the X direction.</figcaption></figure>

a density cavity is created. Surrounding this lower density region, a higher density halo is instantly built up, the size of which is about one electron skin depth (~8 km) thick with up to 5 times the initial solar wind density towards the cusp regions: a mini-magnetosphere has formed. The halo region consists of solar wind particles which are temporarily trapped by the dipole field when the magnetic pressure equals the solar wind pressure. We do not observe a clear shock associated with the mini-magnetosphere structure.

Given the scale size of the anomaly, only the electron population is magnetized, whereas the ions, due to their higher mass, should only feel the dipole field very close to the lunar surface before being scattered nonadiabatically. Ions thus easily penetrate the density halo and create a charge separation between both species. The latter results in the generation of a large normal electric field $\mathbf{E}_n$ directed in the X direction at the subsolar point. Almost instantly after the start of the simulation, $\mathbf{E}_n$ becomes large enough to also deflect the ions. The formation of the mini-magnetosphere is therefore mainly an electrostatic effect.

Making the analogy with Earth’s magnetosphere [8,17], one could refer to the higher density barrier as an electro-sheath rather than the magnetosheath, because the mini-magnetosphere structure is formed under the impulse of the electron dynamics. As indicated by Kallio et al. [17], the analogy is not complete due to the absence of a bowshock.

The mini-magnetosphere shape is asymmetric in the direction perpendicular to the dipole axis (Fig. 2). Indeed,
<figure><img src="image_3.png" /><figcaption>FIG. 2 (color online). Ion charge density profiles, scaled to the initial density $n_0$, along three different planes: the XY plane (upper left) and XZ plane (lower left) through the dipole center and the YZ plane (right) at $x = 0.01d_i$ above the lunar surface. Note that only the lower half of the simulation domain is shown in the X direction.</figcaption></figure>

under influence of the $\nabla B$ and $E \times B$ drifts the trapped electrons are flowing in the +Z direction along the entire length of the electrosheath (see also Fig. 3). The much heavier ions, on the other hand, are deflected on all sides of the dipole structure by the normal electric field $E_n$. Subsequently, an asymmetric density cavity or halo is

created. Note that neither MHD nor hybrid simulations can correctly model this configuration, because the process is initiated by the electrons having highly non-Maxwellian velocity distributions near the mini-magnetosphere structure.

Figure 3 presents the electron, respectively, ion velocity distributions, and pressure profiles along the direction

<figure><img src="image_4.png" /><figcaption>FIG. 3 (color online). Profiles along the direction parallel to the solar wind flow and through the center of the dipole. The upper panels hold the magnetic and kinetic pressure profiles for the electron (left) and ion (right) populations, in code units. The middle panel shows the Z component of the electron velocity distribution (left) and the X component of ion velocity distribution (right), normalized to $v_{sw}$. The lower panels present the density profiles (left), normalized to the initial density $n_0$, and the mirror instability criterium (right) with the area of interest colored in magenta.</figcaption></figure>
parallel to the solar wind flow and through the center of the dipole. The electron population is heated and deflected inside the electrosheath, while most of the ions still penetrate the latter and reach the surface almost unaffected (middle panels). About 10% of the incident solar wind ions are reflected by the density halo, consistent with observations by, e.g., Saito et al. [18] and simulations by, e.g., Kallio et al. [17]. Note that the density halo starts forming at $0.001d_i(\lesssim 13 \text{ km})$ above the lunar surface when the ion dynamic pressure equals the magnetic pressure. The density peaks slightly closer to the surface. In between those points, the ion kinetic pressure is converted into thermal pressure.

The solar wind magnetic field direction $B_{IMF}$ is opposite to the dipole magnetic field along the line $y = L_y/2$, $z = L_z/2$, hence creating a zero point in the total magnetic field configuration at $0.27d_i$ above the surface. Notably, we do not observe any particle flows associated with magnetic reconnection, indicating that the mini-magnetosphere electrosheath currents cannot shield off the dipole field completely, because of either the absorbing surface or the strongly nonadiabatic behavior of the ions. Increasing the magnetic dipole moment (not reported in the present Letter), however, can move the neutral line close to the density halo, producing favorable conditions for electron acceleration by the solar wind electric field and bringing more resemblance with the conventional large-scale magnetosphere.

Finally, we note that the mini-magnetosphere is unstable over time. Occasionally, blobs of plasma are able to enter the density cavity and reach the surface. It is found (Fig. 3, lower right panel) that the configuration is slightly unstable to the mirror instability [31]:

$$ \beta_{\parallel} - \frac{\beta_{\perp}^2}{1 + \beta_{\perp}} < 0. $$

In conclusion, using the implicit particle-in-cell code iPic3D we have produced the first fully kinetic 3D simulation of the solar wind interaction with a crustal lunar magnetic anomaly, forming a mini-magnetosphere. We have shown that the configuration is driven by electron motion, because the LMA scale size is small with respect to the gyroradius of the solar wind ions. We found that a mirror instability develops at the outer edge of the density halo, allowing both electron and ion plasma to reach the lunar surface underneath the dipole field. Hence, the fully kinetic treatment identifies electromagnetic modes that alter the magnetic field at scales determined by the electron physics, and this is a fundamental discovery. Our work now opens new frontiers of research, allowing us to study this process in more detail to evaluate the implications for surface effects, waves, and instabilities, particle acceleration in the region, and its impact on exploration in terms of disturbances of instruments and lunar exploration in

general. Finally, the simulation results are ideally suited to be compared with field or particle observations from spacecraft such as Kaguya (SELENE), Lunar Prospector, or ARTEMIS.

This research has received funding from the European Commission's FP7 Program with the grant agreement SWIFF (Project No. 2633430 [32]) and EHEROES (Project No. 284461 [33]). The simulations were conducted on the computational resources provided by the PRACE Tier-0 Project No. 2011050747 (Curie supercomputer). This research was supported by the Swedish National Space Board, Grant No. 136/11. J. D. has received support through the HPC-Europa2 visitor program (Project No. HPC08SSG85) and the KuLeuven Junior Mobility Program Special Research Fund. Test simulations were performed by A. D. on resources provided by SNIC at PDC Centre for High Performance Computing (PDC-HPC), Grants No. m.2013-1-108 and No. m.2013-1-34.

<footnote>* jan.deca@wis.kuleuven.be</footnote>

<footnote>[1] M. Fuller, Rev. Geophys. 12, 23 (1974).</footnote>

<footnote>[2] P. Dyal, C. W. Parkin, and W. D. Daily, Rev. Geophys. Space Phys. 12, 568 (1974).</footnote>

<footnote>[3] C. T. Russell, P. J. Coleman, Jr., B. R. Lichtenstein, G. Schubert, and L. R. Sharp, in Space Research XIV, edited by M. J. Rycroft and R. D. Reasenberg (Akademie-Verlag, Berlin, 1974), pp. 629-634.</footnote>

<footnote>[4] L. L. Hood, C. T. Russell, and P. J. Coleman, J. Geophys. Res. 86, 1055 (1981).</footnote>

<footnote>[5] R. Lin, K. Anderson, and L. Hood, Icarus 74, 529 (1988).</footnote>

<footnote>[6] D. L. Mitchell, J. S. Halekas, R. P. Lin, S. Frey, L. L. Hood, M. H. Acuña, and A. Binder, Icarus 194, 401 (2008).</footnote>

<footnote>[7] R. P. Lin, D. L. Mitchell, D. W. Curtis, K. A. Anderson, C. W. Carlson, J. McFadden, M. H. Acuna, L. L. Hood, and A. Binder, Science 281, 1480 (1998).</footnote>

<footnote>[8] E. M. Harnett and R. Winglee, J. Geophys. Res. 105, 24997 (2000).</footnote>

<footnote>[9] E. M. Harnett and R. M. Winglee, J. Geophys. Res. 107, 1421 (2002).</footnote>

<footnote>[10] E. M. Harnett and R. M. Winglee, J. Geophys. Res. 108, 1088 (2003).</footnote>

<footnote>[11] M. Kurata, H. Tsunakawa, Y. Saito, H. Shibuya, M. Matsushima, and H. Shimizu, Geophys. Res. Lett. 32, L24205 (2005).</footnote>

<footnote>[12] J. S. Halekas, G. T. Delory, D. A. Brain, R. P. Lin, and D. L. Mitchell, Planet. Space Sci. 56, 941 (2008).</footnote>

<footnote>[13] M. Wieser, S. Barabash, Y. Futaana, M. Holmström, A. Bhardwaj, R. Sridharan, M. B. Dhanya, A. Schaufelberger, P. Wurz, and K. Asamura, Geophys. Res. Lett. 37, L05103 (2010).</footnote>

<footnote>[14] C. Lue, Y. Futaana, S. Barabash, M. Wieser, M. Holmström, A. Bhardwaj, M. B. Dhanya, and P. Wurz, Geophys. Res. Lett. 38, L03202 (2011).</footnote>

<footnote>[15] R. A. Bamford, B. Kellett, W. J. Bradford, C. Norberg, A. Thornton, K. J. Gibson, I. A. Crawford, L. Silva, L. Gargaté, and R. Bingham, Phys. Rev. Lett. 109, 081101 (2012).</footnote>
[16] A. R. Poppe, J. S. Halekas, G. T. Delory, and W. M. Farrell, J. Geophys. Res. **117**, A09105 (2012).

[17] E. Kallio, R. Jarvinen, S. Dyadechkin, P. Wurz, S. Barabash, F. Alvarez, V. A. Fernandes, Y. Futaana, A.-M. Harri, J. Heilimo, C. Lue, J. Mäkelä, N. Porjo, W. Schmidt, and T. Siili, Planet. Space Sci. **74**, 146 (2012).

[18] Y. Saito, M. N. Nishino, M. Fujimoto, T. Yamamoto, S. Yokota, H. Tsunakawa, H. Shibuya, M. Matsushima, H. Shimizu, and F. Takahashi, Earth Planets Space **64**, 83 (2012).

[19] X. Wang, M. Horányi, and S. Robertson, J. Geophys. Res. **117**, A06226 (2012).

[20] X. Wang, C. T. Howes, M. Horányi, and S. Robertson, Geophys. Res. Lett. **40**, 1686 (2013).

[21] S. Markidis, G. Lapenta, and Rizwan-uddin, Math. Comput. Simul. **80**, 1509 (2010).

[22] R. J. Mason, J. Comput. Phys., **41**, 233 (1981).

[23] J. Brackbill and D. Forslund, J. Comput. Phys. **46**, 271 (1982).

[24] G. Lapenta, J. U. Brackbill, and P. Ricci, Phys. Plasmas **13**, 055904 (2006).

[25] P. Ricci, G. Lapenta, and J. U. Brackbill, J. Comput. Phys. **183**, 117 (2002).

[26] M. G. Kivelson and C. T. Russell, in *Introduction to Space Physics*, edited by M. G. Kivelson and C. T. Russell (Cambridge University Press, Cambridge, England, 1995), p. 586.

[27] A. Bret and M. E. Dieckmann, Phys. Plasmas **17**, 032109 (2010).

[28] Initially low-resolution simulations were performed on a larger domain to ensure that the current size does not influence the plasma evolution in the box.

[29] V. M. Gubchenko, Sov. J. Plasma Phys. **14**, 186 (1988).

[30] V. Gubchenko, H. Biernat, and M. Goossens, Radiophys. Quantum Electron. **44**, 72 (2001).

[31] D. J. Southwood and M. G. Kivelson, J. Geophys. Res. **98**, 9181 (1993).

[32] http://swiff.eu.

[33] http://www.eheroes.eu.