Research ArticleSPACE SCIENCES

Was the moon magnetized by impact plasmas?

See allHide authors and affiliations

Science Advances  02 Oct 2020:
Vol. 6, no. 40, eabb1475
DOI: 10.1126/sciadv.abb1475


The crusts of the Moon, Mercury, and many meteorite parent bodies are magnetized. Although the magnetizing field is commonly attributed to that of an ancient core dynamo, a longstanding hypothesized alternative is amplification of the interplanetary magnetic field and induced crustal field by plasmas generated by meteoroid impacts. Here, we use magnetohydrodynamic and impact simulations and analytic relationships to demonstrate that although impact plasmas can transiently enhance the field inside the Moon, the resulting fields are at least three orders of magnitude too weak to explain lunar crustal magnetic anomalies. This leaves a core dynamo as the only plausible source of most magnetization on the Moon.


The Moon presently lacks a core dynamo magnetic field. However, it has been known since the Apollo era that the lunar crust contains remanent magnetization, with localized surface fields reaching up to hundreds of nanoteslas or higher and spanning up to hundreds of kilometers (1). Magnetic studies of Apollo samples and the lunar crust indicate that the magnetizing field likely reached tens of microteslas before 3.56 billion years (Ga) ago (1, 2). The origin of the strongest lunar crustal anomalies and the source of the field that magnetized them have been longstanding mysteries.

Although magnetic fields in rocky bodies are commonly explained by convective dynamos in their metallic cores, a convective dynamo on the Moon may not have had sufficient energy to produce the strongest implied surface paleofields (3, 4). This may imply that a fundamentally different nonconvective dynamo mechanism operated in the Moon or that a process other than a core dynamo produced such magnetization.

Hypothesized nondynamo fields include those of the ancient Earth, the solar nebula, and the solar wind. However, by 3.56 Ga ago, the nebula had already long dissipated (5), while the Earth’s field even at Earth’s Roche limit would only reach ~1 μT, assuming a surface field like that of the Earth today. The solar wind field, also known as the interplanetary magnetic field (IMF), is thought to have only reached up to 10 to 30 nT at the Moon 3.56 Ga ago (fig. S1). However, local compression of the solar wind could amplify the IMF.

In particular, hypervelocity asteroid impacts can vaporize and ionize lunar crustal materials, releasing plasmas directly into the wind (6, 7). Motivated by the observation that some of the strongest and largest lunar crustal anomalies are located at the antipodes of four young large (>600 km in diameter) basins (8), it was hypothesized that as impact plasmas engulf the Moon, they compress the IMF and amplify the induced field inside the Moon, leading to an enhanced crustal field at the antipode. These amplified fields were invoked as the source of the high field intensities needed to explain the antipodal magnetization. Although these crustal anomalies could also be the product of a dynamo field recorded by highly magnetic materials like metal-rich impactor detritus (2) and/or dikes (9), the antipodal fields hypothesis has remained a leading alternative to the dynamo hypothesis for lunar magnetization for the last several decades (1013). Furthermore, variants of this hypothesis, in which fields are generated by transient impact plasmas, have been invoked as the source of magnetization found within smaller craters on the Moon (14, 15), as well as on asteroids and meteorite parent bodies (1619).

However, it has not been conclusively shown that the antipodal impact field mechanism can supply the required paleointensities. Post-impact antipodal magnetization could be acquired as thermoremanent magnetization (TRM) when hot basin ejecta lands and cools at the antipode ~1 to 24 hours after the impact and/or as shock remanent magnetization (SRM) due to shocking of antipodal materials by the landing ejecta (20). For TRM, this would indicate a minimum paleointensity of ~20 μT given the TRM susceptibility of mafic impact-melt breccias and assuming the thickest possible magnetized crustal layers (20 km) (2). SRM susceptibilities are typically about three times lower than TRM susceptibilities (10, 21), which would indicate a minimum paleointensity of 60 μT. If other materials (like mare basalts) and/or thinner magnetized layers (e.g., 2 km) are responsible for the observed surface fields of the lunar anomalies, the paleointensities would be at least an order of magnitude larger than these values (2). As such, we conservatively adopt a range of 20 to 60 μT as the minimum required paleointensity, which is also within the range of laboratory measurements of paleointensities measured from Apollo samples older than 3.5 Ga (1, 22).

Before a basin-forming impact, two coupled magnetic fields would have been present: the IMF and the IMF-induced field in the conducting lunar interior. Following the impact, the expanding, conductive vapor would carve out a region of low to null field within the solar wind (7, 20). This change in the external field would force a change in the coupled internal field. The more resistive outer layers of the Moon (crust and upper mantle) would respond to this change on a time scale of seconds to minutes, while the fields in the highly conductive core and lower mantle would remain almost unchanged over time scales of years and days, respectively (Supplementary Materials). It was suggested (7, 20) that the field lines in the outer layers of the Moon would wrap around the deeper layers and emerge as a concentrated bundle at the antipode, thereby creating a locally amplified field. By invoking conservation of magnetic energy, the field amplification was estimated from the ratio of the lunar cross section to the estimated area of converged field lines (20), giving a proposed antipodal amplification factor of at least 300. Assuming that the steady component of the initial induced field was 30 nT, this would give an antipodal field of at least ~10 μT. This field was proposed to be sustained within vapor at the antipode over a time scale of a day, encompassing the arrival times of most ejecta (20).

