Icarus 286 (2017) 153–192

Contents lists available at ScienceDirect

Icarusjournal homepage: www.elsevier.com/locate/icarus

Magnetization in the South Pole-Aitken basin: Implications for the
lunar dynamo and true polar wander
Michael Nayak a,b,1,∗, Doug Hemingway a,2, Ian Garrick-Bethell a,ca

Department of Earth & Planetary Sciences, University of California Santa Cruz, 1156 High Street, Santa Cruz, CA, 95064, USA
Red Sky Research, LLC, 67 Northland Meadows Dr., Edgewood, NM, 87105, USAc
School of Space Research, Kyung Hee University, Yongin-si, Gyeonggi-do 446-701, Republic of Koreaba r t i c l e

i n f oArticle history:
Received 1 February 2016
Revised 8 September 2016
Accepted 22 September 2016
Available online 14 October 2016
Keywords:
Moon
Magnetic ﬁelds
Impact processesa b s t r a c t
A number of magnetic anomalies are present along the northern edge of the lunar South Pole-Aitken(SPA) basin. A variety of hypotheses for their formation have been proposed, but an in-depth study of
their properties has not been performed. Here we use two different methods to invert for their source
body characteristics: one that completely searches a small parameter space of less than ten uniformstrength dipoles per anomaly, and another that uses grids of hundreds of dipoles with variable magnetization strengths. Both methods assume uniform magnetization directions at each anomaly and with one
exception, produce nearly the same results. We introduce new Monte Carlo methods to quantify errors
in our inversions arising from Gaussian time-dependent changes in the external ﬁeld and the uncertain
geometry of the source bodies. We ﬁnd the errors from uncertainty in source body geometry are almost
always higher. We also ﬁnd a diverse set of magnetization directions around SPA, which we combinewith other physical arguments to conclude that the source bodies were likely magnetized in a dynamoﬁeld. Igneous intrusions are a reasonable explanation (Purucker et al., 2012) for the directional variability, since they could be intruded over different magnetic epochs. However, the directional variability also
implies either surprisingly large amounts of true polar wander or a dynamo not aligned with the lunar
spin axis. We also explore the possibility that true polar wander caused by the SPA impact could allow
iron-rich SPA ejecta to record a diverse set of magnetic ﬁeld directions. Some of this material may havealso become “sesquinary” ejecta and re-impacted across the Moon on 104 –106 year timescales to capture such changes. No completely satisfactory answer emerges, except that the dipole-axis of the lunar
dynamo may have been variable in direction.
Published by Elsevier Inc.

1. Introduction
The South Pole-Aitken (SPA) Basin is the Moon’s largest andoldest well-deﬁned basin (Garrick-Bethell and Zuber, 2009). It isalso the site of the largest grouping of magnetic anomalies on theMoon (Purucker, 2008). Initially discovered by the Apollo 15 subsatellite (Coleman et al., 1972), the origin of these features remains
unknown. Over the last 20 years, three dominant hypotheses fortheir formation have emerged. The ﬁrst is that ejecta from the Im-∗
Corresponding author at: Department of Earth & Planetary Sciences, Universityof California Santa Cruz, 1156 High Street, Santa Cruz, CA, 95064, USA
E-mail addresses: mnayak@ucsc.edu, nayak@redskyresearch.org (M. Nayak).
1
Formerly: NASA Ames Research Center, Planetary Systems Branch (SST), Moffett
Field, CA, 94035, USA, Now at: Maui High Performance Computing Center, Air Force
Research Laboratory (AFRL)
2
Now at: Miller Institute for Basic Research in Science, Department of Earth &Planetary Science, University of California Berkeley, Berkeley, CA, USA

http://dx.doi.org/10.1016/j.icarus.2016.09.038
0019-1035/Published by Elsevier Inc.

brium and Serenitatis basins has accumulated at their antipodes,
which fall within the SPA basin (Hood and Huang, 1991;Lin et al.,
1998; black crosses in Fig. 1). This antipodal ejecta then attains aremanent magnetization from either a dynamo ﬁeld, or a ﬁeld related to the impact event (Hood and Artemieva, 2008). The second
hypothesis is that the anomalies arise from subsurface dikes thatcooled in a dynamo ﬁeld (Purucker et al., 2012). The third is thatiron-rich material derived from the SPA impactor cooled in a dynamo ﬁeld (Wieczorek et al., 2012). Determining which, if any, ofthese hypotheses are true would have implications for understanding the nature and history of the lunar dynamo.
Despite the importance of these anomalies in understanding lunar magnetism, no detailed studies of the source body characteristics have been performed. Purucker et al. (2012) modeled one of
the anomalies within SPA as a series of vertically magnetized dikes,
and estimated their magnetization strength. Global studies of lunar magnetic anomalies have neglected the SPA region (ArkaniHamed and Boutin, 2014; Takahashi et al., 2014), because of the

154

M. Nayak et al. / Icarus 286 (2017) 153–192

6000

5

4000
5
8

−20

8

−25
8

0
5

5

8

2000

5

5

−2000

5

−30

5

8

−35
8

5

5

5

−200

−190

5

5

11
8

5

−180

−170

−160

Elevation (m)

8

−15
5

Latitude (degrees N)

8000
−10

−4000
−6000

−150

−10

12

−15
Area−10Area−2−20−25

Area−5

8

Area−6Area−3Area−1

10

6Area−7−30

4Area−4Area−8−35

2

Total Field Strength (nT)

Latitude (degrees N)

14

Area−9Latitude (degrees N)

−190

−180

−170

−160

−150

−10

12
10

−15
Area−10Area−2

−20
−25

0

8

Area−5

6Area−6Area−3Area−1Area−7

4

Area−82−30
Area−4−35

Total Field Strength (nT)(Mare Ingenii)−200

Area−9(Mare Ingenii)

−200

−190

−180

−170

−160

−150

0

Longitude (degrees E)
Fig. 1. Lunar South Pole-Aitken basin study areas. The SPA rim (thick black line) is from Garrick-Bethell and Zuber (2009). (a) Magnetic contours superimposed over LOLA
topography data (contours are taken from the map in part b). Contours shown are 5, 8 and 11 nT; (b) magnetic ﬁeld map with study areas highlighted (dashed black lines).
Magnetic ﬁeld data are from the Kaguya spacecraft magnetometer measurements taken on day 96 in 2009; at a mean altitude of 39.4 km. Black crosses indicate the locationof the Imbrium (west) and Serenitatis (east) basin antipodes. (c) Magnetic ﬁeld map derived from a spherical harmonic model evaluated at 39.4 km (Purucker and Nicholas,2010). Study areas are highlighted. We refer to areas 1–3 as the northwestern cluster, and areas 5–8 as the eastern stripes.

complicated ﬁeld structure in the region. Using a method they refer to as surface vector mapping, Tsunakawa et al. (2014) calculated
the statistics of the declination of the ﬁeld over the anomalies at
SPA. Based on the distribution of declinations, they found evidence
that the source bodies are horizontally elongated in the east-west
direction. While they also estimated the magnetization direction
and magnetic paleopole from the anomaly centered on the Leibnitz
crater within SPA, they did not estimate the source body characteristics of any of the larger magnetic anomalies that characterize theregion.
Here, we will show that many individual anomalies can be
identiﬁed within SPA. We use software that can readily generateand compare magnetic ﬁeld maps from all available observations,
helping us avoid spurious or poorly deﬁned magnetic anomalies
that may have complicated other efforts to study the region. Comparisons with a spherical harmonic model of the magnetic ﬁeld
(Purucker and Nicholas, 2010) provide a further test of consistency.
We will show that the source body magnetization directions offera key test for their formation hypotheses, and we make a substantial effort to characterize the uncertainty in these directions. The
paper is organized as follows: Section 2 presents our methods, including uncertainty estimation, Section 3 presents our results, andin Section 4, we discuss possible origins of the observed diversity in magnetization directions. We consider impact and dynamo

M. Nayak et al. / Icarus 286 (2017) 153–192

origins for these anomalies, and attempt to explain their directional variability with secular variation, a dynamo not aligned withthe lunar spin axis, and both long and short-timescale true polar
wander. Finally, Section 5 contains our conclusions.
2. Methods
2.1. Data sources
We use magnetometer data from two independent sources: Lunar Prospector (LP-MAG) and SELENE/Kaguya (K-MAG). The Lunar Prospector ﬂuxgate magnetometer measured the vector magnetic ﬁeld at up to 18 Hz and transmitted its measurements at
a reduced resolution of 9 Hz. Level 1B LP-MAG data are obtained from the NASA Planetary Data System (ppi.pds.nasa.gov).
The SELENE/Kaguya spacecraft used a tri-axial ﬂuxgate magnetometer with a sampling rate up to 32 Hz. K-MAG magnetometer data are obtained from the SELENE data archives(l2db.selene.darts.isas.jaxa.jp). The cadence of measurements usedin this study is 0.2 Hz for LP-MAG and 0.25 Hz for K-MAG. Topography data are from the Lunar Orbiter Laser Altimeter (LOLA) (Smithet al., 2010) (pds-geosciences.wustl.edu).
To best capture the Moon’s crustal ﬁeld, all data used for analysis were collected in either the lunar wake or while the Moon
was in the Earth’s magnetotail (hereafter referred to as wake andtail data), avoiding distortions caused by the solar wind noted by
Halekas et al. (2008); Kurata et al. (2005). Tail datasets speciﬁcally exclude epochs during which plasma sheet disturbances werenoted (Halekas et al., 2012). Consecutive orbits are ∼1° in longitude apart. At 0.2 Hz, successive magnetometer measurements are
separated by ∼8 km in the latitudinal direction. All data used inthis study are from the ﬁnal months of Lunar Prospector in 1999
and Kaguya in 2009, when measurements were taken at observation altitudes below 50 km.
2.2. Data processing and anomaly identiﬁcation
Details on the generation of magnetic ﬁeld maps used in thispaper (e.g. Fig. 8) are after Hemingway and Garrick-Bethell (2012),
summarized here for completeness. After subtracting the background ﬁeld (taken to be the mean ﬁeld across each orbit segment
spanning the region of interest), the remaining ﬁelds are assumed
to be due to crustal sources. Data from consecutive orbits are combined and ﬁt to square meshes using Delaunay triangulation, witha grid spacing of 0.25° × 0.25° (7.6 km × 7.6 km equatorial). This
spacing is ﬁner than the spacing between observations (e.g., ∼8 km
latitude, ∼30 km longitude for LP-MAG), ensuring no loss of signal variation during grid generation. Orbit segments that visuallyappear distorted by transient magnetic ﬁeld oscillations, despite
these efforts, are discarded. Gaps in spacecraft coverage, if any, canbe seen in magnetic ﬁeld maps, as the spacecraft measurement locations are shown in all ﬁgures (e.g. Fig. 8, white points).
To analyze the northern SPA region, we divide the strongest
magnetic anomalies into ten study areas (Fig. 1). We choose
anomalies whose Kaguya and LP magnetic ﬁeld maps are consistent with maps from a spherical harmonic model of the ﬁeld
(Purucker and Nicholas, 2010), which can be seen by comparing
Fig. 1b and c. Anomalies that were not consistent with this modelor showed artifacts of external disturbances were not used. Areas
1 and 2 encompass two approximately linear and perpendicularmagnetic features; in their vicinity is area 3, also in proximity to
the Van de Graaff crater (Dyal et al., 1974). Area 1 is very similar
to the area modeled by Purucker et al. (2012). Together, we refer toAreas 1–3 as the northwestern cluster. Area 9 is to the southwestof these areas and is associated with the “swirl” albedo anomalies
at Mare Ingenii (Blewett et al., 2011; Kramer et al., 2011). Area 4

155

