Systematic Uncertainties in the Measurement of Neutron lifetime
Using Lunar Prospector Neutron Spectrometer

LP consists of two geometrically identical neutron detector proportional counter tubes, each of 20 cm in length and 5.7 cm in diameter (Feldman et al., 2004). The tubes are filled with ³He gas at 10 atm pressure and covered with a 0.63 mm thick sheet of either Cd or Sn. The Cd absorbs low energy neutrons, so the “HeCd” detector is sensitive to neutrons of energy $>0.4$ eV (epithermal neutrons) while the “HeSn” detector is sensitive to neutrons of all energies (thermal and epithermal neutrons). The data product from these detectors is a 32-bin pulse amplitude histogram that is produced at a 32 second cadence. We sum the bins in the second half of the pulse height spectra, bins corresponding to higher energy deposition events in the detector to reduce the contribution of non-neutron background interactions, corresponding to energy bins 17 to 32. This produces the neutron counts per 32 seconds data product for the HeCd and HeSn detectors and the corresponding uncertainties in measurements. This choice was made to avoid the background noise and to achieve higher signal-to-noise in the pulse height spectrum (Maurice et al., 2004). To get the thermal neutron counts (energy range of $<$0.4 eV), we subtract the HeCd count rates from the HeSn count rates.

We construct simulated neutron counts in three steps: first by simulating the neutron emission from a given surface composition and temperature, then by mapping these compositions onto the lunar surface and calculating their contribution to the neutron current at a given point in the LP orbit, and finally by evaluating the response of the two NS detectors to the current to produce a count rate that is directly comparable to the detector count rates.

Wilson et al. (2021) derived five lunar compositions from Peplowski et al. (2016) and Jolliff et al. (2000), which are re-binned into a 2^∘ resolution map for this work. These compostions are the procellarum KREEP terrane (PKT), South Pole Aitken terrane (SPA), feldspathic highland terrane (FHT), pure anorthosite terrane (PAN) and non-PKT maria (nPKT). For elemental mass fraction of these composition, refer to Table II in Wilson et al. (2021). Although this map has been used to estimate the value of neutron lifetime, it is a simplified map, and fails to capture the complexities of the lunar terrain. Here, in addition to the five listed compositions, we also model the neutron current from the lunar surface at two more resolutions derived from Prettyman et al. (2006) - (i) a 5^∘ resolution map showing 1790 different compositions and (ii) a 20^∘ resolution map showing 114 different compositions. The 5^∘ binned composition map comes from the NASA Planetary Data System (PDS) (LPP, 2002) and contains weight percent (wt%) or ppm of Mg, O, Al, Si, Ca, Ti, Fe, K, Th, and U. Additional higher resolution (2^∘) binned maps for H and Sm are included from the PDS. To form a complete set, the H and Sm were resampled to 5^∘ and 20^∘ resolutions. Further, Gd, an important neutron absorbing element, was included by scaling the amount of Sm by a constant factor of 1.3, adopted from Hidaka et al. (2000).

To study the effects of latitude-dependent surface temperature, we performed simulations for two different cases using the surface composition maps at 20^∘ resolution reported in Prettyman et al. (2006): (i) constant equatorial temperature of 240 K and (ii) latitude-dependent subsurface lunar temperature reported by Little et al. (2003). The temperature (in $K$) is given by:

The temperature is modified using the TMP card in MCNP, which accounts for temperature effects on the elastic interaction cross sections assuming a free gas model. It also accounts for the thermal motion and modified collision kinematics reported in Werner et al. (2018).

MCNP^® is Los Alamos National Laboratory’s general-purpose Monte Carlo N-Particle code¹¹1https://mcnp.lanl.gov/ that can be used for detector design analysis and particle transport analysis (Werner et al., 2018). We collect neutron current tallies on a spherical shell above the lunar surface. This, combined with the NS detector response gives the neutron count rates. For this work, MCNP version 6.2 was used.