This predicted field strength is consistent with the required paleointensities and the magnetization acquisition times. However, there remain key theoretical gaps that have not been previously addressed. Previous studies focused on hydrodynamic and shock physics simulations (6, 7, 20) and did not explicitly model the magnetic field. Instead, the field amplification was obtained from simplified geometric arguments and conservation of energy as described above. This approach is unsuitable for capturing four important mechanisms. First, the magnetic field and plasmas are coupled by a nonlinear process that cannot be analytically described for complex field topologies. Second, any changes in the induced field direction inside the Moon are due to field induction, which may involve ohmic dissipation, such that magnetic energy is not conserved during the cloud expansion. Third, modeling the motion of the vapor as that of a neutral gas (7, 20) neglects plasma processes that could accelerate it beyond the speeds predicted by hydrodynamics. Fourth, the initial induced field inside the conducting body, previously assumed to equal that of the mean IMF just before the impact, would instead be the vector average over the diffusion time scale of the variable IMF, which is much less than the instantaneous IMF.

This work addresses these gaps by introducing self-consistent modeling of post-impact plasmas and magnetic fields, including field diffusion and dissipation inside the Moon along with revised analytical considerations of the induced field. We combine shock physics simulations of basin excavation and vapor generation with magnetohydrodynamic (MHD) simulations of the coupled interaction of plasmas and fields. Impact basin-forming simulations were performed using the shock physics code iSALE-2D (23), a multimaterial, multirheology code based on the simplified arbitrary Lagrangian-Eulerian (SALE) hydrocode in two dimensions (see Materials and Methods). The resulting vapor was used to drive three-dimensional (3D) MHD simulations that include the interaction of the Moon, the solar wind, and the vapor. For the MHD simulations, we used the Block Adaptive Tree Solar-Wind Roe Upwind Scheme (BATS-R-US) code (24), which is capable of modeling the magnetic field evolution inside resistive bodies (see Materials and Methods) (25). Recent improvements to BATS-R-US enable more efficient coupling between the internal and external fields, permitting us to simulate a diversity of impact scenarios (table S1).

We focus on impact events that form basins the size of Imbrium whose antipodal region contains some of the strongest magnetic anomalies observed from orbit (26). We consider a 60-km-radius impactor vertically striking a 30-km-thick crust at a typical speed of 17 km s−1 (27). The basin formation process is simulated using parameters and methods similar to those that have successfully reproduced recent spacecraft observations of the Moon such as crustal thickness maps (28) and crust and mantle constitutive models (Materials and Methods) (29). The vapor generation and basin excavation are shown in fig. S4. The total vaporized mass produced by the simulations is consistent with that inferred analytically in (7, 30) (Materials and Methods).

In the MHD simulations, we first allow the body and the solar wind flow to reach a steady state and then launch the vapor into the domain, in accordance with the iSALE-2D output, through a spherical cap on the surface approximating the basin (Materials and Methods). The body’s internal conductivity is prescribed based on the measured radial lunar conductivity profile (Supplementary Materials). The incoming solar wind’s properties correspond to those estimated for the Sun at 0.7 to 0.8 Ga after solar system formation, the approximate time of Imbrium’s formation (Supplementary Materials), and remain constant during the simulation (table S1). The initial field inside the body is assumed to be a uniform induced field equal to the surrounding IMF, which is conservatively chosen to be 30 nT, similar to (20).

This field value is an overestimation for two reasons. First, the IMF itself was more likely only ~10 nT at that epoch (Supplementary Materials). Second, the IMF varies on time scales much shorter than the time it takes it to diffuse into the core (tens of years; fig. S2). With respect to the latter, the vector average of the IMF over tens of years amounts to just 1% of the mean IMF magnitude (31), giving 0.1 to 0.3 nT in the core. Furthermore, the lunar rotation, orbital motion, and passage through the Earth’s geomagnetic tail further change the external magnetic field and its orientation with respect to the body, which would further decrease the mean induced field in the core. Nonetheless, we adopt an initial value of 30 nT throughout the body both to enable us to directly compare our results with those of (20) and to conservatively obtain an absolute upper limit on the amplified fields.


The results of our reference simulation, referred to as case 1, are shown in Fig. 1 and movie S1. As predicted in (7, 20), the expanding impact plasma creates a magnetic cavity and enhances the IMF at its periphery. This field compression occurs due to the deceleration of the solar wind as it encounters the dense expanding vapor, causing the IMF carried by the wind to pile up against the vapor. This process converts the kinetic energy stored in the bulk flow of the wind into magnetic energy (fig. S3). While the amplified IMF slightly diffuses into the crust, three important counteracting effects emerge from these self-consistent 3D simulations that were not previously considered.

Fig. 1 Plasma flow and magnetic field evolution following a basin-forming impact on the Moon.

Snapshots are extracted at 10, 50, 150, and 300 s after impact in the plane containing the impact vector (−z direction), solar wind flow (+z direction), and the IMF (+x direction). The impact location is at (x, y, z) = (0, 0, 1) Rm. The left panels show the plasma density (color contours) and velocity (white arrows, scaled to the speed and pointing in flow direction). The middle panels show the magnetic field magnitude (color contours) and vector (black arrows, scaled to magnitude and pointing in field direction). The right panels show diagrams highlighting the factors controlling the field evolution at each snapshot. The arrows marked by U and B are the solar wind velocity and IMF direction, respectively. The full evolution can be viewed in the Supplementary Materials (movie S1).