is an isolated anomaly to the southeast of areas 1–3. Areas 5–8
are roughly east-west trending linear features found to the east ofthese study areas, and together we refer to them as the eastern
stripes. The anomalies at areas 1 and 8 appear connected in some
maps of the magnetic ﬁeld (Purucker et al. 2012). Finally, we alsodeﬁne area 10, interesting because it falls outside the basin’s outer
topographic rim (Garrick-Bethell and Zuber, 2009).
We divide all data in this paper by the day of observation, suchthat each dataset for a given anomaly is collected at approximately
the same altitude. The datasets may be referenced to the dayon which they were collected, e.g., 1999 day 172 (LP-MAG; readas 99172) or 2009 day 123 (K-MAG; read as 09123). Figs. 2 and3 show an overview of the total magnetic ﬁeld maps created for allobservation days used for all study areas; Appendix A shows details of magnetic inversions for all these areas (see Sections 2.3 and2.4). Supplementary Appendix B shows maps of the total ﬁeld for
all available days that were considered, but not used.
2.3. Inversion algorithm 1: deﬁned dipoles, constant magnetization(DD-CM)To assess the robustness of our results, we use two different algorithms to invert for source body characteristics. Both are regressions to ﬁnd the least squared error. The ﬁrst completely searchesa parameter range using manually placed dipoles as the source
bodies. The regression parameter space varies burial depth in km,
magnetic dipole moment in Am2 , magnetic dip (inclination) in degrees downward from the horizontal, and declination in degrees
clockwise from north. All dipoles in a given study area are constrained to be at the same depth, moment and direction. Source
dipole positions (latitude θ s and longitude φ s ) are placed manually by the user at locations where the magnetic ﬁeld strength isgreatest according to our maps. We refer to this henceforth as thedeﬁned dipoles, constant magnetization (DD-CM) algorithm. As wewill show, the ability to completely search a parameter space at
relatively fast computation speeds gives this algorithm some advantages when uncertainties must be estimated.
The depth is allowed to vary between 1 km and 99 km in stepsof 0.25 km, and the dipole magnetic moment between 1011 Am2
and 1014 Am2 in steps of 2.5 × 1011 Am2 . Depth solutions are further constrained to be no shallower than the deepest negative topography in a given study area, ensuring that all solutions lie beneath the lunar surface; depth is measured against a reference
sphere of 1737.4 km (Smith et al., 2010). The magnetic inclinationis allowed to vary from −90° to + 90° in steps of 1° and declinationfrom −180° to + 180° in steps of 1°.
For n ﬁeld measurements inside a given study area, the difference between the model and the data is computed. We minimize
the root mean square (RMS) total error of the ith east (δ Beast ), north(δ Bnorth ) and radial (δ Bradial ) component errors, δ Btotal , as:(δ Btotal )RMS =
1 δ Bi, east 2 + δ Bi, north 2 + δ Bi, radial 2
n
n

(1)

i=1

This quantity is then ranked in decreasing order of total error
to ﬁnd the best-ﬁt magnetic characteristics.
2.4. Inversion algorithm 2: gridded dipoles, variable magnetization(GD-VM)
We test the validity of our inversion results using an alternative, more objective approach, at the expense of exploring a muchlarger parameter space. Instead of placing dipoles manually, weestablish a 0.25° × 0.25° grid of dipoles covering the entire study
area (Hemingway and Garrick-Bethell, 2012; Nicholas et al., 2007).

156

M. Nayak et al. / Icarus 286 (2017) 153–192

Fig. 2. Overview of all datasets used for SPA areas 1–5. Maps show total ﬁeld strength in units of nanotesla, as a function of latitude and longitude (degrees). The number
in parentheses denotes a study area from Fig. 1. Mean measurement altitudes are shown. Data collection days may be read in the format YY-DDD, e.g., 09123 is 2009 day
123, and so on.

All dipoles are again constrained to be at the same depth, inclination and declination. Unlike the ﬁrst algorithm, however, the magnetic moment is allowed to vary among the dipoles. This createsa large m + 3 parameter space, where m is the number of dipoles;
Table 1 lists the number of gridded dipoles for each study area,
which ranges from 320–1408. The solution is then found via a genetic search algorithm that minimizes the RMS error (for details,
see Hemingway and Garrick-Bethell, 2012). The algorithm evolves
through “generations” to progress toward a better ﬁt (smaller error) to the data. However, the evolutionary nature of the algorithm does not guarantee optimality of the solution. We refer tothis henceforth as the gridded dipoles, variable magnetization (GDVM) algorithm. The grid of best-ﬁt magnetic moments found by
this method, ranging across several orders of magnitude at eachanomaly, can be seen in Supplementary Appendix C, which shows
representative results for each study area.
2.5. Uncertainty estimation for magnetization directions
Error in our regressions arises from: (1) time-variable contributions to the measured ﬁeld due to non-crustal ﬁelds or instrument
noise, (2) the ideal assumption that the source is uniformly magnetized in the same direction, and (3) the simpliﬁed geometry of
our source models, even if the assumption of unidirectional magnetization were completely true.
2.5.1. Time-variable ﬁelds
One method to account for the ﬁrst source of error is to report the set of solutions with an RMS error equal to or less thanthe uncertainty in the magnetic ﬁeld measurements (Parker, 1991).

The range of magnetization directions in this set then deﬁnes theerror in the best-ﬁt direction. The drawback of this method is thatthe use of the ﬁeld’s uncertainty is arbitrary. For example, at different observation altitudes, the magnitude of the RMS error forthe best-ﬁt solution will vary, such that the measurement uncertainty will have a different impact on the size of the set of solutions that represent the error. More importantly, as we show below, the set of allowable solutions suggested by propagating Gaussian measurement noise (or external ﬁeld oscillations) of a given
magnitude through the inversion is often small compared with thesize of the set of solutions contained by an RMS error of equalmagnitude. That is, using the uncertainty in the ﬁeld measurement
to deﬁne the allowable solution set can overestimate the error in
the regression.
To demonstrate this effect, we simulate the effect of Gaussian measurement noise in each of the 45 regressions we perform(all maps shown in Figs. 2 and 3). Before doing so, we ﬁrst estimate the characteristic noise magnitude. For both LP-MAG and
K-MAG, the instrument noise is < 0.1 nT, such that the dominantsource of error arises from ﬂuctuations in the interplanetary magnetic ﬁeld (IMF) rather than the instrument. To estimate a characteristic amplitude of the oscillations in the IMF over the timescaleof observation at one magnetic anomaly, while in the wake ormagnetotail, one would ideally examine the time-dependent oscillations just prior to ﬂying over the anomaly of interest. However,there is no simple way to separate oscillations arising from smallscale crustal ﬁelds and time dependent ﬁelds. As a rough proxy,
we can examine the ﬁeld variations over the Mare Imbrium region
while the Moon is in the wake. Since the Imbrium region is known
to have the lowest amount of crustal magnetism anywhere on theM. Nayak et al. / Icarus 286 (2017) 153–192

157

Fig. 3. Overview of all datasets used for SPA areas 6–10. Maps show total ﬁeld strength in units of nanotesla, as a function of latitude and longitude (degrees). The number
in parentheses denotes a study area from Fig. 1. Mean measurement altitudes are shown. Data collection days may be read in the format YY-DDD, e.g., 09123 is 2009 day
123, and so on.

Moon (Mitchell et al., 2008; Tsunakawa et al., 2014), magnetic ﬁeld
oscillations there are most likely to represent only time dependent
effects. We ﬁnd that the standard deviations of the ﬁeld at Imbrium are 0.07 nT and 0.13 nT for two representative days (Fig. 4).
However, to be conservative, we use an order of magnitude higher
standard deviation of 1 nT as the characteristic value. Takahashi etal. (2014) make a similar (1 nT) assumption about the magnitude
of external ﬁeld contributions.
Next, to simulate the effect of typical time-dependent ﬁelds in
our DD-CM regressions, we added random Gaussian noise with astandard deviation of 1 nT to each ﬁeld component (east, north,radial) for all 45 datasets. An example is shown in Fig. 5. Note that
we are assuming that the noise at each measurement point is uncorrelated with the previous measurement point. Our assumption
is only valid if oscillations in the IMF occur on timescales fasterthan the time between measurements (5 seconds for LP data, and4 seconds for Kaguya data). Future work will determine the dependence of our error estimates on a variety of correlation timescales,
and better determine the typical timescale of time-dependent ﬁeld
oscillations.
For each of the 45 datasets (Figs. 2 and 3), we generated 20
sets of noise-added data. Assuming the best-ﬁt directions follow
a Fisher distribution (the Fisher distribution gives the probabilityof ﬁnding a direction within an angular area from the true mean;

Fisher, 1953), their angular standard deviations and Fisher distribution precision parameters k, were then obtained from each setof 20 simulations. The mean k value across all altitudes within astudy area was then used as the ﬁnal estimate of the uncertainty
from time-dependent ﬁelds for that study area (for example, forarea 1, which uses data from four different altitudes, four different k values were averaged together). The ability to simulate theeffect of noise in this manner represents an advantage of the DDCM method, which has a faster computation time compared to the
GD-VM method.
We ﬁnd that the angular standard deviation in best-ﬁt direction associated with 1 nT noise is always smaller than the set of
directions permitted when using an error threshold of 1 nT. For
example, the best-ﬁt RMS error for area 1, 09096 (Fig. 5) is 2.4 nT,
which is already greater than the ∼1 nT uncertainty we estimate
for the ﬁeld (this inequality is true for nearly all inversions in thisstudy). The set of solutions that are allowed by considering 1 nT
additional RMS error (total residual error 3.4 nT) produces ∼30°of uncertainty in the best-ﬁt magnetization direction (Fig. 6). In
contrast, Fig. 6 also shows that adding 1 nT Gaussian noise to theobservations via the above method only produces an angular standard deviation of 5° in the best-ﬁt direction. The smaller size ofthe effect of simulated Gaussian noise is true in all of our regressions.

158

M. Nayak et al. / Icarus 286 (2017) 153–192

Table 1
Lunar Prospector (LP) and Kaguya (KG) datasets used for inversions, with best-ﬁt parameters found using both the DD-CM and GD-VM algorithms. Mean observational
altitude, data day and site latitude and longitude are shown. GD-VM dipoles column indicates the number of gridded dipoles used for the GD-VM simulation. Estimatedmagnetizations from both algorithms are shown (see text for details). The merged row for each study area shows the best-ﬁt parameters obtained from merged data at all
altitudes using the GD-VM algorithm. Highlighted rows correspond to ﬁgures shown in the main text (Figs. 8–17); GD-VM dipole maps for these datasets are also shown in
Supplementary Appendix C. Best-ﬁt ﬁgures for non-highlighted rows are in Appendix A.
DD-CM algorithm

GD-VM Algorithm

Study Area

Spacecraft

Year

Day

Latitude(deg N)Longitude(deg E)

Alt
(km)

Depth
(km)

Moment(1012 Am2 )

Inc
(deg)

Dec
(deg)

Mag
(A/m)

Dipoles

Depth
(km)

Inc
(deg)

Dec
(deg)

Mag
(A/m)

1

LP-MAG
LP-MAG
K-MAG
K-MAG
MergedLP-MAG
K-MAG
K-MAG
K-MAG
MergedLP-MAG
K-MAG
K-MAG
K-MAG
MergedLP-MAG
LP-MAG
K-MAG
K-MAG
MergedK-MAG
K-MAG
LP-MAG
K-MAG
MergedLP-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG
MergedLP-MAG
LP-MAG
K-MAG
K-MAG
MergedK-MAG
LP-MAG
K-MAG
LP-MAG
K-MAG
MergedLP-MAG
LP-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG
MergedLP-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG
Merged

1999
1999
2009
2009

142
61
96
123

−25
–25
–25
–25

–184
–184
–184
–184

22.4
32.3
39.4
33

25
24
36
37

8.4
6.4
11.0
6.4

13
32
42
20

–123
–116
–132
–114

0.12
0.09
0.06
0.11

1999
2009
2009
2009

33
69
96
123

–20
–20
–20
–20

–188
–188
–188
–188

32.8
44.7
39.4
33

46
37.5
40
40

21.0
26.0
15.0
20.0

36
29
26
32

21
–129
–135
–141

0.06
0.06
0.07
0.07

1999
2009
2009
2009

61
69
96
123

–27
–27
–27
–27

–192
–192
–192
–192

32.3
44.7
39.4
33

47
55
49.5
48

19.0
62.0
18.0
17.0

31
32
49
49

1
–8
4
2

0.08
0.10
0.07
0.06

1999
1999
2009
2009

115
142
96
151

–33
–33
–33
–33

–183
–183
–183
–183

32.3
22.4
39.4
13.7

28
40.5
28
23.5

6.4
12.0
6.2
5.7

–64
–73
–69
–82

–46
–89
–23
86

0.06
0.07
0.06
0.06

2009
2009
1999
2009

96
151
142
123

–19
–19
–19
–19

–165
–165
–165
–165

39.4
13.7
22.4
33

32.5
55
46
38.5

11.0
34.0
21.0
7.9

–33
–24
–29
–56

–16
–12
–23
6

0.06
0.12
0.09
0.06

1999
1999
1999
2009
2009

142
33
115
96
123

–24
–24
–24
–24
–24

–165
–165
–165
–165
–165

22.4
32.8
32.3
39.4
33

46
55
55
46
55

23.0
33.0
30.0
22.0
21.0

47
49
41
49
43

–38
24
–8
31
–14