To simulate the GCRs we used the particle flux equations in Usoskin et al. (2017) and the reported mean solar modulation during the orbital insertion of LP (Jan 1998) of 505 MV, which results in a GCR proton flux of 3.575 particles/cm²-s and $\alpha$ particle flux of 0.323 particles/cm²-s. We simulate only protons, but account for the contribution of $\alpha$ flux using the scaling from McKinney et al. (2006) of a factor of 3.8 more neutrons produced per incident $\alpha$ than per incident proton. Thus the total equivalent omni-directional proton flux is $3.575+3.8\times 0.323=4.8$ p/cm²-s. Since our MCNP simulations start with a $2\pi$ surface instead of total omni-directional flux, the incident flux is a factor of 4 lower (McKinney et al., 2006, Mesick et al., 2018), resulting in an incident simulated proton flux of 1.2 particles/cm²-s.

We tally the neutrons emitted after the interaction of GCRs with the lunar surface using MCNP current tallys and normalize for source flux (1.2 particles/cm²-s) and tally area. We also tested the angular distributions of a few compositions and found the neutron current to be independent of the cos($\theta$)^1.25 relationship, $\theta$ being angle of neutron travel relative to the tally surface. This is unlike the flux tally which needs additional appropriate solid angle normalization. To keep the simulations computationally minimal, we do not use the angle-dependent outputs. For each of the compositions tested, the normalized neutron current integrated over the thermal peak ($<$0.4 eV, when neutrons have slowed to an energy/velocity equal to a temperature dependent thermal equilibrium) is noted to be $\approx$0.2 particles/cm²-s.

Figure 1 shows the simulated neutron current for the five compositions reported in Wilson et al. (2021) at an altitude of 10 km above the lunar surface. The neutrons initially generated by the GCR interactions have energies on the order of 10s-100s of MeV (Mesick et al., 2018). The neutrons lose energy through elastic and inelastic scattering, eventually ending up in thermal equilibrium with the surface. The spectrum of neutrons emitted from the surface depends on the surface composition through the likely energy loss per scatter and the likelihood of absorption. The population of neutrons within the thermal, epithermal and fast regions of the neutron energy spectrum are counted to give insights into the surface composition for planetary science measurements. We observe a 20% variation in the epithermal and fast neutron current; however, these differences are not easily discernible in Figure 1 due to the large range on the y-axis.

MCNP is designed to compute neutron current by creating surfaces through which neutrons cross after they are released from a planetary surface, such as the Moon. This approach assumes that the neutron escapes the gravitational field of the planetary body. However, low-energy neutrons are gravitationally-bound and can return, interact with the surface, and be re-emitted, which increases the neutron current at low energies(Feldman et al., 1989). It may therefore be necessary to account for the low-energy neutrons with an energy less than 0.0295 eV, that do not escape lunar gravity. These effects can be accounted for using the FIELD card in MCNP6.2+ and MCNPX, which were implemented based on particle physics discussed in Feldman et al. (1989). However, we tested the effects of gravitationally bound downward travelling neutrons using the 5^∘ map and found that MCNP 6.2 uses the incorrect value of $\tau_{n}$ resulting in an over-prediction of the gravity effects. We estimate the gravitational effect to be less than 0.5% on the total thermal neutron flux for lunar compositions with the low water content (tens of ppm) used in this study. We therefore disabled the FIELD card in our simulations.

It is established that the neutron spectrum produced using the compositions as measured by GRS does not match with the neutron spectrometer data (Peplowski et al., 2013). To account for this, we scale the elemental mass fraction of oxygen such that the weighted sum of the elemental microscopic neutron absorption cross sections matches the measured macroscopic neutron absorption cross section. The data for macroscopic absorption cross-section from the science-phase of the LP mission is obtained from Feldman et al. (2000) and GRS measurements are obtained from the PDS (LPP, 2002). The microscopic neutron absorption cross sections for each of the abundance-weighted sum of the elemental constituents are obtained from Prettyman et al. (2006). We use the correction technique described in Equation 2 of Peplowski et al. (2013):

where $N_{A}$ is Avogadro’s number, $w_{i}$ is the elemental weight fraction, and $A_{i}$ is the atomic weight. This correction factor for some compositions was higher than others, but less than 30% overall. The corrected compositions, after scaling elemental mass fraction of oxygen were used as inputs to the MCNP simulations.

Figure 2 shows the integrated thermal neutron current for the 20^∘ and 5^∘ resolution maps before and after correction. An overall increase in integrated thermal neutrons is noted.