First, the Moon’s resistive outer layers (crust and upper mantle) destroy magnetic flux at a rate comparable to the expansion rate of the vapor. As discussed below, the loss rate of the magnetic field is consistent with theoretical estimations of the dissipative power of the crust, which removes magnetic energy from the system. Second, 3D diffusion of the field in the mantle and crust allows the field to slip around the core rather than it being anchored in the core. As a result of these two effects, at 300 s after the emergence of vapor from the impact basin (Fig. 1, bottom), when the vapor has almost completely engulfed the Moon, the field in the crust and upper mantle is actually weaker than the initial field induced by the IMF (rather than stronger as proposed by the impact plasmas hypothesis). While (7, 20) proposed that the field lines inside the Moon converge without loss at the antipode, the simulations show that this is prevented by the combination of dissipative losses and 3D diffusive transport. As such, our results do not support the arguments used to calculate amplification factors of >300 in (20) because there is no conservation of magnetic energy or field convergence. Instead, we find a localized and transient amplification factor of <2, which occurs just ~50 s after the impact (Fig. 2). Therefore, the amplified field is too weak to explain the antipodal anomalies and it reaches its peak value long before the estimated >~1-hour arrival time (20) of the bulk of the ejecta at the antipode.

Fig. 2 Magnetic field at the time of maximum field for the simulation in Fig. 1.

(A) 3D view at 50 s after impact. The spherical surface at the center is the lunar surface. The transparent yellow surface is an isosurface of density of 107 cm−3, approximating the shape of the cloud periphery. The color contours show the magnetic field on the lunar surface and in the x-z and y-z planes, and the black contours show the Moon-centric distance in lunar radii, Rm. The point of view was chosen to overlook the area antipodal to the impact (red cross). (B) Magnetic field as a function of time. (Top) Mean field inside the Moon as a function of time. (Bottom) Maximum field found inside the crust (upper 5% of radius of Moon) as a function of time.

Third, the strongest compressed fields at the vapor periphery do not actually reach the crust but are rather pushed away from the body. At the same time, the incoming solar wind flow is blocked and the IMF becomes disconnected from the body (fig. S3). After this occurs, the total magnetic energy in the body can only decrease because it is starved from advective resupply by the IMF (Fig. 2). This provides yet another reason that plasma amplified fields cannot account for the crustal magnetization: the strongest amplification occurs far above the surface of the Moon.

Apart from the three factors limiting the field listed above, a fourth key factor that would further decrease the resulting field in a realistic setting is that the initial (pre-impact) field induced by the IMF in the core and lower mantle would be lower due to the constant changes in the IMF’s direction and would instead equal the mean value of the IMF over their magnetic diffusion time scales (months to years). As such, in reality, we expect the initial mean IMF to be at most 10% and 1% of the 30-nT value assumed here in the mantle and core, respectively. This would further reduce the final mantle and crustal fields by at least an order of magnitude from the values obtained in these simulations.

An additional loss mechanism that might have limited the antipodal field is magnetic reconnection. In (7, 20), the field geometry envisioned to exist at the antipode would have given rise to antiparallel field lines in the plasma just outside the antipode. In actuality, we found in our simulations that magnetic reconnection of these fields does not take place because the antiparallel field geometry never forms in the first place; any magnetic flux pushed toward the antipode either dissipates inside the Moon or is advected away by the vapor.

Parameter space study of different impact scenarios

The simulation in Figs. 1 and 2 is for only one possible choice for the IMF direction, solar wind speed, impact location, and impact cloud physical properties. We modeled seven additional endmember scenarios with different combinations of parameter values (table S1, Fig. 3, and fig. S6). The maximum crustal fields for all cases were calculated as in Fig. 2B (table S1).

Fig. 3 Plasma flow and magnetic field evolution following four different impact scenarios (cases 2, 4, 6, and 7).

Snapshots from 50 s after launch of the vapor into the MHD simulations (table S1) are shown. Layout, color coding, and symbols are as in Fig. 1, except that the right column depicts the initial conditions, where U and B are the solar wind velocity and IMF direction, respectively. (A) Impact on upwind side (case 2). (B) IMF parallel to that of the solar wind flow (case 4). (C) Lunar crust and mantle with enhanced conductivities (case 6). (D) Colder vapor and faster wind (case 7).

Cases 2, 4, and 8 explored alternative impact locations and relative orientations of the IMF and solar wind velocity. In case 2 (Fig. 3A), the impact is on the upwind side of the Moon. The compressed IMF is pushed away from the body, never reaching the crust. Because the crust is now shielded from the IMF, the initial induced field decays. In case 4 (Fig. 3B), the IMF is parallel to the solar wind flow direction (the vapor temperature is also decreased; see details in case 3 below), and the compressed field again does not occur at the surface. In case 8 (fig. S6C), the ambient plasma is static, leading to a substantially weaker compression at the vapor periphery. This approximates an impact occurring while the Moon is inside the Earth’s magnetotail. These cases demonstrate that the crustal field amplification magnitude and its time of occurrence after impact depend on the geometry of the Sun-Moon-Earth-impactor-IMF system. Among cases 1, 2, 4, and 8, the highest amplification occurs in case 1.