0.10
0.08
0.08
0.07
0.08

1999
1999
2009
2009

33
115
96
123

–27
–27
–27
–27

–165
–165
–165
–165

32.8
32.3
39.4
33

55
41.5
24.5
37

22.0
12.0
8.0
11.0

35
45
46
40

–160
–180
152
131

0.09
0.07
0.08
0.07

2009
1999
2009
1999
2009

123
61
151
142
96

–31
–31
–31
–31
–31

–164
–164
–164
–164
–164

33
32.3
13.7
22.4
39.4

25.5
46
25
55
28

6.9
21.0
5.9
20.0
7.8

2
1
–17
–17
–1

146
–172
138
–128
–162

0.05
0.03
0.05
0.04
0.03

1999
1999
1999
1999
2009
2009

61
88
115
142
96
123

–36
–36
–36
–36
–36
–36

–198
–198
–198
–198
–198
–198

32.3
32.2
32.3
22.4
39.4
33

29.5
34
25
22.5
26
23.5

26.0
37.0
30.0
17.0
23.0
15.0

1
–12
–6
–2
–5
–5

119
146
139
144
125
135

0.16
0.13
0.14
0.18
0.16
0.19

1999
1999
1999
2009
2009

61
88
142
123
151

–14
–14
–14
–14
–14

–153
–153
–153
–153
–153

32.3
32.2
22.4
33
13.7

36
20.5
21
19
16

25.0
8.3
17.0
15.0
14.0

34
32
29
21
25

–180
166
–175
–172
–168

0.05
0.10
0.06
0.06
0.07

580
580
580
580
580
1408
1408
1408
1408
1408
432
432
432
432
432
320
320
320
320
320
850
850
850
850
850
576
576
576
576
576
576
341
341
341
341
341
1100
1100
1100
1100
1100
1100
672
672
672
672
672
672
672
506
506
506
506
506
506

5
5
26
14
6
21
39
38
26
12
21
27
23
23
11
24
13
24
8
9
24
25
26
28
12
23
31
21
32
34
12
26
25
21
24
14
39
45
13
36
46
13
7
9
7
4
7
5
4
15
12
5
5
4
4.5

18
22
30
20
19
11
5
7
7
8
20
13
27
21
25
–55
–52
–70
–62
–85
–23
–34
–30
–36
–30
50
39
47
31
38
43
28
38
34
34
39
–3
3
0
–6
–3
–1
–4
–11
–8
–5
–8
–3
–7
15
11
14
10
10
13

–128
–128
–126
–116
–123
–71
–54
–52
–53
–59
3
–1
–8
2
–8
–56
–50
–28
–42
–73
–23
–33
–33
–25
–39
–28
53
18
–21
–25
15
–132
159
155
143
165
–131
134
123
–119
–132
–119
145
162
141
153
141
143
146
156
–167
–174
–172
–155
–160

0.70
0.62
0.12
0.24
0.42
0.20
0.10
0.11
0.16
0.14
0.16
0.16
0.13
0.13
0.15
0.09
0.20
0.08
0.30
0.17
0.15
0.17
0.14
0.11
0.14
0.16
0.13
0.16
0.12
0.11
0.14
0.15
0.12
0.14
0.12
0.13
0.08
0.06
0.31
0.10
0.07
0.12
0.52
0.43
0.49
0.95
0.41
0.65
0.58
0.20
0.23
0.56
0.50
0.87
0.47

2

3

4

5

6

7

8

9

10

2.5.2. Non-uniform magnetization direction
The second contribution to error is the non-uniform magnetization directions in the source bodies. However, it is impossible
to make any inferences about the magnetization if we permit theinﬁnite number of possible source magnetizations with mixed directions. Therefore we must at least assume that the source is uniformly magnetized. We attempt to mitigate this effect by selecting
small, well-deﬁned anomalies (Fig. 1). Of course, if non-uniformity
dominates the source of error then these regression results are lessmeaningful.

2.5.3. Source body geometry
Finally, we address what is likely the dominant source of uncertainty: the complex geometry of the magnetic source bodies compared to our simpliﬁed models. We do this in two ways. The ﬁrst
is to compute the best-ﬁt magnetization separately for different
altitude datasets at each anomaly, instead of ﬁtting to a merged
altitude dataset. Data from different altitudes display a diversityof magnetic ﬁeld morphologies, as ﬁeld strength decays with altitude (which ranges in this study from 13–46 km). Hence, each
altitude will lead to different choices of dipole locations for theM. Nayak et al. / Icarus 286 (2017) 153–192

(b)

28
1

26

0.8

24

0.6
0.4

22

0.2

100

50
0

6

0

8

28

2

26

1.5

24

1

22

0.5

0.5

1

1.5

Mean: 0.5 nT, Standard deviation: 0.13 nT

100

50

0

0

2
4
6
8
Longitude (degrees E)0(d) 150

2.5

Total Field Strength (nT)

Latitude (degrees N)(c)

4

Observations

2

Mean: 0.41 nT, Standard deviation: 0.07 nT

150
Observations

1.2

Total Field Strength (nT)

Latitude (degrees N)

(a)

159

0

0.5
1
1.5
Magnetic field strength (nT)Fig. 4. Magnetic ﬁeld maps and histograms of ﬁeld strength at eastern Mare Imbrium in the lunar wake. Panels (a) and (c) are magnetic ﬁeld maps at Imbrium using
LP-MAG observation days 118 and 131 in 1999, respectively. White dots indicate observation locations. The mean observation altitude for day 118 is 31.3 km and 29.8 km forday 131. Panels (b) and (d) are corresponding histograms of the strength of the total magnetic ﬁeld.

−182

−186
4

−23
−24
−25
−26
−27

2
0
−2
−186

−182

−4
6
4
2
0
−2
−4

−23
−24
−25
−26
−27
−186

−186

−23
−24
−25
−26
−27
−186

−182

4
2
0
−2
−4
−6
−8

4
2
0
−2
−4

−23
−24
−25
−26
−27

5

−23
−24
−25
−26
−27

−5
−182

−23
−24
−25
−26
−27

1

0
−5
−182

−182
2
1
0
−1
−182
3
2
1
0
−1

−186
5

−186

2

−186

0

−186

3

−186

−182

−23
−24
−25
−26
−27

−182

−23
−24
−25
−26
−27

−182

4

−182

−23
−24
−25
−26
−27

0
−2
−4
−186

TOTAL (nT)

2

−23
−24
−25
−26
−27

EAST (nT)

4

Recovered model
10
8
6
4
2

−23
−24
−25
−26
−27

NORTH (nT)

6

−186

Latitude (degrees N)

Data + Noise
8

RADIAL (nT)

Data
−23
−24
−25
−26
−27

−182

Longitude (degrees E)
Fig. 5. Error estimation method for time-variable (non-crustal) magnetic ﬁelds. Kaguya magnetic ﬁeld data for area 1, 2009 day 96 (left column), the effect of random
noise with a standard deviation of 1 nT added to that data (middle column) and one of our Monte Carlo trial regressions for the depth, moment and direction values fordata + noise with the DD-CM algorithm (right column). Dipole source geometry is identical to the placement in Fig. 8 (black dots).

DD-CM algorithm; for example, compare Fig. 8 to Figures A-1 to
A-3 (Appendix A). This variability is an advantage, as it allows us
to probe the sensitivity of our results to source geometry, withthe spread in the best-ﬁt magnetization directions representing themeasure of uncertainty (this advantage will be particularly obvious
for area 2, which exhibits bimodal error minima in its error space).
This is analogous to the practice in paleomagnetism of samplinga single rock formation multiple times at different sites (Irving,1964). Because each measured direction at each altitude is independently calculated we can also assign a 95% conﬁdence interval
to the mean direction from all altitudes (Butler, 1998). In practice,

we combine the directions obtained from both the DD-CM and GDVM algorithms to calculate the ﬁnal mean and conﬁdence interval,even though the direction measurements for a given day are notstrictly independent across the two algorithms (see Section 3). In
calculating the ﬁnal mean, we also include one best-ﬁt direction
obtained from the merged data at all altitudes, using the GD-VM
algorithm.
The second way we account for uncertainty due to source bodygeometry is to explore the range in directions that are returnedby displacing the nominal DD-CM algorithm dipoles for all 45
datasets. We calculate the effect of placing dipoles anywhere ran-

160

M. Nayak et al. / Icarus 286 (2017) 153–192

nT
7.5
7.0
6.5
6.0
5.5
5.0
4.5
4.0
3.5
3.0
2.5
Fig. 6. RMS error map for the magnetization direction for area 1, 2009 day 96, centered at 0° inclination and 0° declination (whole-sphere Mollweide projection, in which
the southern hemisphere is positive inclination, same as in Fig. 18). The small white circular contour indicates one angular standard deviation of dispersion, from Monte
Carlo simulations of the addition of 1 nT Gaussian noise to observations (Fig. 5). The larger, outer white dashed contour indicates the minimum error solution plus 1 nT.

Latitude (degrees N)

TOTAL FIELD (nT)

Finally, it is linearly added to the angular 95% conﬁdence interval obtained from the mean of the variable altitude results withina given area (conservatively assuming the Monte Carlo error does
not add in quadrature, since it may be correlated).
The following is a summary of our uncertainty estimation
methods:

-186

-184

-182

Longitude (degrees E)
Fig. 7. Error estimation method for uncertainty in source body geometry. Black
points are the nominal dipole locations for a DD-CM regression at area 1. Circles
represent the allowable locations of dipoles placed randomly in our Monte Carlosimulations and are 0.5° in radius. K-MAG total ﬁeld observations in the lunar wakeon 2009 day 96. Compare to Fig. 8.

domly on a 0.5°-radius circle from their nominal location (Fig. 7).
The value of 0.5° represents the approximate error in longitudinal
uncertainty in the anomaly peak ﬁeld location (orbits are spaced
by approximately 1° for a given constant altitude data set). For
each data set, we use 100 random placements of dipoles at any
location on the deﬁned circles and calculate the 100 best-ﬁt magnetization directions and k values. The k values are then averaged across altitudes within an area to obtain the uncertainty fromsource geometry for that area. Again, this represents an advantageof the DD-CM method compared to the GD-VM method, due tothe ability to modify model dipole placements and a faster computation time. We will see this advantage manifest when exploring ambiguous regression results for area 2. Finally, we note thatfor area 4 we use 0.25°-radius circles, due to the close proximityof the model dipoles there.
The dispersions in best-ﬁt directions obtained from our Monte
Carlo simulations for time-dependent ﬁelds (Section 2.5.1) and
0.5°-radius dipole displacements are then combined using their kvalues, averaged together for each area (Irving, 1964):

1/ktotal = 1/ktime−dependent, mean + 1/ksource−geometry, mean

(2)

An angular standard deviation θ 63 is then estimated by (Irving,1964):

θ63 ≈ 81/k0.5(3)

1. Account for time-dependent external ﬁeld contributions andinstrument noise: Perform Monte Carlo simulations for every
anomaly and altitude with the addition of 1 nT noise, and obtain mean k for the best-ﬁt directions within an anomaly area,
using the DD-CM algorithm (Fig. 5, Table 2).
2. Account for uncertainty in source geometry: Perform Monte
Carlo simulations for every anomaly and altitude while altering
the dipole placement, and obtain mean k for the best-ﬁt directions within an anomaly area, using the DD-CM algorithm (Fig.
7, Table 2).
3. Combine the k values from steps 1 and 2 and obtain the angular standard deviation via Eqs. 2 and 3.
4. Account for uncertainty in source geometry by using variable
altitude data, and hence variable dipole placement. Obtain a
95% conﬁdence interval on the mean direction from the combined best-ﬁt directions from the DD-CM and GD-VM algorithms (Table 3).
5. Linearly add the angular dispersion from step 3 to the angular
95% conﬁdence interval from step 4 to obtain the total angular
uncertainty estimate in the best-ﬁt direction.
3. Results
3.1. Regression results
Minimum error magnetization directions are found for each Lunar Prospector and Kaguya dataset, using both the DD-CM (Section2.3) and the GD-VM (Section 2.4) algorithms. Table 1 compiles thevalues of these best-ﬁt directions. The best-ﬁt models are shownin Figs. 8 –17 for study areas 1–10. The ﬁgures show 10 representative examples, one for each study area; Appendix A shows all
results for every dataset listed in Table 1. Angular standard deviations and k for the noise-added and dipole displacement simulations are shown in Table 2 for all datasets. RMS error maps for
all datasets (similar to Fig. 6) are compiled in Figs. 18 and 19;

M. Nayak et al. / Icarus 286 (2017) 153–192

161