We start with neutrons from the spallation on the lunar surface, and use particle kinematics that take into account the ballistic trajectories due to gravity. The influence of lunar gravity is significant for low-energy neutrons used in the measurement of neutron lifetime. Neutrons with energies below 0.0295 eV that exit the lunar surface will become gravitationally bound. If these neutrons do not undergo beta decay during flight, they will return to the surface. On the contrary, higher-energy neutrons will follow hyperbolic trajectories, trading kinetic energy for potential energy as a function of radial distance. The orbits of these neutrons originate at all points on the lunar surface and intersect with a given satellite location. The closed-form solutions for neutron flux at a given satellite position due to neutron propagation in a spherically-symmetric gravitational field were presented by Feldman et al. (1989) (Feldman formalism). The solutions use the incident energy, angle, flux density, and transit time as a function of emission energy, angle, and radius. These are briefly described in Appendix VIII.1, including the correction of a sign error as noted in Mesick et al. (2020). This approach, commonly used in planetary science neutron analyses, does not account for the location of neutron emission. It assumes that the neutrons reaching a given satellite location originate from a single surface composition. This assumption is valid for planetary science analyses focused on probing surface composition and reflects the spatial resolution of those studies.

However, in this analysis, the satellite’s altitude is high enough to sample multiple composition units on the surface, thus supplementing the Feldman formalism is essential, such that the neutron flux can be mapped to its emission location. Specifically, we apply the Feldman formalism to generate a grid of neutron energies, angles, and transit times originating from the surface, along with their corresponding values at a given satellite location, for both upward-traveling neutrons and gravitationally-bound neutrons returning to the surface. This provides an analytical treatment of the changing neutron flux density from the surface-to-point problem. We then describe the neutron trajectories on that grid as orbital motion (following e.g. Curtis (2005)), calculating the orbit eccentricity and the true anomaly at both the satellite and the surface. The true anomaly for neutrons with inward trajectories intersecting the spacecraft is described in Appendix VIII.2. Since the true anomaly difference is simply an angular difference, this allows us to trace the neutron orbit from the satellite location back to its intersection with the surface.

Since the satellite’s velocity is comparable to that of the neutrons, it must be subtracted from the neutron velocity to obtain the neutron velocity in the rest frame of the detectors. We note that the spacecraft velocity provided by the mission ephemeris is in a global cartesian coordinate system with the z-axis aligned with the rotation axis of the Moon. To account for this, the neutron velocities are transformed from the coordinate system of the Feldman formalism, in which all angles are relative to local zenith, to the global coordinate system. This is described in VIII.3. After subtracting the satellite velocity, we perform a final coordinate transformation to the spacecraft’s reference frame. The spacecraft spins at approximately 12 rotations per minute, so only the spin axis is relevant. During most of the main science phase of the mission, the spin axis was roughly aligned (or anti-aligned) with the Moon’s rotation axis. However, this alignment was not in place during the orbital insertion period used in this analysis. We use the spacecraft attitude information from the ephemeris, which specifies the direction of the spacecraft spin axis as a right ascension (196-197^∘) and declination (58-59^∘) in the ECLIPJ2000 reference frame. We retrieve the coordinate axes of the Moon in the same reference frame using the SpiceyPy (Annex et al., 2020) wrapper of the SPICE toolkit (Acton Jr, 1996, Acton et al., 2018), and transform the neutron velocity into the reference frame of the LP spacecraft.

The detector response was simulated using Geant4 version 10.07 (Allison et al., 2016, 2006, Agostinelli et al., 2003), with a detector geometry derived from the original LPNS instrument description (Feldman et al., 1999). More recent measurements (Peplowski et al., 2020) have shown that similar cylindrical He³-filled gas proportional counters have a position-dependent response, shortening the region that detects incident neutrons. With this correction, the best current estimate of the active region of the LPNS detectors is 16.72 cm long rather than the original 20 cm specification (David Lawrence, email message to co-author, July 6, 2022). A 0.54 MeV energy deposited threshold was used in the simulation to correspond to energy bins 17 to 32 in the measurement data. This choice ensures that the pulse height spectra in the 32-bin pulse amplitude histogram is captures and the contribution of non-neutron background interactions is minimal. Further, only the proportional counter tubes and their immediate cladding were included in the simulation geometry.