The largest overall amplifications in the crust were found in cases 3, 5, and 7. These have the same impact location and relative orientation of the IMF and solar wind velocity as in case 1 but with different choices of physical parameters for the solar wind and vapor. In case 3 (fig. S6A), the initial vapor temperature in the MHD model was reduced from that specified by the iSALE-2D output (2000 K) to 500 K. This does not change the ionization state of the plasma: the plasma is still perfectly conducting and described by the ideal MHD approximation. This leads to slower expansion and gives the compressed IMF more time to diffuse into the Moon. Although idealized, artificially changing the initial vapor temperature mimics the effects of friction with impact melt, which may slow down the vapor (20). In case 5, the solar wind speed is increased from 400 to 1000 km s−1, increasing the energy available for IMF compression at the vapor periphery. Both cases increase the enhanced field induced in the crust (table S1). Case 7 is constructed with both the faster wind speed and the cooler vapor, which yields the largest amplified field in the crust (107 nT). Case 6 assumes outer layers with higher conductivity, exhibiting less amplification than in case 1. In this case, there is less flux destruction by the crust but the field lines are more inhibited from being wrapped around the moon.

The crustal field amplification for all cases lasts for tens to hundreds of seconds, which does not coincide with the time needed to acquire TRM. Therefore, in Fig. 4, we compare our results to the minimum required SRM paleointensity (60 μT). Even in the most favorable scenario (case 7), in which the maximum transient field is 3.6 times higher than the initial field, this field is only 0.0025 of the required paleointensity.

Fig. 4 The maximum predicted crustal amplified field compared to the paleointensities of the fields that magnetized the Moon.

Red arrows mark the maximum enhanced fields for each of the eight simulation cases, each of which differs by one or two parameters from the baseline (case 1). From left to right, these are baseline simulation (case 1), impact location on upwind side of Moon (case 2), colder impact vapor (case 3), IMF parallel to solar wind velocity (case 4), faster solar wind (case 5), higher conductivity of crust and mantle (case 6), faster solar wind and colder impact vapor (case 7), and no solar wind flow (case 8). Detailed descriptions of the cases are given in table S1. The blue solid line marks the minimum required paleointensities. The black solid line marks the initial induced internal field used in the simulations (30 nT; an extreme upper limit). The black dashed line marks the more plausible initial value (1 nT) based on the vector mean of a realistic IMF at 3.9 Ga ago.


Level of field enhancement due to the vapor expansion into the solar wind

The MHD simulations (Figs. 1 and 3 and fig. S6) demonstrate that the vapor expansion causes the IMF carried by the solar wind to be enhanced because the vapor constitutes an obstacle to the wind, causing it to decelerate and pile up. The source of the compressed IMF magnetic energy is the bulk kinetic energy of the upstream wind. Initially, a pile-up region is created in which the magnetic energy density equals the kinetic energy density in the undisturbed wind flow (for the baseline case; fig. S3, top). This gives a field of ~60 nT inside the pile-up region (Fig. 1, top). As the vapor continues to expand (fig. S3, middle), it sweeps up increasingly more of the incoming IMF, leading to an amplification factor of up to 3 to 4. This level of amplification of the IMF is consistent with observed factors of 2 to 6 amplification [depending on the solar wind pressure (32, 33)] in pile-up regions ahead of comets and Venus’ ionosphere and is lower than the IMF compression ratio of 15 estimated for impact plasmas on the Moon in (7).

Examination of the simulations in our parameter space study further demonstrates that it is the solar wind kinetic energy that controls the level of compression. In case 8 (fig. S6), the outside plasma is not flowing and the field enhancement there is small and diffuse. In cases 5 and 8, the wind speed is increased compared to the baseline (table S1) and the enhancement increases to up to a factor 4.6 in the wind pile-up region (140 nT). In addition, in cases 3 and 7, in which the vapor speed is reduced, the compressed field is similar in magnitude to that obtained in the baseline case. Nonetheless, even the amplified fields obtained at the cloud periphery do not reflect the field enhancement induced in the crust because of the crustal resistivity, the effects of the magnetic cavity, and the transient nature of cloud expansion.

Magnetic dissipation counteracts crustal field enhancements

The main factor inhibiting field enhancement inside the Moon is the resistivity of the crust. While the induced field inside the Moon is being reconfigured due to the conducting cloud expansion (Fig. 1), it is also being dissipated in the upper mantle and crust due to ohmic dissipation. This mechanism operates whenever the magnetic field in a resistive material is not curl free (see derivation in the Supplementary Materials). Ohmic dissipation converts magnetic energy into heat with an energy density loss rate of W = J2/σ, where J is the current density and σ is the conductivity. In turn, μJ = ∇ × B, where B is the magnetic field vector and μ is the local magnetic permeability (Supplementary Materials). Once the Moon is exposed to the weaker fields of the expanding vapor, the current density inside the Moon increases due to the change in the field direction between the interior and exterior of the Moon and dissipation starts to dominate over field transport. The resulting magnetic field evolution takes on a spatially complex structure that is reflected in the simulations, leading to removal of flux from the crust and upper mantle (Fig. 1).

Estimating the magnetic dissipative power of the crust