Table 2
Uncertainty estimates for regression results from Monte Carlo simulations. Minimum RMS error recovered using the DD-CM and GD-VM
algorithms for all areas and altitudes. Fisher precision parameter k, angular standard deviation θ 63 (≈ 81/k1/2 ), obtained from Monte Carlosimulations for the effects of time-dependent ﬁelds (“noise”) and 0.5°-radius displaced dipoles, for the DD-CM algorithm. Highlighted rows
correspond to ﬁgures shown in the main text (Fig. 8–17).
Study Area

1

2

3

4

5

6

7

8

9

10

Spacecraft

LP-MAG
LP-MAG
K-MAG
K-MAG
LP-MAG
K-MAG
K-MAG
K-MAG
LP-MAG
K-MAG
K-MAG
K-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG
K-MAG
K-MAG
LP-MAG
K-MAG
LP-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG
K-MAG
LP-MAG
K-MAG
LP-MAG
K-MAG
LP-MAG
LP-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG
LP-MAG
LP-MAG
LP-MAG
K-MAG
K-MAG

Year

1999
1999
2009
2009
1999
2009
2009
2009
1999
2009
2009
2009
1999
1999
2009
2009
2009
2009
1999
2009
1999
1999
1999
2009
2009
1999
1999
2009
2009
2009
1999
2009
1999
2009
1999
1999
1999
1999
2009
2009
1999
1999
1999
2009
2009

Day

142
61
96
123
33
69
96
123
61
69
96
123
115
142
96
151
96
151
142
123
142
33
115
96
123
33
115
96
123
123
61
151
142
96
61
88
115
142
96
123
61
88
142
123
151

Min. RMS error

Min. RMS error

Noise simulation

Displaced dipole simulation(DD-CM)(nT)(GD-VM)(nT)

θ 63

Precision
parameter kθ 63Precision
parameter k

7.4
5.6
2.4
6.3
5.9
3.7
4.0
5.7
4.1
3.1
1.7
2.8
1.5
4.4
1.1
8.1
2.0
9.1
4.6
2.2
5.7
3.8
3.3
2.2
2.7
3.1
2.2
1.5
2.4
2.6
2.5
7.5
4.4
1.6
6.3
5.6
5.7
9.0
4.0
4.5
3.0
3.2
4.7
1.9
8.0

3.9
2.2
1.1
2.3
3.1
2.3
2.2
3.2
2.0
1.6
0.9
1.3
1.1
2.0
0.8
3.7
1.0
4.0
1.9
1.1
3.1
1.6
2.0
1.1
1.5
1.6
1.3
0.9
1.4
1.7
1.5
3.9
2.6
1.1
3.2
2.3
2.1
3.5
1.5
1.5
1.5
1.7
2.3
1.3
4.2

1.6
2.1
5.0
2.9
1.5
3.3
2.1
1.7
3.1
3.1
4.3
5.3
6.4
4.2
7.2
1.9
5.3
1.3
2.0
5.2
2.4
4.8
4.8
7.8
4.6
5.0
6.0
8.9
5.2
4.4
6.0
2.4
3.7
7.2
1.4
1.8
1.9
1.0
2.8
1.8
3.0
7.9
2.2
3.9
1.2

2526
1474
264
774
3120
620
1542
2318
676
674
353
234
159
375
127
1761
232
3959
1648
243
1150
288
283
109
308
262
184
83
244
347
184
1121
476
128
3496
1999
1837
6997
853
2004
730
104
1315
431
4571

7.4
7.3
4.2
3.8
4.5
5.7
44.1
14.3
5.7
7.2
6.1
6.4
4.7
5.3
5.5
8.6
3.7
3.6
3.3
4.1
15.1
11.5
5.3
6.9
4.4
4.6
6.3
8.3
6.7
28.0
3.9
36.6
3.2
5.3
4.3
8.1
9.6
8.5
8.1
5.4
7.7
11.9
10.2
7.3
15.5

120
124
381
446
332
203
3
32
202
128
177
160
300
232
219
89
472
501
601
383
29
50
234
140
343
311
164
95
147
8
443
5
625
235
352
99
71
91
101
228
112
46
63
124
27

Table 3
Uncertainty estimates derived from using different altitude datasets. The paleolatitude, paleolongitude, value of k, and the 95% conﬁdence circle (∝95 ) are obtained from the combined set ofbest-ﬁt directions from the DD-CM and GD-VM algorithms, across all altitudes. The ∝95 value isadded to the angular uncertainty estimated by Monte Carlo simulations (equations 2 and 3, derived from mean Table 2 values) to obtain the ﬁnal circular error, also shown. The distorted version of this ﬁnal error circle at the paleopole is also shown (Figs. 21 and 22): the semi-axis of theellipse along the great circle path from site to pole (dp) and semi-axis of the ellipse perpendicular
to that path (dm).
Area

1
2
3
4
5
6
7
8
9
10

Altitude averaged
k

α−95

58.5
3.3
38.7
27.8
33.0
14.5
11.2
3.2
50.0
25.4

6.8
33.9
8.4
9.9
9.1
12.4
16.1
30.9
5.9
9.2

Paleolatitude

Paleolongitude

Final

dp

dm

7.0
24.0
9.6
24.5
8.6
16.6
18.2
18.7
6.3
10.3

13.1
44.3
17.4
28.4
15.1
25.7
29.9
37.3
12.6
19.8

circular error
–34.8
8.5
46.9
51.8
68.1
38.0
–78.5
–54.2
–37.0
–83.3

84.2
100.4
165.1
39.1
104.5
–165.1
–62.0
37.8
–67.6
82.9

12.2
41.0
15.7
16.4
13.3
19.9
24.5
37.1
12.6
18.9

162

M. Nayak et al. / Icarus 286 (2017) 153–192

Fig. 8. Area 1 best-ﬁt results for both inversion algorithms. (Left column) K-MAG observations in the lunar wake in 2009, day 96, at 38.5 km, compared to (middle column) amodel using several manually placed dipoles of equal magnetization (DD-CM), and (right column) a grid-based model using many dipoles of variable magnetization (GD-VM).
Source dipoles in the DD-CM algorithm are approximately located at maxima in the observed total ﬁeld (black dots). White dots in the east panel are locations of spacecraft
magnetometer measurements.

Fig. 9. Area 2 best-ﬁt results with both inversion algorithms. Figure details are as in Fig. 8, except that magnetometer measurements are taken from K-MAG, 2009 day 123,
altitude 32.5 km.

these illustrate the difference between using an arbitrary uncertainty threshold and the Monte Carlo methods employed in thiswork.
Fig. 20 compiles the values of the best-ﬁt directions from the
DD-CM and GD-VM algorithms. For the DD-CM algorithm, we show
one standard deviation of dispersion from 1 nT noise (Fig. 20a) andone standard deviation of dispersion from displacing the model
dipoles by 0.5° (Fig. 20b) (see Methods). Overall, we ﬁnd the error

from the uncertainty in source geometry (Fig. 20b) is usually larger
than the error from the effects of time-dependent ﬁelds (Fig. 20a).
Interestingly, we ﬁnd source geometry errors at areas 2 and8 are substantially larger than for the other anomalies. The large
error in area 2 arises from the existence of two nearly equal
local minima in its DD-CM error map, which is not seen at anyother anomaly (see Fig. 18). The two minima are approximately
110° apart. Small displacements of the dipoles from their nom-

M. Nayak et al. / Icarus 286 (2017) 153–192

163

Fig. 10. Area 3 best-ﬁt results with both inversion algorithms. Figure details are as in Fig. 8, except the altitude here is 37.6 km.

Fig. 11. Area 4 best-ﬁt results with both inversion algorithms. Figure details are as in Fig. 8, except the altitude here is 38.3 km.

inal locations ﬂip the best-ﬁt solution into the other minimum,
producing its large uncertainty ellipses, particularly for day 09096(Fig. 20b). One of our four best-ﬁt solutions (day 99033) falls
in this secondary minimum. Because area 1 and 2 have slightoverlap (Fig. 1), we tested if the double minima could be a resultof area 1 data inﬂuencing the regression at area 2. We ran the
DD-CM algorithm against data for area 2, but this time excluded

spacecraft measurements over area 1. The results are not significantly different; the two minima are still seen for area 2. Using
the GD-VM algorithm, one ﬁnds the best-ﬁt directions for area 2
are approximately between the two minima found by the DD-CM
algorithm (Fig. 20c), ∼60° from either one. This difference between
the two algorithms is the largest for any area. Further work will
be required to better understand the nature of the error space at

164

M. Nayak et al. / Icarus 286 (2017) 153–192

Fig. 12. Area 5 best-ﬁt results with both inversion algorithms. Figure details are as in Fig. 8, except that magnetometer measurements are taken from K-MAG, 2009 day 123,
altitude 33.3 km.

Fig. 13. Area 6 best-ﬁt results from both inversion algorithms. Figure details are as in Fig. 8, except that magnetometer measurements are taken from K-MAG, 2009 day 123,
altitude 32.4 km.

area 2, and its source body characteristics; we will investigate, forexample, oppositely magnetized blocks of magnetization in thisregion (cf. Parker, 1988).
For area 8, the large error arises from the ﬂatness of the error
space, instead of the existence of multiple minima (Fig. 18). In
particular, the declination of this nearly equator-pointing magnetization vector is poorly constrained, with results spanning a rangeof ∼100° of arc (a result found to be true using either the DD-CM
or GD-VM algorithm). Again, without the error maps provided by
the DD-CM algorithm, we would not have been able to isolate
and understand these different sources of uncertainty in areas 2
and 8.

Next, we calculate the mean magnetization strength for all
study areas using the best-ﬁt depths and magnetic moments recovered by both the DD-CM and GD-VM algorithms (Table 1). For
this calculation, the area of the magnetic anomalies is taken as thearea across which the magnetic ﬁeld B ≥ Bmax /4, where Bmax is thepeak magnetic ﬁeld in a given study area. We choose to approximate the thickness of the magnetic material as twice the average of the best-ﬁt magnetic source depths. Most values cluster in
the 0.1–0.3 A/m range found by Purucker et al. (2012), who studied a region similar to area 1. A few isolated datasets at the lowest
altitudes show higher magnetizations, but with the exception ofarea 9, average magnetizations for a particular study area across all

M. Nayak et al. / Icarus 286 (2017) 153–192

165

Fig. 14. Area 7 best-ﬁt results from both inversion algorithms. Figure details are as in Fig. 8, except that magnetometer measurements are taken from K-MAG, 2009 day 123,
altitude 31.7 km.

Fig. 15. Area 8 best-ﬁt results from both inversion algorithms. Figure details are as in Fig. 8, except that magnetometer measurements are taken from K-MAG, 2009 day 123,
altitude 31.4 km.

datasets are close to the 0.1–0.4 A/m range suggested by Wieczorek
and Weiss (2010) for mare basalts and maﬁc impact melt breccias.
Even for area 9, only the GD-VM average magnetization is outside
this range.
In most cases the GD-VM algorithm returns higher values formagnetizations than the DD-CM algorithm. This is a result of theshallower depths returned by the GD-VM algorithm, which arise
because of the sheet-like nature of the magnetic source, comparedto the dipolar model in the DD-CM algorithm. In general the depthof the magnetic source bodies is not well constrained from ourresults, due to the possibility of non-uniquely trading depth with
moment.

Finally, best-ﬁt paleopoles for all datasets are calculated fromthe best-ﬁt magnetization directions (Butler, 1998) (Fig. 21,
Table 3). The uncertainties derived from the methods described in
the Methods section yield a circular error ellipse for each area that
becomes distorted when its paleopole is calculated. Hence we calculate the two different semi-axes of the ellipse along the great
circle path from the site to pole (dp) and the semi-axes of the ellipse perpendicular to that path (dm) (Butler, 1998). For area 2, we
also show the paleopoles from the two local error minima foundin the DD-CM method, particularly because one falls very close to
the paleopole for area 1 (within its uncertainty), while the otherpaleopoles from areas 3 and 6 (within the uncertainty of 6).

166

M. Nayak et al. / Icarus 286 (2017) 153–192

Fig. 16. Area 9 best-ﬁt results from both inversion algorithms. Figure details are as in Fig. 8, except that magnetometer measurements are taken from K-MAG, 2009 day 123,
altitude 28.3 km.

Fig. 17. Area 10 best-ﬁt results from both inversion algorithms. Figure details are as in Fig. 8, except that magnetometer measurements are taken from K-MAG, 2009 day
123, altitude 34.2 km.

M. Nayak et al. / Icarus 286 (2017) 153–192

167