Lunar Prospector was spin-stabilized, rotating around its axis with a period of roughly 5 seconds, which is much shorter than the 32 second integration time of the primary data set. The detector response was thus computed spin-averaged around the spacecraft rotation axis. The spacecraft bus itself shadows the detectors in some of the geometries. This effect was corrected empirically in the LPNS primary mission data (Maurice et al., 2004) rather than being modeled. However, those results are not applicable to the orbital insertion data used in this study. Thus, this effect is not corrected in our analysis. The results of Maurice et al. (2004) suggest that the impact of this effect is $\approx$ 30% shadowing for HeCd at the equator and reduces steadily as the spacecraft moves towards the poles. HeSn is expected to experience similar effects with an 180^∘ offset. The final spin-averaged detector responses are shown in Figure 3 for each detector, expressed as the detector effective area, i.e., the size of an equivalent perfectly efficient detector. Neutron flux was convolved with the detector response to obtain the detector count rates. The simulated detector count rate, for each of the two detectors, is thus the product of the incident neutron flux and effective area integrated over all energies and angles.

We subtract the HeCd count rates from the HeSn count rates to obtain the thermal neutrons (E$<$ 0.4 MeV), which are used as the model parameter to compute the neutron lifetime. The thermal neutron residual is computed through a simple subtraction of thermal neutrons using the measured and modeled detector count rates. Figure 4 shows simulated and measured HeCd and HeSn detector count rates from LP data using surface composition of 20^∘ resolution map, and the corresponding thermal neutron residuals between data and model. The computed detector count rates are assumed to be Poisson distributed, and the uncertainties ($\sigma$) are computed using 16 channel count rates containing the pulse height spectra of the 32 channel LP detector counts (Maurice et al., 2004). We note that the residuals have low variance ($\sim$ 1.7 $\sigma$ for the 20^∘ map), however retains structures when the spacecraft is in an active orbit maneuver.

To quantify the effects of the lunar sub-surface composition, the neutron current was modeled in MCNP using three different resolutions of composition maps, as described in Section II.2. Figure 5 shows the residuals in thermal neutrons between LP data and each of the three composition maps. All three resolutions of the maps have low statistical RMS in thermal neutron residuals. The composition map adopted from Wilson et al. (2021) that uses five compositions re-binned at 2^∘ resolution, produces the smallest RMS. These variances are reported in Table 1. We note that, although a simpler composition map has the least RMS in residuals, both 5^∘ and 20^∘ maps capture the variations in lunar terrain better than the five compositions. We also note that the composition maps are measurements made by LP during its mission period in a circular orbit, which could imply biases in this analysis. During the orbit insertion period considered in this analysis, the closest spacecraft approach is $\sim$85 km, while the farthest approach is $\sim$16,900 km. Using 1.5 times the altitude as an estimate for the spatial footprint, we calculate the spatial footprint to be 127.5 km, or 4.2^∘, at the closest approach. At this distance, the half-angle to the limb of the Moon is 17.6^∘, suggesting a significant reduction in contributions due to the non-uniform neutron emission angles. However, at the farthest distance, the half-angle increases to 84.7^∘, indicating that neutron emission is sourced from nearly half of the Moon’s surface. A more detailed and independent understanding of the lunar terrain and neutron emission is necessary to fully understand the effects of composition maps on the neutron lifetime.

To quantify the effects of the sub-surface temperature on neutron lifetime, two different cases are used in MCNP modeling - a (i) constant equatorial temperature at all latitudes and (ii) a latitude dependent temperature where poles are cooler than the equator. Figure 6 shows the neutron spectra for the two cases, using one of the compositions picked from the 20^∘ resolution composition map. For a constant equatorial temperature case, the poles are assumed to be hotter than they actually are. This produces a lower particle current and the thermal resonance is reached at a slightly higher energy. When the temperature is modeled as a function of the latitude, the thermal resonance is achieved at a relatively lower energy but with a relatively a higher particle current. Overall, higher number of thermal neutrons are seen when the effects of latitude dependent temperature is modeled.

We note that in the low energy regimes ($<0.2~{}eV$), the neutron current is relatively lower at the poles for the temperature varying case compared to the constant equatorial temperature due to lower subsurface temperatures. To quantify the effects of subsurface temperatures, we calculate the correlation coefficients between these two cases, given by:

where $f_{C}$ and $f_{V}$ represent normalized neutron current for constant equatorial temperature and latitude dependent subsurface temperature MCNP outputs respectively, and $\sigma_{C}$ and $\sigma_{V}$ represent corresponding standard deviations.

Figure 7 shows the difference in integrated thermal neutrons, correlation coefficients for neutrons of energy $<0.4~{}eV$, and the corresponding residuals in thermal neutrons for each of the two temperature conditions. For the calculation of correlation coefficients, a nominal value of $\tau_{n}$ = 880 s is assumed. Table 2 gives the RMS in the residuals shown in Figure 7. We note that, although a simple case of constant lunar equatorial temperature results in lower RMS, a latitude dependent temperature condition is realistic and such models are therefore necessary to capture the effects of sub-surface lunar temperatures in the measurement of neutron lifetime.

We take neutrons emitted due to spallation on the lunar surface and add complexities as it travels upwards towards the spacecraft and interacts with the HeCd and HeSn detectors. The piece of the analysis is referred to as propagation in this paper. The purpose of propagation is to compute the detector count rates using the neutron current from MCNP, which constitutes as the model. The first step in this analysis is to account for the ballistic trajectories of the neutrons due to lunar gravity. For a neutron of energy $K$ and mass $m$, leaving the surface, the binding energy (V) is given by

where, G is gravitational constant, $R_{m}$ and $M$ are radius and mass of the Moon respectively. The final energy of the neutron when it hits the detector at a radius of $R$ is given by

If the neutron leaves the surface at an angle $\theta$ and $\mu=cos\theta$, the incident angle of the particle is given by

The second step is to account for the lifetime of the neutron. There could be two cases depending on whether the kinetic energy of the particle is greater or lesser than the binding energy. When the kinetic energy of the particle is less than the binding energy ($\frac{K}{V}<1$):

where $A=\sqrt{4\frac{K}{V}\mu^{2}(1-\frac{K}{V})}$ and
$B=2\frac{K}{V}-1$

When the kinetic energy of the particle is more than the binding energy ($\frac{K}{V}>1$)

The third step is to apply effective area of the detector and account for the geometric factor correction caused due to solid angle effects at the altitude of the spacecraft,

The location of neutron emission is calculated by treating the neutron motion as orbital motion. We calculate the angle between the neutron periapsis and current location, that is, the true anomaly, at both the satellite position and the lunar surface. The difference between the true anomalies at these locations is the difference in central angle between neutron emission and the satellite location. For example, for a neutron travelling due north, the difference in true anomaly simply gives the difference in latitude between the satellite and emission location. For other neutron directions, this angle is rotated.

We calculate the true anomaly following the notation of Curtis (2005). Note that some symbols are reused from the previous section with different meanings. Defining the gravitational parameter $\mu=GM$ for the gravitational constant $G$ and mass of the Moon $M$ and the magnitude of the specific relative angular momentum $h=rv_{\perp}$ for radius $r$ and perpendicular component of the velocity $v_{\perp}$, the eccentricity is:

$v_{r}$ is the radial component of the velocity. $e$ is a constant and can be evaluated at any point in the orbit; in practice we evaluate it at the point of emission. Then the true anomaly at any radius $r$ is:

For gravitationally-bound neutrons on inward trajectories when they intersect the satellite, the true anomaly is $2\pi-\theta$ as calculated in the equation above.

Following the kinematics described in Section II.4, we get the neutron energy in a frame of reference where $\hat{z_{r}}$ is towards the local zenith at the spacecraft. To get the thermal velocity of the neutron in the rest frame of the spacecraft, we should transform the coordinate system where $\hat{z}$ is aligned with the rotation axis of the Moon and the spacecraft.

Assuming $\hat{y_{r}}$ to lie entirely in the x-y plane in the standard Moon frame, we can define the spacecraft position by $(x_{s},y_{s},z_{s})$, such that radius of the spacecraft is given by $R=\sqrt{x^{2}_{s}+y^{2}_{s}+z^{2}_{s}}$ and the component in the x-y plane is $r=\sqrt{x^{2}_{s}+y^{2}_{s}}$.

The velocity of the spacecraft is thus defined by:

We then subtract spacecraft velocity to get neutron velocity in the spacecraft rest and rotation frame, and recalculate energy and angle.