We verified that the amount of simulated magnetic energy lost by the crust in the simulations is realistic by comparing the simulations with the total dissipative power estimated analytically from the crustal thickness and resistivity (Supplementary Materials). We estimate a lower bound for the dissipated power by taking the crust to constitute the outer 5% of the Moon and assuming a resistivity of 106 ohm·m. This should underestimate the resistive losses because this value is 100 times lower than the measured lunar crustal resistivity; this conservative value was adopted in the simulations for the purpose of numerical stability (see Supplementary Materials). To obtain a lower bound on induced currents, we further assume that the spatial change in the field occurs smoothly and over the largest possible length scale, from 30 nT at the core-mantle boundary to 3 nT at the surface. This would give rise to a dissipative power of ~7 × 108 W (eq. S6 in the Supplementary Materials). By comparison, the initial magnetic energy contained inside the Moon is of the order of 8 × 109 J, assuming a uniform field of 30 nT integrated over the entire lunar volume. Thus, the crust of the Moon can very effectively reduce the magnetic energy when exposed to a magnetic cavity. When totally engulfed by the cavity, it could remove magnetic energy that amounts to 10% of the initial energy every second.

The fact that the crust of the Moon can serve as such an effective sink of magnetic energy may seem unexpected because it is known that the IMF is free to diffuse through its outer layers and into the lunar wake without loss (31). This is because magnetic dissipation only occurs when the magnetic field is not curl free; for a uniform field flowing past the body, only lossless magnetic diffusion would take place on the upwind side (Supplementary Materials). In addition, without an impact event, magnetic flux would continue to be constantly replenished by the incoming wind. In contrast, after the impact, vapor expansion causes the incoming IMF to change direction such that the Moon is gradually magnetically isolated from the IMF (Supplementary Text and fig. S3).

Interplay of field enhancement and field removal during vapor expansion

The field removal over the time observed in the simulations is consistent with the above theoretical lower bound. In case 1 (Fig. 1), the vapor-induced cavity is initially located near the impact site and, by 10 s, there is significant reduction in the field inside the crust adjacent to the impact. As the vapor expands toward the antipode, more and more magnetic energy is being destroyed in the crust. At the same time, magnetic energy from the deeper layers diffuses into the now low-field crust, and this energy too is progressively removed. Although the enhanced field of the compressed solar wind at the cloud periphery diffuses into the crust (Figs. 1 to 3 and fig. S6), this is not enough to counteract the dissipative losses in the crust. In particular, the diffusion time of the crust is so short (~10−5 s; fig. S2) that as the leading edge of the crustal enhancement advances with the advancing outer edge of vapor cloud sliding along the surface, the enhanced crustal field in the trailing edge disappears soon after the arrival of the cavity.

Lastly, by 300 s, the vapor engulfs the Moon and cuts off any new solar wind fields from reaching the lunar surface. Only the most conducting layers (e.g., the lower mantle and core) still retain the initial field magnitude (Fig. 1, bottom). This comparison shows that exposing the body to a magnetic cavity and blocking the solar wind lead to removal of most of the initial magnetic energy stored inside the outer layers of the Moon within minutes along with any enhanced field (up to a factor of 3.5) due to the compression of the solar wind.

Implications for the lunar dynamo hypothesis and magnetization of other impact basins

Theoretical studies of fields generated directly by impacts (1415) have proposed that craters are magnetized while being excavated due to charge separation and electric currents within the vapor, which may induce transient magnetic fields. However, these studies also did not consider the prohibitive effects of the resistive body. Because the predicted induced fields from these studies should be subject to magnetic field dissipation similar to that described for impact amplified fields in the present study, the former existing models may also overestimate the generated paleointensities. This is consistent with the fact that numerous paleomagnetic investigations of impact craters on the Earth have found that impact-heated rocks record the background field and found no evidence of an amplified or locally generated transient field (3437).

The impact-amplified fields hypothesis has been a leading alternative to a core dynamo origin of crustal magnetization in the Moon and other small planetary bodies. Our analysis indicates that any such fields are too weak to explain the strong lunar crustal anomalies and paleointensities of >3.5-Ga-old Apollo samples, supporting the proposal that lunar paleomagnetism is likely largely a record of dynamo action on the Moon. Impact plasmas may still be a viable mechanism for magnetizing some regions of the crust if they form in the presence a preexisting core-dynamo field on the Moon (8), but such an interaction has not yet been studied using MHD simulations. Nonetheless, the recognition that strong and long-duration fields cannot have been formed by impact-generated plasmas in the absence of a dynamo field highlights the need for reconciling the low magnetic field intensities implied by current scaling laws for convective dynamos with high paleointensities observed from the early Apollo samples. Solutions may include energetic dynamos powered by precession (38, 39), a basal silicate magma ocean dynamo (40), or perhaps a core convection dynamo that only transiently generated strong fields (4).


Impact simulations of an Imbrium-sized basin

The iSALE-2D code. The impact simulations were performed using the iSALE-2D shock physics code (23), which is based on the SALE hydrocode solution algorithm (41). To simulate hypervelocity impact processes in solid materials, SALE was modified to include an elastoplastic constitutive model, fragmentation models, various equations of state, and multiple materials (42, 43). Other improvements to iSALE include a modified strength model (44), a porosity compaction model (45), and a dilatancy model (46).