Fig. 18. RMS error map (nT) for the magnetization direction for all datasets for SPA study areas 1–5, obtained using the DD-CM algorithm. Southern hemisphere is positive
inclination; the whole-sphere Mollweide projection is centered at 0° inclination and 0° declination. The solid white contour indicates one angular standard deviation ofdispersion from the best-ﬁt solution (red dots, Table 1), from Monte Carlo simulations of the addition of 1 nT Gaussian noise to observations, using the DD-CM algorithm.
The larger, dashed white contour indicates the minimum error solution plus 1 nT, a measure of uncertainty using an arbitrary threshold equal to the measurement uncertainty(not used in our ﬁnal analysis). The error from Monte Carlo simulations of the effects of displacing the nominal source dipole locations by 0.5° is not shown, and nor isthe 95% conﬁdence interval from averaging different altitude data sets. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web
version of this article.)

3.2. Comparison with other magnetic paleopoles
In Fig. 22b-d we show paleopoles from isolated magnetic anomalies studied by Takahashi et al. (2014) andArkani-Hamed and Boutin (2014), as well as paleopoles inferredfrom the magnetization of the crust at craters (Arkani-Hamed
and Boutin, 2014). The magnetized crust paleopoles from ArkaniHamed and Boutin (2014) are the means of the values listed intheir Table 2. None of the anomalies studied by these groups were
inside SPA. To compare our results with these three datasets, wereverse all of our paleopoles into the same southern hemisphere(Fig. 22a), as in Takahashi et al. (2014). Here we have dropped
the paleopole found from the mean of all data sets at area 2, andonly show the paleopoles from the two local error minima in the
DD-CM algorithm.
Takahashi et al. (2014) found paleopoles that cluster into twogroups, one near the present pole, and another at mid-latitudes.
We ﬁnd a wider dispersion in the distribution of our paleopoles,
but we do ﬁnd some clustering at the present pole (areas 5, 7 and10), as well as some near their mid-latitude cluster (areas 3, 6,
and 8). The remaining poles are not easily assigned to either ofthese clusters. However, overall, the paleopoles we ﬁnd do seem
to avoid longitudes on the farside (the bottom half of the sphere in
Fig. 22a), and latitudes < 30°. None of the paleopoles from ArkaniHamed and Boutin (2014) show obvious correlation with any of theclustering found here or in Takahashi et al. (2014).
Finally, we note that many of the paleopoles are substantiallyseparated by their error ellipses, such that it is unlikely that only

one or two paleopoles could be used to explain all of the anomalies, without severely affecting their RMS error values. This canalso be visualized by examining the error spaces in Figs. 18 and19; because most anomalies are at similar latitudes and longitudes(most within ∼30°), the diversity in paleopole locations is mostly
determined by the diversity in magnetization directions. Hence,the diversity in RMS error minima in Figs. 18/19 also graphicallyillustrates the diversity in paleopole locations.
In the next section, we address possible origins for the dispersion seen in our study, and test some of the hypotheses foranomaly formation.
4. Discussion
The wide variation in magnetic paleopoles derived from SPA’s
magnetic anomalies presents a puzzle. In Sections 4.1–4.4, we discuss the implications of the diverse paleopole locations for theformation of these anomalies, the history of true polar wander, andthe nature of the lunar dynamo.
4.1. Magnetization by impact processes
A long-standing hypothesis is that the strongest lunar magnetic
anomalies are genetically related to the antipodes of the Imbrium,
Orientale, Serenitatis and Crisium basins (Hood and Williams,
1989; Hood et al., 2001; Lin et al., 1998). Compression and ampliﬁcation of the IMF by impact-produced plasma may be strongest
at the basin antipode, where impact ejecta may also preferentially168

M. Nayak et al. / Icarus 286 (2017) 153–192

Fig. 19. RMS error map (nT) for the magnetization direction for all datasets for SPA study areas 6–10, obtained using the DD-CM algorithm. Description is as in Fig. 18.

collect (Hood and Artemieva, 2008). The antipodes of the Imbrium
and Serenitatis basins are close to our study areas (Fig. 1), suggesting this process may be responsible for forming the anomalies wehave examined here. If true, antipodal ejecta may become magnetized via either thermo-remanent magnetization (TRM) or shockremanent magnetization (SRM). Below we assess these two possibilities.
4.1.1. TRM in an impact-produced ﬁeld
If the magnetization were produced by a TRM, hot ejecta
would cool in the presence of transient IMF-ampliﬁed ﬁelds thatwould last, at most, for one day. Assuming a thermal diffusivityof 10−6 m2 /s, the thermal cooling timescale for one day is <1 m.
Therefore, TRM would be restricted to <1 m of material. The magnetic moments we ﬁnd across all SPA study areas range between
1013 Am2 and 1014 Am2 . Using source body horizontal extents from
Fig. 1 (black dashed boxes), depths of ∼1 m lead to high TRM magnetizations ranging between 103 and 104 A/m. These values are 3–
4 orders of magnitude higher than samples recovered by Apollo
missions. Additionally, the top ∼1 m would have been completely
overturned and demagnetized in the time since the antipodal impact ∼4 billion years ago (Arnold, 1975). Therefore, TRM from

M. Nayak et al. / Icarus 286 (2017) 153–192

impact-related ﬁelds is implausible, in agreement with Hood andArtemieva (2008).
4.1.2. SRM in an impact-produced ﬁeld
Alternately, impact shock pressures from the deposition ofejecta at the antipode may create SRM in the ejecta deposit. However, to allow unidirectional SRM to be stably imparted to a largescale (e.g. > 40 0 0 km2 for areas 1 and 2; > 30,0 0 0 km2 for areas
5–8) geologic formation, the shock waves must pass through therock when the entire region is at rest. Shock waves begin propagating through the rock instantly upon impact, yet at that moment
the ejecta is still traveling at 2.0–2.4 km/s (Hood and Artemieva,
2008). Therefore, the ejecta at the antipode will change orientation
as it comes to rest, after any SRM has formed, such that the total
remanence of the ejecta deposit will be randomized and nulled.
Hood and Artemieva (2008) also rule out SRM in the ejecta, using
a different line of reasoning.
It is possible that the underlying rock could be shocked by
the impacting ejecta and acquire a unidirectional SRM (Hood andArtemieva, 2008). However, we suggest that there are two problems with this hypothesis. Firstly, one would expect nearly vertical local magnetization if the SRM producing ﬁeld was compressedby solar wind plasma (see, for example, ﬁeld lines in Figs. 8–9 ofHood and Artemieva, 2008). Instead, only one magnetic anomaly,
area 4, has a nearly vertical magnetization, and the rest show apreference for low inclinations, if anything (the rest are at least> 33° away from the vertical, Fig. 20).
Secondly, deposition of impact ejecta would take place overa short period of time, over which the ambient ﬁeld direction
is likely to be nearly constant. There are only two basin antipodes within SPA, which implies that the magnetic anomalies
should have one of two magnetization directions. Roughly, areas 1–
3 (northwest cluster), 4, and 9 are antipodal to Imbrium, while areas 5–8 (eastern stripes), and 10 are antipodal to Serenitatis. First,
we consider the ﬁve anomalies at the Imbrium antipode. Here weﬁnd that areas 1 and 3 have magnetization vectors 109° apart(Fig. 21). The magnetization at area 4 is 95° from that at area 1 and
70° from area 3. The magnetization at area 9 (Mare Ingenii) is 104°
from area 1, 140° from area 3, and 128° from area 4. All of theseseparations are well outside the error ellipses. The only two clusters are the directions of areas 1 and 2, and possibly area 3 withthe second error minimum obtained at area 2. In sum, there areat least four widely separated magnetization directions near the
Imbrium antipode, all well separated by their error ellipses. It is
plausible the similar directions at areas 1 and 2 suggest this pair
was magnetized by deposition of impact ejecta, but it would not
explain the magnetization of the other areas.
At the ﬁve anomalies near the Serenitatis antipode, areas 7, 8,
and 10 have similar directions, but 5 and 6 are well separated fromeach other and from 7, 8, and 10 (Fig. 21). Therefore, we conclude
again that it is not likely that the Serenitatis impact is responsible
for magnetizing all of these anomalies.
In sum, the diversity of directions argues against the SRM hypothesis, or at least allows only a subset of the geographically clustered anomalies to be due to SRM. Similar arguments can be applied to ruling out SRM from surface seismic waves. Considering
all of the observations above, we conclude that at least some of
the magnetic source bodies in SPA were magnetized in a lunar dynamo ﬁeld, rather than a ﬁeld associated with impact events.
4.2. Magnetized South Pole-Aitken basin ejecta or volcanic bodies?
Wieczorek et al. (2012) proposed that material from the SPA
impactor might be the source of many of the SPA basin’s magnetic anomalies, and even other anomalies across the Moon. Under this hypothesis, hot iron-enriched material from the SPA im-

169

pactor acquired a TRM in a dynamo ﬁeld. However, if SPA ejecta
landed hot and cooled in a dynamo ﬁeld, the resulting anomalies should all have the same magnetization direction (but see
Section 4.3.4). Instead, the diverse magnetization directions (seeSection 4.1.2, above) suggest that they were magnetized at different epochs.
The last remaining viable hypothesis is that magnetic anomalies in the northern SPA basin formed as a result of magnetized
sub-surface dikes (Purucker et al., 2012). Dikes forming over longtime periods would permit different magnetization directions during different magnetic epochs. The cluster of three paleopoles close
to the present pole (areas 5, 7 and 10, Fig. 22a) would suggest
a traditional axial-aligned dynamo magnetized these dikes when
the Moon was in its present orientation. However, the diversity
in paleopoles seen is still enigmatic (see Section 4.3). Further, if
the dike hypothesis is true, this implies that the dikes near SPA
are special in some way, since the nearside of the Moon, coveredmuch more extensively by volcanism (and presumably associatedwith subsurface dikes), shows no magnetic structures like those at
SPA. Andrews-Hanna et al. (2014) reported linear gravity-gradientanomalies that may be dikes, but so far, no obvious correlation between these structures and magnetic anomalies have been found.
We note that Gong and Wieczorek (2016) ﬁnd a correlation between magnetization and gravity anomalies in some locations.
4.3. Possible explanations for diverse magnetic ﬁeld directions
The above arguments suggest the magnetizing ﬁelds for manyof the SPA anomalies arose from a dynamo. However, they do notoffer an obvious explanation for the diversity in ﬁeld directions.
Due to the small size of the lunar core, and the rapid decay ofmagnetic quadrupole and higher terms as a function of distance,
it is likely that the ancient lunar dynamo was dominantly dipolarat the surface (Weiss and Tikoo, 2014). If the dipole was alignedwith the Moon’s spin axis, the magnetization directions contain information about the Moon’s paleopole. The diverse paleopoles inFig. 21 seem to imply large amounts of true polar wander(Goldreich and Toomre, 1969; Runcorn, 1983). While the diversity
in Fig. 21 is surprising, in comparison with the results of (ArkaniHamed and Boutin, 2014; Takahashi et al., 2014) (Fig. 22a), the diversity of paleopoles seen in our ﬁndings is still greater. Below we
consider some possible explanations for the diverse paleopole locations.
4.3.1. Non-axially aligned dipoles and impact-induced dynamos
Currently, our understanding of the nature of the lunar dynamo
is limited. Paleomagnetic studies favor a dynamo that existed fromapproximately 4.2–3.3 Ga ago. Mechanisms for sustaining such along-lived core dynamo are uncertain, with recent proposals for
dynamos driven by mechanical stirring from impacts and precession (Dwyer et al., 2011; Le Bars et al., 2011; Weiss and Tikoo,2014). However, it is not known if these dynamos produce thesame ﬁeld organization as the Earth’s dynamo. It may be possible that these exotic dynamos exhibit more variable dipole axis directions, which could explain the diversity of paleopoles seen. Unfortunately, more work modeling small dynamos like that of theMoon is needed before evaluating this hypothesis further.
4.3.2. Secular variation
Another possible origin for the diverse paleopole locations isgeomagnetic secular variation of a dynamo that is on averagealigned with the lunar spin axis. The Earth’s spin axis is presently
11° away from its magnetic dipole axis, and it is plausible that the
Moon’s ancient dipole axis was also not exactly aligned with its
spin axis. It is also possible that secular drift might be larger on abody with a small core such as the Moon. However, some of the170

M. Nayak et al. / Icarus 286 (2017) 153–192