Target and impactor properties. The constitutive model assumed in the simulations is a combination of an analytical equation of state (ANEOS) and a material model (strength and damage of target and projectile rocks and thermal softening and acoustic fluidization as a response to the impact). These parameters and their chosen values for lunar impacts are explained in (29). In the latter, a basalt ANEOS was used for the crust, while here a granite ANEOS was used (47) because it is closer to anorthosite in its properties. Nonetheless, both simulations resulted in a similar basin formation process.

The impact assumed in the simulations is vertical. The cratering efficiency is not sensitive to the exact value of the impact angle and speed as long as the angle is far from the horizontal (>30o) such that a vertical impactor would create a similar basin size as a slightly faster impactor with a slightly lower impact incidence (48). This accommodates for some variation in the impact velocity, which is not exactly known.

Numerical grid setup. The computational resolution is determined via the cell resolution of the projectile, as the projectile effectively initiates the impact calculation. To achieve good accuracy, we used 30 CPPR (cells per projectile radius) in a cylindrical symmetry along the vertical axis. The grid was composed of cells that were 4.5 km by 4.5 km in two dimensions. However, to convert the iSALE-2D grid to the full 3D representation required for driving the MHD simulations, the cells were referred to by their volumes, which grow with the distance from the symmetry axis. The data from iSALE-2D were outputted every 0.1 s and used to derive the input to the MHD calculation.

Impactor size and simulated vapor production. The mass and thermal energy of the vapor emitted from the basin determine the expansion speed and the pressure that the cloud exerts on the solar wind. We performed an iSALE-2D simulation using a dunite impactor with a radius of 60 km striking the Moon at 17 km s−1. These impact parameters would produce an Imbrium-sized basin according to impact simulations performed by (29) and impactor speed estimations from (28). Figure S5 shows the density, temperature, and vertical and horizontal speeds of the vapor from the simulation at 162, 262, 362, and 462 s after the impact.

The impactor properties chosen here differ from those presented in (20), which instead used a projectile of 120-km radius with speed of 18 km s−1. On the basis of updated thermal structure profiles for the Moon, a 120-km impactor could produce almost a South Pole–Aitken–sized basin (49) and therefore cannot be used to create an Imbrium-sized basin. We further justify our choice of impactor size by comparing the total mass of vapor produced by the impact to that appearing in previous studies. The total mass of vapor produced by our simulation can be calculated by integrating over the simulation domain and converting the mass in a 2D slice to the mass obtained over a 3D volume. Each cell in iSALE-2D represents a volume in cylindrical coordinates of size dV = r dφ dr dz for azimuth φ, radius r, and height z and where dφ = 1 rad (23). The total mass over the domain is therefore given by Mtot=ρ(r,z)r dφ dr dz=2πρ(r,z)rdr dz(1)

Using Eq. 1, we integrate the mass over the entire domain at 520 s after impact (the last time step in the iSALE simulation). Regions with density >50 kg m−3 may have more than two phases (e.g., vapor, melt, and solid), in which case iSALE-2D does not give an accurate solution for each of the phases. We therefore discard regions with these densities from all the vapor-related quantities derived from iSALE-2D.

The results of the integration give a total mass of 2.5 × 1019 kg. Although the total vapor mass produced by an impactor with a radius of 120 km in (20) was not reported, we can compare our results to the amount of vaporized mass derived in (7) using analytical calculations of shock vaporization rates from (30). The total vaporized mass obtained by (7) (see table 1 therein) equals 2.77 × 1019 kg for an impactor of radius 68 km hitting the Moon at 15 km s−1, where both impactor and the lunar surface are made of gabbroic anorthosite. We conclude that although the target and impactor materials in (7, 30) and our iSALE-2D simulations differ, the size and velocity of the impactors are very close and produce similar total vaporized mass. The iSALE-2D simulations for the dunite impactor with a 60-km radius are therefore appropriate for testing the impact-amplified fields hypothesis for an Imbrium-sized basin.

MHD simulations

Governing equations and numerical schemes for resistive bodies. The BATS-R-US code (24) enables solving several systems of equations pertaining to MHD, which correspond to different physical approximations of coupled plasmas and electromagnetic fields. Here, we numerically solve the first-order approximation of ideal single-fluid MHD equations outside the body, which control the evolution of mass, momentum, energy, and the magnetic field [equations 6 to 9 in (24)]. Although the single-fluid approximation is used here, the ionized vapor is made of a mixture of elements and the solar wind is made predominantly of ionized hydrogen. Nonetheless, we can correctly capture the overall dynamics by setting the correct mass density of the vapor and the wind.

Inside the body, only magnetic field evolution takes place and there are no plasma dynamics. For the plasma outside the Moon, the induction equation is solved assuming zero resistivity so that only the advective term is solved. Inside the body, there is no advection and only the magnetic diffusion term of the induction equation needs to be solved. However, the magnetic fields inside and outside the body should be coupled. We achieve this by coupling the cells immediately inside and outside the lunar surface so that they pass magnetic field information such that the magnetic field state inside the body is used as a boundary condition to the cells just outside the body, and vice versa. This is done in all cells touching the surface to yield a 3D evolution. The magnetic field diffusion inside the body is solved over a variable-resistivity profile described in section S2. Any numerical scheme that solves for the magnetic field evolution will introduce numerical errors that would violate the condition that the divergence of the magnetic field must be zero. BATS-R-US offers several methods for controlling for these errors; here, we use the hyperbolic cleaning method as described in (24).

The full model equations are solved using a second-order accurate numerical scheme with a Courant stability condition of 0.8. Magnetic field diffusion, which is controlled by a second-order spatial derivative, can severely limit the time step due to the stability condition in an explicit scheme. We therefore use the semi-implicit scheme such that all the MHD operators are solved explicitly, while the magnetic field diffusion term is solved implicitly.

Grid structure. The equations are solved on a 3D spherical grid extending over a range of radial distances r = [0.2, 7] Rm, where Rm is one lunar radius (1737 km). The radial spacing of grid cells varies logarithmically, starting from a cell size of 0.01 Rm at the inner boundary and 0.337 Rm at the outer boundary. The radial cell size at the lunar surface is 0.04 Rm. This way, sufficiently small cells are located where magnetic diffusion and the cloud evolution take place and computational resources needed to resolve the outer regions are reduced.

The grid has a uniform angular spacing of 11° in both latitude and longitude. Adjacent grid points on the surface are approximately 0.09 Rm apart. The spherical cap on the lunar surface approximating the impact basin (from which vapor is emitted according to the iSALE-2D simulations) has a radius of 0.2 Rm. Thus, this grid resolution means that the vapor cloud is emitted from 65 surface grid cells.

We verified that the simulations have sufficient resolution to capture the magnetic field amplification. For this purpose, we repeated the simulation in case 8 while doubling the number of cells in the radial and angular directions and examined the change over time of the maximum field in the crustal layer. We found no appreciable change in the peak field or the timing of its occurrence. The average field over the lunar body volume had only a 1% difference between the two simulations during the peak field time, and no more that 7.5% at any other time. For the purposes of computational expediency and to enable the production of many simulations exploring parameter space, the grid used in all simulations was that of the standard resolution listed above.

Boundary conditions and stopping criterion. For a spherical grid, the boundary condition in the azimuthal and latitudinal directions is periodic. In the radial direction, there are three spherical boundaries: the inner and outer boundaries and an additional intermediate boundary between the lunar body and the surrounding plasma. The inner boundary of the grid (at 0.2 Rm) marks the outer edge of the lunar core. The exclusion of core from the numerical domain significantly improves the computational performance. First, the singularity at the center of the spherical grid is avoided and the smallest grid cell size is limited by omitting the center. Second, the core’s conductivity is orders of magnitudes higher than that of the mantle and crust, having a diffusion time scale of the order of years (Supplementary Materials). This means that the core field does not effectively change over the time scale of vapor expansion (hours). Therefore, magnetic field changes inside the core do not affect the dynamics other than providing an effectively fixed inner field. We therefore apply a fixed inner boundary condition for the magnetic field at the core-mantle boundary, equaling that of the initial IMF vector. With this choice, the time step of the simulation is not limited by the size of the core cells and time steps of the order of 0.03 to 0.1 s are achieved by a typical simulation.

The intermediate boundary condition between the body and the outside plasmas (at r = 1 Rm) is such that the magnetic field is allowed to smoothly vary across it, while three types of plasma boundary conditions are dynamically applied in different regions: (i) Boundary cells exposed to the impinging solar wind fully absorb it, (ii) boundary cells exposed to the dense vapor plasma obey a free-slip boundary condition, and (iii) boundary cells approximating the area of the impact site, described in the next section, emit vapor (until the vapor is turned off, at which point the cells automatically obey any of the previous two boundary conditions).

For the outer spherical boundary (at r = 7 Rm), an inflow/outflow boundary condition is applied such that some cells set the inflow of the solar wind (depending on its flow direction for that simulations), while everywhere else an outflow condition is applied.

The outflow boundary condition requires that all information must flow out of the domain. This means that the simulation is stopped if unphysical inflows are created (this can occur once the vapor blocks out the solar wind from the outflow side of the outer boundary, causing the magnetic field topology to be much different from the initial one). The size of the domain can be increased as to delay this from occurring, effectively allowing the cloud to expand freely without getting close to the boundaries. To avoid unnecessarily large domains (and computationally expensive simulations), we only require that the simulation extends up to the point in time when the total magnetic energy inside the body has reached its peak and started decreasing and the compressed IMF is already being driven away by the vapor. This typically occurs within 50 to 200 s after launching the vapor into the MHD domain. We found that a radial extent of 7 Rm is sufficiently large such that the cloud intersecting the outer boundary does not introduce numerical inflows.

Using impact simulations to drive the MHD boundary conditions at the impact site

The iSALE-2D simulations run for up to 520 s after the impact. Up to the first ~100 s, the vapor is mostly inside the basin (figs. S4 and S5). The shape of the region occupied by vapor is highly irregular, and a significant part of the vapor does not have an upward vertical speed and is trapped in the basin. The basin grows and continues to generate vapor until at ~150 s most of the vapor inside of it has a vertical velocity. At the same time, the walls of the basin grow in the vertical direction, preventing the vapor from expanding in the horizontal direction on the surface of the Moon. The peak flux out of the basin occurs ~250 to 300 s and then decays on a time scale of 250 s (fig. S5). The iSALE-2D simulations were stopped here because we are only interested in the mass flux relevant to calculating the input into the MHD simulations.