Fig. 20. Best-ﬁt magnetization directions for study areas deﬁned in Fig. 1, for all datasets in Figs. 2 and 3. Directions derived from different altitudes and different spacecraftat the same study area are represented by the same color. Positive inclinations are in the lower hemisphere (ﬁlled circles). (A & B) Results for the DD-CM algorithm.
Ellipses in A represent one angular standard deviation of dispersion from 1 nT noise simulations. Ellipses in B represent one angular standard deviation of dispersion fromsimulations of source dipoles displaced from their nominal positions (0.5°), (C) Results from the GD-VM algorithm, including results from regressions to merged data at all
altitudes shown in Fig. 2 and Fig. 3, within a given study area (stars).

paleopole locations at SPA are over ∼75° apart (after accounting
for the possibility of reversals, Fig. 22a), and over the timescale ofmagnetic anomaly formation, the mean dipole orientation mightaverage to be aligned with the spin axis. Presently, we have little
information available about the lunar dynamo to further evaluate
the role of secular variation in explaining the diversity of paleopole
directions.

4.3.3. True polar wander over long timescales
The large amount of polar wander implied by the paleopolesis diﬃcult to reconcile with other geophysical constraints for theorientation history of the Moon. Currently there are two comprehensive studies that derive estimates of the degree of polar wander on the Moon due to long-term changes in the Moon’s moments of inertia (Garrick-Bethell, 2016). The ﬁrst uses the shapeM. Nayak et al. / Icarus 286 (2017) 153–192

171

Fig. 21. SPA magnetic paleopoles. Paleopoles are the mean of the combined results from the DD-CM and GD-VM algorithms, from all directions (altitudes) shown in Fig. 20.
Ellipses represent the 95% conﬁdence interval obtained when calculating the mean direction from the combined DD-CM and GD-VM results, plus the dispersion from Monte
Carlo simulations of time-dependent ﬁelds and dipoles displaced from their nominal locations (see Methods). Error ellipses are not calculated for the two different minimum
error solutions for area 2 from the DD-CM algorithm.

and gravity of the Moon exterior to large basins to establish theearliest orientation of the Moon (Garrick-Bethell et al., 2014). The
authors ﬁnd the present lunar pole has changed by ∼36° fromits earliest axis, but it does not coincide with any of the paleopole clusters we ﬁnd (Fig. 22a). The second study uses polarhydrogen deposits as a constraint on the history of polar wander, and infers that up to ∼10° of paleolatitude change has occurred (Siegler et al., 2016). However, they cannot readily produce the large paleopole changes implied by the anomalies studied here. In summary, the limited number of available geophysicalmodels that estimate lunar polar wander cannot produce the diversity and large magnitude of paleopole changes required to explain our observations. Interestingly, we ﬁnd that all of our paleopoles are > 35° from the present equator (Fig. 22a, accounting forreversals), which may be due to the diﬃculty in producing thelarge required changes in the Moon’s moments of inertia.

4.3.4. True polar wander due to SPA formation
Very large initial changes in the Moon’s moments of inertia due
to SPA’s crater might have produced large amounts of polar wander. Eventually, these changes must have subsided over millions
to billions of years, since the gravity signature observed today ismuted (Zuber et al., 2013). If hot material was cooling in the presence of a dynamo ﬁeld throughout the Moon’s reorientation, its
magnetization could in principle capture multiple lunar orientations, thereby producing the diversity of paleopole locations. Someof this hot material could be iron-rich material from the SPA impactor (Wieczorek et al., 2012). Following impact, the free precession damping time of the Moon is ∼2 × 105 years (Peale, 1976;
Williams et al., 2001), assuming a semimajor axis of 30 Earth radii,
dissipation quality factor Q = 50 and degree-2 love number k2 = 0.1(which we assume to be representative of the Moon when SPA
formed). The length scale for cooling in ∼2 × 105 years is ∼3 km(assuming a thermal diffusivity of 10−6 m2 s-1 ), which is small
compared to the scale of the anomalies we observe, but perhaps
not so much as to preclude recording a measurable magnetization.

However, using a simple model for reorientation, we ﬁnd thatdensity anomalies produced by SPA’s crater do not produce paleopoles that overlap with those of its magnetic anomalies. In
our model, we replace SPA’s gravity potential inside the outer topography rim with values between −3 times the maximum, and+ 2 times the maximum of the present day potential, in increments of 0.25 times the maximum potential. The inertia tensor
and paleopoles were then calculated for each of these cases using the resulting globally calculated degree-2 spherical harmonics(cf. Garrick-Bethell et al., 2014). These positive and negative gravity
values represent very different models of SPA’s effects, but illustrate the range of paleopoles that are possible. We ﬁnd the maximum extent of true polar wander from SPA’s formation passes
close to the paleopoles from areas 7 and 10 (Fig. 22a, the negative
anomaly path approaches the center of SPA as the anomaly size
grows, as expected). However, these anomalies are already close to
the present pole, and do not need to be explained by polar wander.
We do ﬁnd that the limits of polar wander approach areas 4 and 8,
such that this process could plausibly explain some of these magnetic paleopoles. Unfortunately, the paleopoles from areas 1, 2, 3, 5,
and 9, and to a lesser extent area 6, do not come close to the paleopoles produced by SPA. Thus, cooling of SPA ejecta deposits overshort timescales (∼105 years), or even dikes that formed when theMoon resided at any of the SPA-produced paleopoles, cannot fully
explain the observed diversity of magnetic paleopoles we ﬁnd.
There are many unknowns in modeling the paleopoles allowed by SPA’s formation. For example, reorientation and magnetic anomaly formation depends on the SPA impact’s effect on dynamo operation (Arkani-Hamed and Olson, 2010a, 2010b), the inertia tensor of the Moon just before SPA formation and the lengthscale (and thereby cooling timescale) of the materials making up
the anomalies. Furthermore, the spin vector of the Moon will befreely precessing around its angular momentum vector during thefree precession damping timescale. If the dynamo after SPA impact
retained its alignment along the Moon’s angular momentum vector, this precession could broaden the range of paleopoles permitted (essentially accessing a range of paleopoles around the pale-

172

M. Nayak et al. / Icarus 286 (2017) 153–192

Fig. 22. Paleopoles from SPA anomalies compared with published paleopoles from other anomalies. All ﬁgures show the southern hemisphere (all points are southernlatitudes). (A) Our results from Fig. 21, with all north poles reversed into the southern hemisphere. The paleopole obtained from the mean of all data at area 2 (see Fig.
21) has been omitted in this ﬁgure. Red and blue triangles are possible lunar paleopoles resulting from SPA’s gravity anomaly immediately after its formation, assuming
an extreme range of density anomaly models; red (blue) points are positive (negative) density anomalies (see Section 4.3.4). SPA’s center (a limiting paleopole for very
negative density models of SPA) is shown at (−53.2° S, 191° E) (black triangle). The magenta triangle shows the paleopole dervied from the tidal component of the Moon’s
topography, outside of large basins (Garrick-Bethell et al., 2014). (B) Paleopoles from the 11 anomalies reported by Takahashi et al. (2014). No anomalies are within SPA. Blue(red) points represent inversions from Kaguya (Lunar Prospector) data. In many cases, multiple points represent paleopoles from single sites. Stars represent the means of
the two clusters. (C) Paleopoles from 10 isolated magnetic anomalies reported by Arkani-Hamed et al. (2014). (D) Paleopoles from crustal magnetization at 20 impact craters
reported by Arkani-Hamed et al. (2014). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)M. Nayak et al. / Icarus 286 (2017) 153–192

opoles shown in Fig. 22a, with the range depending on the precession angle). Therefore, we cannot deﬁnitively rule out polar wander
processes as the origin for some of the diverse paleopole locations
at SPA.
A variant of the hypothesis of iron-rich SPA ejecta (Wieczorek
et al., 2012) is the cooling of iron-rich “sesquinary” (Zahnle et al.,2008) impactors formed by the SPA impact. Ejected into orbit immediately after impact, studies on various planetary bodies showthat these impactors can return to a body from 104 years (Nayaket al., 2016) to 106 years (Gladman et al., 1995) post-impact, atapproximately escape velocity (Nayak and Asphaug, 2016). With
reimpacts spread across this timescale, reaccreting iron-rich material originally from SPA, either still hot from the impact or heating
upon re-impact, could record a diverse set of orientations as theMoon reorients in response to moment of inertia changes. Acrossthis timescale, impact locations become random; such sesquinary
magnetism caused by these iron-rich impactors would be widelydistributed around the Moon.
5. Conclusions
Using two different inversion methods, we ﬁnd diverse directions of magnetization among magnetic anomalies in the northern SPA basin. The diverse directions help rule out impact-relatedﬁelds as their only origin. Intrusive bodies, such as the dikes proposed by (Purucker et al., 2012), are a plausible explanation. The
diverse paleopole locations could imply large amounts of true polar wander, but true polar wander inferred independently fromgravity (Garrick-Bethell et al., 2014) and hydrogen deposits (Siegler
et al., 2016) implies more modest changes in the Moon’s orientation. The diverse directions argue against the hypothesis that they
all formed from iron-rich SPA ejecta that cooled in a dynamo ﬁeld
(Wieczorek et al., 2012). A simple gravity anomaly model for large
amounts of SPA-produced reorientation fails to explain at least ﬁveof the paleopoles observed, but many unknowns remain in modeling this process. Some SPA ejecta may have produced iron-rich

173

sesquinary impactors that landed across the Moon and recorded
orientation changes as they cooled in a dynamo ﬁeld. A dynamothat was not aligned with the lunar spin axis remains a plausible explanation for all observations, but gaps remain in our ability
conﬁrm this hypothesis. The wide variety of viable hypotheses andlarge number of unknowns highlight the complexity of interpretingthe origins of lunar magnetic anomalies and their paleopoles.
Acknowledgments
This research was supported by the NASA Discovery Data Analysis Program (grant NNX16AJ07G), and the BK21 plus program
through the National Research Foundation (NRF), funded by the
Ministry of Education of Korea. MN acknowledges support fromthe National Defense Science and Engineering Graduate (NDSEG)
Fellowship, 32 CFR 168a, funded by the Air Force Oﬃce of Scientiﬁc Research, Department of Defense. Partial support was also
provided by (1) The Alfred P. Sloan Foundation, (2) The UC Santa
Cruz / NASA Ames University Aﬃliated Research Center, Aligned
Research Program and (3) Red Sky Research, LLC. Thanks to Rob
Coe and Gary Glatzmaier for helpful comments. Opinions, interpretations and recommendations expressed are those of the authors
and are not necessarily endorsed by the US Air Force or the Department of Defense.
Appendix A. Comparison of magnetic ﬁeld observations with
best-ﬁt models for all datasets
We use 45 datasets for the analysis in this work, the total ﬁeld
plots of which are shown in Fig. 2 (areas 1–5) and Fig. 3 (areas
6–10). In Figs. 8–17 we presented example magnetic ﬁeld maps forthe best-ﬁt source magnetizations, from both inversion algorithms,
for areas 1–10, respectively. Here we present the magnetic ﬁeld
maps for the best-ﬁt source magnetizations for the 35 remaining
datasets.

M. Nayak et al. / Icarus 286 (2017) 153–192

−23
−24
−25
−26
−27

10

−23
−24
−25
−26
−27

10

0

0
−5

TOTAL (nT)55EAST (nT)

−23
−24
−25
−26
−27

10

0
−10

0

−186 −184 −182

−186 −184 −182

−186 −184 −182

−10

NORTH (nT)

Lat (deg N)
Lat (deg N)

GD−VM

Lat (deg N)

DD−CM

Lat (deg N)

KG Data
−23
−24
−25
−26
−27

RADIAL (nT)174

Fig. A1. Area 1 best-ﬁt results for both inversion algorithms. (Left column) K-MAG observations in the lunar wake on 2009 day 123, compared to (Middle column) a modelwith manually placed dipoles of equal magnetization (DD-CM) and (Right column) a grid-based model with variable magnetizations (GD-VM). Source dipoles in the DD-CM
algorithm are approximately located at maxima of the observed total ﬁeld (black dots). White dots in the east panel are locations of spacecraft magnetometer measurements.
The mean measurement altitude is 31.6 km. Longitude (degrees east) is indicated on the x-axis.

5
0

5
0
−5

5
0
−5

TOTAL (nT)

10

EAST (nT)

−23
−24
−25
−26
−27

GD−VM

NORTH (nT)

Lat (deg N)

−23
−24
−25
−26
−27

DD−CM

10
0
−10
−186 −184 −182

−186 −184 −182

−186 −184 −182

RADIAL (nT)

Lat (deg N)
Lat (deg N)

−23
−24
−25
−26
−27

Lat (deg N)

LP Data
−23
−24
−25
−26
−27

Fig. A2. Area 1 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 31.1 km and the observations are from1999 day 61.

M. Nayak et al. / Icarus 286 (2017) 153–192

100

TOTAL (nT)

20

10
0
−10

10
0
−10

10
0
−10
−186 −184 −182

−186 −184 −182

−186 −184 −182

EAST (nT)

−23
−24
−25
−26
−27

GD−VM

NORTH (nT)

Lat (deg N)

−23
−24
−25
−26
−27

DD−CM

RADIAL (nT)

Lat (deg N)
Lat (deg N)

−23
−24
−25
−26
−27

Lat (deg N)

LP Data
−23
−24
−25
−26
−27

175

Fig. A3. Area 1 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 21.4 km and the observations are from1999 day 142.

DD−CM

5

−22
−24

0

−18

5

−20

0

−22
−24

−5

−18

5

−20

0

−22
−24

−5

−18

5

−20

0

−22
−24
−194

−5
−190

−186

−194

−190

−186

−194

−190

−186

EAST (nT)−20

NORTH (nT)

10

−18

TOTAL (nT)GD−VM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

KG Data

Fig. A4. Area 2 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 43.5 km and the observations are from2009 day 69.

M. Nayak et al. / Icarus 286 (2017) 153–192

Lat (deg N)

DD−CM

−20

5

−22
−24

0
5

−18
−20

0

−22

−5

−24

5

−18
−20

0

−22

−5

−24

EAST (nT)10−18

TOTAL (nT)GD−VM

NORTH (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

KG Data

10
−18
−20

0

−22
−24
−194

−10
−190

−186

−194

−190

−186

−194

−190

−186

RADIAL (nT)176

Fig. A5. Area 2 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 38.5 km and the observations are from2009 day 96.

DD−CM

10

−20
−22
−24

0
10

−18
−20

0

−22

−10

−18

5
0
−5

−20
−22
−24

−18

10

−20

0

−22
−24
−194

NORTH (nT)

−24

−10
−190

−186

−194

−190

−186

−194

−190

−186

EAST (nT)−18

TOTAL (nT)GD−VM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

Fig. A6. Area 2 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 31.1 km and the observations are from1999 day 33.

M. Nayak et al. / Icarus 286 (2017) 153–192

−25

GD−VM
8
6
4
2
0

−26
−27
−28
−29

TOTAL (nT)

Lat (deg N)

−25

DD−CM

EAST (nT)

Lat (deg N)
Lat (deg N)−25

Lat (deg N)

KG Data
−25

177

2
0
−2

−26
−27
−28

NORTH (nT)

−29

5

−26

0

−27
−28

RADIAL (nT)

−5

−29

5
0
−5

−26
−27
−28
−29
−195

−190

−195

−190

−195

−190

Fig. A7. Area 3 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 42.8 km and the observations are from2009 day 69.

−27
−28
−29

2

−26
−27

0

−28

−2

−29

5

−26

0

−27
−28

−5

−29
−25

5

−26
−27

0

−28

−5

−29
−195

−190

−195

−190

−195

−190

EAST (nT)

8
6
4
2

−26

TOTAL (nT)GD−VM

NORTH (nT)

−25

DD−CM

RADIAL (nT)

Lat (deg N)
Lat (deg N)−25

Lat (deg N)−25

Lat (deg N)

KG Data

Fig. A8. Area 3 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 37.6 km and the observations are from2009 day 96.

M. Nayak et al. / Icarus 286 (2017) 153–192

5

−25

5

−28

−26
−27

0

−28

−5

−29

TOTAL (nT)

−25

−27

EAST (nT)

0

10

−26

10

−26

0

−27
−28
−29

−10

−25

10

−26

0

−27
−28

−10

−29
−195

−190

−195

−190

−195

−190

NORTH (nT)−29

Lat (deg N)

Lat (deg N)

GD−VM

Lat (deg N)

DD−CM

Lat (deg N)

LP Data−25

RADIAL (nT)178

Fig. A9. Area 3 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 29.6 km and the observations are from1999 day 61.

40

−32
−33

20

−34
−35

−32
−34
−35

−20

NORTH (nT)

10
0
−10

−33

40
20
0
−20
−40

20

−32
−33

0

−34
−35

−32
−33
−34
−35
−186

TOTAL (nT)

GD−VM

−184

−182

−186

−184

−182

−186

−184

−182

EAST (nT)DD−CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

KG Data

Fig. A10. Area 4 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 12.1 km and the observations are from2009 day 151.

M. Nayak et al. / Icarus 286 (2017) 153–192

−33

5

−34
−35

0
4
2
0
−2
−4

−32
−33
−34
−35

5

−32
−33

0

−34
−35

−5
10

−32
−33

0

−34
−35
−186

EAST (nT)10

TOTAL (nT)GD−VM

NORTH (nT)

DD−CM

−32

−184

−182

−186

−184

−182

−186

−184

−182

−10

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data179

Fig. A11. Area 4 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 30.9 km and the observations are from1999 day 115.

10

−34
−35

0

−32

5

−33

0

−34
−35

−5

−32

10

−33

0

−34
−35

−10

−32

10

−33

0

−34

−10

−35
−186

EAST (nT)−33

TOTAL (nT)

GD−VM

−184

−182

−186

−184

−182

−186

−184

−182

NORTH (nT)DD−CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data−32

Fig. A12. Area 4 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 21.8 km and the observations are from1999 day 142.

M. Nayak et al. / Icarus 286 (2017) 153–192

GD−VM

0

Lat (deg N)

−16

−22

4
2
0
−2
−4

EAST (nT)−16

5

NORTH (nT)

−16

−18

5

−20

−18
−20

−18

0

−20

−5

−22

5
−18

0

−20
−22

−5
−165

−160

−165

−160

−165

−160

RADIAL (nT)

Lat (deg N)−22

Lat (deg N)

DD−CM

−16

Lat (deg N)

KG Data

TOTAL (nT)180

Fig. A13. Area 5 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 39.1 km and the observations are from2009 day 96.

20

−20

10

−22

0

TOTAL (nT)

−18

10

−18

0

−20

−10

−22

10

−18

0

−20

−10

−22

20

−18

0

−20
−22

−20
−165

−160

−165

−160

−165

−160

EAST (nT)

−16

GD−VM

NORTH (nT)

Lat (deg N)

−16

DD−CM

RADIAL (nT)

Lat (deg N)
Lat (deg N)−16

Lat (deg N)

KG Data−16

Fig. A14. Area 5 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 13.2 km and the observations are from2009 day 151.

M. Nayak et al. / Icarus 286 (2017) 153–192

−16

−16

−22

−10

−16

10

5

−22

0
5

−18

0

−20
−22

−5
10

−18

0

−20

−18

0

−20
−22

−10
−165

−160

−165

−160

−165

−160

NORTH (nT)−20

RADIAL (nT)

10

−18

EAST (nT)

Lat (deg N)
Lat (deg N)

TOTAL (nT)GD−VM

Lat (deg N)

DD−CM

−16

Lat (deg N)

LP Data

181

GD−VM
6
4
2
0

−24

−22

TOTAL (nT)

DD−CM

−22
−24

−22
−24
−168 −166 −164 −162

−168 −166 −164 −162

−168 −166 −164 −162

4
2
0
−2
−4

NORTH (nT)

2
0
−2

−24

EAST (nT)KG Data
−22

5
0
−5

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Fig. A15. Area 5 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 22.3 km and the observations are from1999 day 142.

Fig. A16. Area 6 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 39.1 km and the observations are from2009 day 96.

5

24

0

4
2
0

22
24

2
4

5

22

0

24

5

22

5
0

24

5
168

166

164

162

168

166

164

162

168

166

164

162

EAST (nT)GD VM

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)
Lat (deg N)
Lat (deg N)

LP Data22

TOTAL (nT)M. Nayak et al. / Icarus 286 (2017) 153–192

Lat (deg N)

182

Fig. A17. Area 6 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 32.8 km and the observations are from1999 day 33.

Fig. A18. Area 6 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 32.5 km and the observations are from1999 day 115.

M. Nayak et al. / Icarus 286 (2017) 153–192

183

28

6
4
2

26

2

27

0

27

28

TOTAL (nT)

GD VM

2

4
2
0

26
27
28

26

2
4

5
0

27
28

5
168

166

164

162

168

166

164

162

168

166

164

162

EAST (nT)DD CM

NORTH (nT)KG Data
26

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Fig. A19. Area 6 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 22.6 km and the observations are from1999 day 142.

Fig. A20. Area 7 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 39.0 km and the observations are from2009 day 96.

8
6
4
2

27
28

26

2
0

27

2

28

26

5
0

27

5

28

26

5
0

27

5

28
168

166

164

162

168

166

164

162

168

166

164

162

EAST (nT)GD VM

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)
Lat (deg N)
Lat (deg N)

LP Data26

TOTAL (nT)M. Nayak et al. / Icarus 286 (2017) 153–192

Lat (deg N)

184

5

27
28

0

26

TOTAL (nT)

GD VM

2
0

27

2

28

26

5
0

27

5

28

26

5
0

27

5

28
168

166

164

162

168

166

164

162

168

166

164

162

EAST (nT)DD CM

NORTH (nT)LP Data
26

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Fig. A21. Area 7 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 32.6 km and the observations are from1999 day 33.

Fig. A22. Area 7 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 32.4 km and the observations are from1999 day 115.

32

4
2
0

30

2
0

32

TOTAL (nT)

GD VM

2

5

30

0

32

5

2
0

30
32
170

2
165

160

170

165

160

170

165

160

EAST (nT)

DD CM

30

NORTH (nT)KG Data

185

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

M. Nayak et al. / Icarus 286 (2017) 153–192

Lat (deg N)

Lat (deg N)

Lat (deg N)

KG Data

DD CM

GD VM

30

20

32

10
0

TOTAL (nT)

Lat (deg N)

Fig. A23. Area 8 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 39.4 km and the observations are from2009 day 96.

10

30

0

32

10

10
0
10

30
32

20

30

0

32

20
170

165

160

170

165

160

170

165

160

Fig. A24. Area 8 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 13.0 km and the observations are from2009 day 151.

30

5

32

0

5
30

0

32

5

5
0

30
32

5

5

30

0

32

5
170

165

160

170

165

160

170

165

160

EAST (nT)GD VM

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)
Lat (deg N)
Lat (deg N)

LP Data

TOTAL (nT)M. Nayak et al. / Icarus 286 (2017) 153–192

Lat (deg N)

186

DD CM

GD VM

30

10

32

5
0

TOTAL (nT)LP Data

10

30

0

32

5
0

30
32

5

5
0

30
32
170

5
165

160

170

165

160

170

165

160

NORTH (nT)10RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Fig. A25. Area 8 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 32.7 km and the observations are from1999 day 61.

Fig. A26. Area 8 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 23.6 km and the observations are from1999 day 142.

M. Nayak et al. / Icarus 286 (2017) 153–192

10

36

5
38

0
34

5

36

0

38

5
10

34

0

36
38

10

34

10

36

0

38
200

EAST (nT)34

TOTAL (nT)

GD VM

10
195

200

195

200

195

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

KG Data187

Fig. A27. Area 9 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 37.0 km and the observations are from2009 day 96.

36
38

10
34

0

36
38

10
34

10

36

0

38

10

34

10

36

0
10

38
200

195

200

195

200

195

EAST (nT)

20
15
10
5

34

TOTAL (nT)GD VM

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

Fig. A28. Area 9 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 28.7 km and the observations are from1999 day 61.

M. Nayak et al. / Icarus 286 (2017) 153–192

Lat (deg N)

34

20

36

10

38

0
34

5
0

36

5

38

TOTAL (nT)GD VM

EAST (nT)

DD CM

20

34

0

36
38

20
34

10
0
10

36
38
200

195

200

195

200

195