We use the results of the impact simulations to derive a time-dependent boundary condition at the impact site inside the MHD simulations. The iSALE-2D density, temperature, and velocity can be directly interpolated into the BATS-R-US grid. However, the irregular shape of the basin and vapor cloud in the initial stages resulted in numerical instabilities in the MHD simulations. In addition, the evolving basin wall cannot be simulated in an MHD code. At the final time step of the iSALE-2D simulation, there is still a lot of vapor mass that is only slowly leaving the basin with low vertical speeds and is prevented from expanding horizontally on the surface by the crater ejecta/wall (fig. S4). Thus, direct coupling would require making some arbitrary choices as to how to treat that vapor and its velocity.

Instead of using the iSALE-2D results directly, we use the rate at which vapor is produced and released from the basin (fig. S5) to create a time-dependent emitting boundary on the surface of the Moon inside the 3D MHD simulations. The emitting region is a circular cap of radius 0.2 Rm, which is the approximate radius of the basin between 200 and 500 s after impact. The MHD simulations are first allowed to achieve a quasi–steady-state solution for a steady wind flow around the resistive body. At this point, the MHD impact simulation is restarted by setting the boundary conditions at the impact region according to the plasma parameters from the iSALE-2D simulation, taken at the time that the vapor starts to leave the basin (around 250 s). Although this occurs around 250 s into the iSALE-2D simulation, the launch of the vapor into the solar wind marks the initiation time of the MHD impact simulations. For simplicity, we set the speed of the vapor to zero at the emitting boundary so as not to impose a specific direction given that the velocity in iSALE-2D is constrained by the evolving impact basin wall. However, the large pressure gradient between the dense vapor and the rarified solar wind quickly accelerates the vapor in the MHD simulation such that it reaches similar speeds high above the basin in both the iSALE-2D and BATS-R-US simulations. In accordance with the temporal evolution of the mass flux out of the basin, the density in the emitting region is decreased exponentially over 250 s after launch, starting from a value of 10−2 g cm−3, corresponding to the density of vertically moving vapor that emerges from the basin starting at around 250 to 260 s after impact (fig. S5). The temperature of the vapor remains fixed at the emitting surface, with a baseline value of 2000 K, which is consistent with the temperature of the vertically moving vapor in the iSALE-2D simulations (this value is reduced in cases 3, 4, and 7 to mimic the effects of vapor-melt friction; see the main text and table S1). After 250 s into the MHD simulation, the boundary condition in the basin region is replaced by the same boundary condition applied elsewhere on the lunar surface (see the previous section).

The emitting boundary approach differs from the vapor expansion simulations in (7), in which the initial condition was defined by having the entire vaporized mass placed inside the domain (above the lunar surface). However, the latter approach, which was implemented in a hydrodynamic simulation, is not suitable for an MHD simulation where the magnetic field is coupled to the plasma. This is because creating an initial region with vapor and no magnetic field violates the divergencelessness property of the magnetic field. Emitting the vapor from the boundary, as was done here, is necessary to avoid this violation; the newly introduced vapor emerges from the boundary and pushes the overlying solar wind and IMF, thus self-consistently creating the cavity. Nonetheless, the shape of the cloud periphery in our simulations (Figs. 1 and 2) is consistent with that produced in (7).

The choice to drive the MHD simulations using only an interval of 250 s out of the iSALE-2D simulation is justified as follows. The vapor that emerges out of the basin during this time is sufficiently massive to ultimately engulf the Moon and sweep the available magnetic flux from inside and outside the Moon to the antipode. From that point on, the creation of more vapor is just replenishing the cloud. The initial field has already been swept up and compressed, while the incoming IMF does not increase the body’s magnetic field because it cannot reach the body (Fig. 1). The IMF is magnetically disconnected from the induced field inside the Moon after cloud convergence. We conclude that the maximum field can only occur while the cloud first expands toward the antipode and later vapor emission can be ignored.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank L. Hood for insightful discussions during the course of this study. Computational 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. Part of this work was carried out using the SWMF (Space Weather Modeling Framework) and BATS-R-US tools developed at University of Michigan Center for Space Environment Modeling (CSEM) and made available through the NASA Community Coordinated Modeling Center (CCMC). We gratefully acknowledge the developers of iSALE-2D, including G. Collins, K. Wünnemann, D. Elbeshausen, T. Davison, B. Ivanov, and J. Melosh. The scientific color maps Vik and Hawaii by F. Crameri were used in this study to prevent visual distortion of the data and exclusion of readers with color-vision deficiencies (50). Funding: We thank the NASA Solar System Workings Program (grant NNX15AL62G), the NASA Solar System Exploration Virtual Institute (grant NNA14AB01A), and the Skoltech Faculty Development Program for support. K.M. research is fully funded by the Australian Research Council. Author contributions: R.O. adapted the BATS-R-US code to perform impact-driven simulations with a resistive body, with assistance from G.T. and K.M. R.O. performed the BATS-R-US simulations. K.M. performed the iSALE-2D simulations. B.P.W. conceived the project, and R.O. and B.P.W. designed the study and overall MHD and paleomagnetism approach and analyzed the results. B.P.W. and Y.S. supervised the project. All authors were involved in writing the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional files can be requested from the authors. The BATS-R-US source code is available to download as part of the SWMF through University of Michigan ( Access to the iSALE code can be requested on a case-by-case basis (

Stay Connected to Science Advances

Navigate This Article