NORTH (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

RADIAL (nT)188

Fig. A29. Area 9 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 28.5 km and the observations are from1999 day 88.

36

10

38

0
34

10

36

0

38

10

34

10

36

0

38

10

34

10
0
10

36
38
200

195

200

195

200

195

EAST (nT)20

34

TOTAL (nT)GD VM

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

Fig. A30. Area 9 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 28.7 km and the observations are from1999 day 115.

M. Nayak et al. / Icarus 286 (2017) 153–192

30
20
10
0

36
38

34

10

36

0
10

38

34

20

36

0

38

20

34

20

36

0
20

38
200

195

200

195

200

195

EAST (nT)34

TOTAL (nT)GD VM

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

189

Fig. A31. Area 9 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 19.6 km and the observations are from1999 day 142.

16

0
12

20

14

0

16

20
50

12
14

0

16

50
12

50

14

0

16

50
156 154 152

156 154 152

156 154 152

EAST (nT)50
14

TOTAL (nT)GD VM

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

KG Data
12

Fig. A32. Area 10 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 13.9 km and the observations arefrom 2009 day 151.

M. Nayak et al. / Icarus 286 (2017) 153–192

Lat (deg N)

12
10
8
6
4
2

12
14
16

4
2
0

12
14

2
4

16

TOTAL (nT)GD VM

EAST (nT)

DD CM

10

12
14

0

16

10
12

10

14

0

16

10
156 154 152

156 154 152

156 154 152

NORTH (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

RADIAL (nT)190

Fig. A33. Area 10 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 32.3 km and the observations arefrom 1999 day 61.

10

14

5
16

0
12

5

14

0

16

5

EAST (nT)12

TOTAL (nT)

GD VM

10

12
14

0

16

10
12

10

14

0

16

10
156 154 152

156 154 152

156 154 152

NORTH (nT)DD CM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

Fig. A34. Area 10 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 32.0 km and the observations arefrom 1999 day 88.

M. Nayak et al. / Icarus 286 (2017) 153–192

DD CM

10
16

0
12

10

14

0

16

10
20

12
14

0

16

20
12

20

14

0

16

20
156 154 152

156 154 152

156 154 152

EAST (nT)20

14

NORTH (nT)12

TOTAL (nT)GD VM

RADIAL (nT)

Lat (deg N)

Lat (deg N)

Lat (deg N)

Lat (deg N)

LP Data

191

Fig. A35. Area 10 best-ﬁt results for both inversion algorithms. Details are as in Fig. A1, except that the mean measurement altitude is 21.9 km and the observations arefrom 1999 day 142.

Supplementary materials
Supplementary Appendices B and C can be found, in the online
version, at doi:10.1016/j.icarus.2016.09.038.
References
Andrews-Hanna, J.C., Besserer, J., Head III, J.W., Howett, C.J.A., Kiefer, W.S., Lucey, P.J.,
McGovern, P.J., Melosh, H.J., Neumann, G.A., Phillips, R.J., 2014. Structure andevolution of the lunar Procellarum region as revealed by GRAIL gravity data.
Nature 514, 68–71.
Arkani-Hamed, J., Boutin, D., 2014. Analysis of isolated magnetic anomalies andmagnetic signatures of impact craters: evidence for a core dynamo in the early
history of the Moon. Icarus 237, 262–277. doi:10.1016/j.icarus.2014.04.046.
Arkani-Hamed, J., Olson, P., 2010a. Giant impact stratiﬁcation of the Martian core.
Geophys. Res. Lett. 37. doi:10.1029/2009GL041417.
Arkani-Hamed, J., Olson, P., 2010b. Giant impacts, core stratiﬁcation, and failure of
the Martian dynamo. J. Geophys. Res. 115, E07012. doi:10.1029/2010JE003579.
Arnold, J.R., 1975. Monte Carlo simulation of turnover processes in the lunar regolith. In: Proceedings of the Lunar Planetary Science Conference 6th,
pp. 2375–2395.
Blewett, D.T., Coman, E.I., Hawke, B.R., Gillis-Davis, J.J., Purucker, M.E., Hughes, C.G.,
2011. Lunar swirls: examining crustal magnetic anomalies and space weathering
trends. J. Geophys. Res. 116, E02002. doi:10.1029/2010JE003656.
Butler, R., 1992. Paleomagnetism: Magnetic Domains to Geologic Terranes. Blackwell
Scientiﬁc Publications, Boston, p. 319.
Coleman, P.J., Schubert, G., Russell, C.T., Sharp, L.R., 1972. Satellite measurements of
the Moon’s magnetic ﬁeld: a preliminary report. Moon 4, 419–429. doi:10.1007/
BF0 05620 08.
Dwyer, C.A., Stevenson, D.J., Nimmo, F., 2011. A long-lived lunar dynamo driven bycontinuous mechanical stirring. Nature 479, 212–214. doi:10.1038/nature10564.
Dyal, P., Parkin, C.W., Daily, W.D., 1974. Magnetism and the interior of the Moon.
Rev. Geophys. Space Phys. 12. doi:10.1029/RG012i0 04p0 0568.
Fisher, R.A. Dispersion on a sphere. Proceedings of the Royal Society of London Series A. 217, 295–305, http://dx.doi.org/10.1098/rspa.1953.0064
Garrick-Bethell, I, 2016. A simple history of lunar true polar wander. 47th Lunar
Planetary Science Conference, abstract 2874.
Garrick-Bethell, I., Perera, V., Nimmo, F., Zuber, M.T., 2014. The tidal–rotational shapeof the Moon and evidence for polar wander. Nature 512, 181–184. doi:10.1038/
nature13639.
Garrick-Bethell, I., Zuber, M.T., 2009. Elliptical structure of the lunar South PoleAitken basin. Icarus 204, 399–408. doi:10.1016/j.icarus.2009.05.032.

Gladman, B.J., Burns, J.A., Duncan, M.J., Levison, H.F. The dynamical evolution of lunar impact ejecta. Icarus. 118, 302–321. http://dx.doi.org/10.1006/icar.1995.1193.
Goldreich, P., Toomre, A., 1969. Some remarks on polar wandering. J. Geophys. Res.
74, 2555–2567. doi:10.1029/JB074i010p02555.
Gong, S., Wieczorek, M.A., 2016. Is the lunar magnetic ﬁeld correlated with gravity
or topography? In: 47th Lunar and Planetary Science Conference, abstract 2290.
Halekas, J.S., Brain, D.a., Lin, R.P., Mitchell, D.L., 2008. Solar wind interaction withlunar crustal magnetic anomalies. Adv. Space Res. 41, 1319–1324. doi:10.1016/j.
asr.20 07.04.0 03.
Halekas, J.S., Poppe, A., Delory, G.T., Farrell, W.M., Hor, M., 2012. Solar wind electron interaction with the dayside lunar surface and crustal magnetic ﬁelds: Evidence for precursor effects. Earth Planets Space 64, 73–82. doi:10.5047/eps.2011.
03.008.
Hemingway, D., Garrick-Bethell, I., 2012. Magnetic ﬁeld direction and lunar swirl
morphology: insights from Airy and Reiner Gamma. J. Geophys. Res. 117, E10012.
doi:10.1029/2012JE004165.
Hood, L., Williams, C., 1989. The lunar swirls: distribution and possible origins. In:
Proceedings of the 19th Lunar Planetary Science Conference, pp. 99–113.
Hood, L., et al., 2011. Initial mapping and interpretation of lunar crustal magnetic
anomalies using lunar prospector magnetometer data. J. Geophys. Res. 106,27825–27839.
Hood, L.L., Artemieva, N.A., 2008. Antipodal effects of lunar basin-forming impacts:
initial 3D simulations and comparisons with observations. Icarus 193, 485–502.
doi:10.1016/j.icarus.2007.08.023.
Hood, L.L., Huang, Z., 1991. Formation of magnetic anomalies antipodal to lunar impact basins: two dimensional model calculations. J. Geophys. Res. 96,
9837–9846.
Irving, E., 1964. Paleomagnetism and Its Application to Geological and Geophysical
Problems, 1st ed.. John Wiley & Sons, Inc, New York, p. 399.
Kramer, G.Y., et al., 2011. Characterization of lunar swirls at Mare Ingenii: a model
for space weathering at magnetic anomalies. J. Geophys. Res. 116, E04008.
doi:10.1029/2010JE003669.
Kurata, M., Tsunakawa, H., Saito, Y., Shibuya, H., Matsushima, M., Shimizu, H., 2005.
Mini-magnetosphere over the Reiner Gamma magnetic anomaly region on theMoon. Geophys. Res. Lett. 32, L24205. doi:10.1029/2005GL024097.
Le Bars, M., Wieczorek, M.A., Karatekin, Ö., Cébron, D., Laneuville, M., 2011. An
impact-driven dynamo for the early Moon. Nature 479, 215–218. doi:10.1038/
nature10565.
Lin, R.P., et al., 1998. Lunar surface magnetic ﬁelds and their interaction with thesolar wind: results from lunar prospector. Science 281, 1480–1484. doi:10.1126/
science.281.5382.1480.
Mitchell, D.L., et al., 2008. Global mapping of lunar crustal magnetic ﬁelds by Lunar
Prospector. Icarus 194, 401–409. doi:10.1016/j.icarus.2007.10.027.

192

M. Nayak et al. / Icarus 286 (2017) 153–192

Nayak, M., Nimmo, F., Udrea, B., 2016. Effects of mass transfer between Martiansatellites on surface geology. Icarus 267, 220–231. doi:10.1016/j.icarus.2015.12.
026.
Nayak, M.V, Asphaug, E., 2016. Sesquinary Catenae on the Martian satellite phobos from reaccretion of escaping ejecta. Nat. Commun. 7, 12591. doi:10.1038/
NCOMMS12591.
Nicholas, J.B., Purucker, M.E., Sabaka, T.J., 2007. Age spot or youthful marking: Origin
of Reiner Gamma. Geophys. Res. Lett. 34, L02205. doi:10.1029/2006GL027794.
Parker, R.L., 1988. A statistical theory of seamount magnetism. J. Geophys. Res. 93,
3105–3115. doi:10.1029/JB093iB04p03105.
Parker, R.L., 1991. A theory of ideal bodies for seamount magnetism. J. Geophys. Res.
96, 16101–16112. doi:10.1029/91JB01497.
Peale, S.J., 1976. Excitation and relaxation of the wobble, precession, and librationof the Moon. J. Geophys. Res. 81, 1813–1827. doi:10.1029/JB081i011p01813.
Purucker, M.E., 2008. A global model of the internal magnetic ﬁeld of the Moonbased on lunar prospector magnetometer observations. Icarus 197, 19–23.
doi:10.1016/j.icarus.2008.03.016.
Purucker, M.E., Head, J.W., Wilson, L., 2012. Magnetic signature of the lunar South
Pole-Aitken basin: character, origin, and age. J. Geophys. Res. 117, E05001.
doi:10.1029/2011JE003922.
Purucker, M.E., Nicholas, J.B., 2010. Global spherical harmonic models of the internal
magnetic ﬁeld of the Moon based on sequential and coestimation approaches.
J. Geophys. Res. 115, E12007. doi:10.1029/2010JE003650.
Runcorn, S.K., 1983. Lunar magnetism, polar displacements and primeval satellites
in the Earth–Moon system. Nature 304, 589–596. doi:10.1038/304589a0.

Siegler, M.A., et al., 2016. Lunar true polar wander inferred from polar hydrogen.
Nature 531, 480–484. doi:10.1038/nature17166.
Smith, D.E., et al., 2010. Initial observations from the lunar orbiter laser altimeter(LOLA). Geophys. Res. Lett. 37, 1–6. doi:10.1029/2010GL043751.
Takahashi, F., Tsunakawa, H., Shimizu, H., Shibuya, H., Matsushima, M., 2014. Reorientation of the early lunar pole. Nat. Geosci. 7, 409–412. doi:10.1038/ngeo2150.
Tsunakawa, H., Takahashi, F., Shimizu, H., Shibuya, H., Matsushima, M., 2014. Regional mapping of the lunar magnetic anomalies at the surface: method and its
application to strong and weak magnetic anomaly regions. Icarus 228, 35–53.
doi:10.1016/j.icarus.2013.09.026.
Weiss, B.P., Tikoo, S.M., 2014. The lunar dynamo. Science 346, 1198. doi:10.1126/
science.1246753.
Wieczorek, M.A., Weiss, B.P., 2010. Testing the lunar dynamo hypothesis using global
magnetic ﬁeld data. 41st Lunar and Planetary Science Conference, abstract 1625.
Wieczorek, M.A., Weiss, B.P., Stewart, S.T., 2012. An impactor origin for lunar magnetic anomalies. Science 335, 1212–1215. doi:10.1126/science.1214773.
Williams, J., Boggs, D., Yoder, C.F., Ratcliff, J.T., Dickey, J.O., 2001. Lunar rotational
dissipation in solid body and molten core. J. Geophys. Res. 106, 27933–27968.
doi:10.1029/20 0 0JE0 01396.
Zahnle, K., Alvarellos, J.L., Dobrovolskis, A., Hamill, P., 2008. Secondary andsesquinary craters on Europa. Icarus 194, 660–674. doi:10.1016/j.icarus.2007.10.
024.
Zuber, M.T., et al., 2013. Gravity ﬁeld of the Moon from the gravity recovery andinterior laboratory (GRAIL) mission. Science 339, 668–671. doi:10.1126/science.
1231507.