# Optimal Control of Transient Flows in Pipeline Networks with Heterogeneous
Mixtures of Hydrogen and Natural Gas

Luke S. Baker1, Saif R. Kazi2, Rodrigo B. Platte1, and Anatoly Zlotnik2 *This
study was supported by the U.S. Department of Energy’s Advanced Grid Modeling
(AGM) project “Dynamical Modeling, Estimation, and Optimal Control of
Electrical Grid-Natural Gas Transmission Systems”, as well as LANL Laboratory
Directed R&D project “Efficient Multi-scale Modeling of Clean Hydrogen
Blending in Large Natural Gas Pipelines to Reduce Carbon Emissions”. Research
conducted at Los Alamos National Laboratory is done under the auspices of the
National Nuclear Security Administration of the U.S. Department of Energy
under Contract No. 89233218CNA000001.1Luke Baker and Rodrigo Platte are with
the School of Mathematical and Statistical Sciences at Arizona State
University, Tempe, Arizona, 85281<EMAIL_ADDRESS>Kazi and
Anatoly Zlotnik are in the Applied Mathematics & Plasma Physics Group at Los
Alamos National Laboratory, Los Alamos, New Mexico, 87545;
<EMAIL_ADDRESS>
###### Abstract

We formulate a control system model for the distributed flow of mixtures of
highly heterogeneous gases through large-scale pipeline networks with time-
varying injections of constituents, withdrawals, and control actions of
compressors. This study is motivated by the proposed blending of clean
hydrogen into natural gas pipelines as an interim means to reducing end use
carbon emissions while utilizing existing infrastructure for its planned
lifetime. We reformulate the partial differential equations for gas dynamics
on pipelines and balance conditions at junctions using lumped elements to a
sparse nonlinear differential algebraic equation system. Our key advance is
modeling the mixing of constituents in time throughout the network, which
requires doubling the state space needed for a single gas and increases
numerical ill-conditioning. The reduced model is shown to be a consistent
approximation of the original system, which we use as the dynamic constraints
in a model-predictive optimal control problem for minimizing the energy
expended by applying time-varying compressor operating profiles to guarantee
time-varying delivery profiles subject to system pressure limits. The optimal
control problem is implemented after time discretization using a nonlinear
program, with validation of the results done using a transient simulation. We
demonstrate the methodology for a small test network, and discuss scalability
and potential applications.

##  I Introduction

Transportation of natural gas through networks of large-scale transmission
pipelines has been studied in steady-state [1, 2, 3, 4, 5] and transient [6,
7, 8] operations with applications to the optimal control of compressor
actuators. In steady-state, the flow of gas in the network is balanced, so
that inflows from processing plants and supply stations and outflows from
withdrawal stations sum to zero. Steady-state pipeline flows are described
using simple time-invariant algebraic equations that relate pressure drop in
the direction of flow to mass flow along each pipeline. In the transient
regime, computational complexity increases significantly because the flow in
each pipeline cannot be modeled with simple algebraic equations but rather
requires a system of nonlinear partial differential equations (PDEs) [9, 10].
Model reduction methods have been proposed to reduce the complexity of
optimizing gas flows in networks [11, 12]. Although natural gas is projected
to be a primary fuel source through the year 2050 [13], worldwide economies
have invested in transitioning from fossil fuels such as natural gas to more
sustainable resources. Hydrogen is a potential candidate, which, because it
does not produce carbon dioxide when burned, is considered to have the
potential to address climate change [14]. Natural gas pipeline operation and
management protocols may be modified to transport mixtures of natural gas and
hydrogen. Recent studies indicate that natural gas pipelines can safely and
effectively transport mixtures of up to 20% hydrogen by volume [15, 16].
However, the complexity of modeling steady-state and transient flows, and thus
designing and operating pipelines, is compounded with the injection of
hydrogen [17, 18].

Natural gas and hydrogen have significantly different physical and chemical
properties. Hydrogen is less dense than natural gas, and the speed of sound
through hydrogen is roughly four times as large as that of natural gas.
Viscosity, velocity, density, pressure, and energy of the gas mixture all vary
with varying hydrogen concentration [19, 20], and these directly affect the
transportation of the mixture [21]. Numerical simulations have been performed
to demonstrate various effects on steady-state and transient-state flows of
mixtures of hydrogen and natural gas in pipeline networks [22, 23, 24, 25, 26,
27, 28]. The method of characteristics was used in the numerical simulation of
transient flows on cycle networks with homogeneous flow mixtures [25]. Another
recent study investigates gas composition tracking using a moving grid method
and an implicit backward difference method [23]. It was shown that both
methods of tracking perform well, but the implicit difference method may lose
some finer detail in the response due to numerical diffusion. A finite element
method using COMSOL Multiphysics was also developed [26]. That study considers
the effects of hydrogen concentration on the compressibility factor of the
mixture and its relation with pressure. Moreover, there the authors
demonstrate that pressure may exceed pipeline limitations in the transient
evolution of flow and that the likelihood of this happening increases
proportionally with increasing hydrogen concentration.

In contrast to pure natural gas, few studies have examined optimization of
steady-state and transient operations of mixtures of hydrogen and natural gas
in networks. To our knowledge, there are no results on the optimal control of
transient flows of heterogeneous mixtures of gases in pipelines or networks of
pipelines. Optimal control of compressor actuators for transport of pure
natural gas typically seeks to minimize the cost of running compressors while
being subjected to PDE flow dynamics, nodal pressure and nodal flow balance
constraints, and inequality box constraints that limit the pressure throughout
the network [29]. Other formulations may use an objective function that
maximizes economic value [30]. When transients are sufficiently slow, a
friction-dominated approximation may be made [31], and this was shown to be
valid in the regime of normal pipeline operations [32]. We use friction-
dominated modeling to simplify the reduced modeling in the heterogeneous gas
setting.

In this study, we formulate a control system model for transporting
heterogeneous mixtures of gases through pipeline networks of general form, and
extend optimal control problems for gas pipeline flow to this setting. Our key
advance is modeling the mixing of constituents in time throughout the network,
which requires doubling the state space needed for a single gas and increases
numerical ill-conditioning. This enables the formulation and solution of
optimal control problems in which constituent gases may be injected at
different points in the network at varying concentrations, e.g., the addition
of 100% hydrogen at certain nodes. An algorithm is implemented to obtain
solutions, and the results are demonstrated on a small test network that
includes a cycle.

The remainder of the manuscript is organized as follows. The governing
equations for the flow of mixtures of gases in a network are presented in
Section II. In Section III, an endpoint discretization method is employed to
reduce the system of PDEs to a system of ordinary differential equations
(ODEs). There, we show that the discretization method is consistent and
results in the equations for natural gas only in the case of zero hydrogen
injection, and yields the steady-state equations when supply and withdrawal
are held constant. Section IV describes time-discretization of the optimal
control problem that yields a nonlinear program (NLP). The NLP is solved for a
test network in Section V, and we discuss applications in Section VI.

##  II Network Flow Control Formulation

We begin by defining notation that will be used in the study. A gas network is
modeled as a connected and directed graph $(\mathcal{E},\mathcal{V})$ that
consists of edges (pipelines) $\mathcal{E}=\\{1,\dots,E\\}$ and nodes
(junctions) $\mathcal{V}=\\{1,2,\dots,V\\}$, where $E$ and $V$ denote the
cardinalities of the sets. It is assumed that the nodes and edges are ordered
within their sets according to their integer labels. The symbol $k$ is
reserved for identifying edges in $\mathcal{E}$ and the symbols $i$ and $j$
are reserved for identifying nodes in $\mathcal{V}$. Supply nodes
$\mathcal{V}_{s}\subset\mathcal{V}$ and withdrawal nodes
$\mathcal{V}_{w}\subset\mathcal{V}$ are assumed to be disjoint sets that
partition $\mathcal{V}$, i.e.,
$\mathcal{V}_{s}\cup\mathcal{V}_{w}=\mathcal{V}$ and
$\mathcal{V}_{s}\cap\mathcal{V}_{w}=\emptyset$. It is assumed that supply
nodes are ordered in $\mathcal{V}$ before withdrawal nodes so that $i<j$ for
all $i\in\mathcal{V}_{s}$ and $j\in\mathcal{V}_{w}$. The graph is directed by
judiciously assigning a positive flow direction along each edge. It is assumed
that gas flows in the positive oriented direction of an edge so that the mass
flux and velocity of the gas are positive quantities. The notation $k:i\mapsto
j$ means that edge $k\in\mathcal{E}$ is directed from node $i\in\mathcal{V}$
to node $j\in\mathcal{V}$. For each node $j\in\mathcal{V}$, we define
(potentially empty) incoming and outgoing sets of pipelines by
${}_{\mapsto}j=\\{k\in\mathcal{E}|k:i\mapsto j\\}$ and
$j_{\mapsto}=\\{k\in\mathcal{E}|k:j\mapsto i\\}$, respectively.

The transportation of the mixture of hydrogen and natural gas is modeled as a
simplification of the isothermal Euler equations. For each pipe
$k\in\mathcal{E}$, the flow variables are natural gas density
$\rho_{k}^{(1)}(t,x)$, hydrogen density $\rho_{k}^{(2)}(t,x)$, and mass flux
$\varphi_{k}(t,x)$ of the mixture, with $t\in[0,T]$ and $x\in[0,\ell_{k}]$,
where $T$ denotes the time horizon and $\ell_{k}$ denotes the length of the
pipe. Assuming zero inclination and sufficiently slow transients, the flow
through edge $k$ is governed by the friction-dominated PDEs

$\displaystyle\partial_{t}\rho_{k}^{(m)}+\partial_{x}\left(\frac{\rho_{k}^{(m)}}{\rho_{k}^{(1)}+\rho_{k}^{(2)}}\varphi_{k}\right)\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!0,$ (1)
$\displaystyle\partial_{x}\left(\sigma^{2}_{1}\rho_{k}^{(1)}+\sigma^{2}_{2}\rho_{k}^{(2)}\right)\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!-\frac{\lambda_{k}}{2D_{k}}\frac{\varphi_{k}|\varphi_{k}|}{\rho_{k}^{(1)}+\rho_{k}^{(2)}},\quad$
(2)

where (1) is defined for each constituent $m=1$ and $m=2$. Superscripts “1”
and “2” attached to a gas variable are conserved for identifying natural gas
and hydrogen variables, respectively. The parameters for each pipeline
$k\in\mathcal{E}$ are diameter $D_{k}$, friction factor $\lambda_{k}$, speed
of sound through natural gas $\sigma_{1}$, and speed of sound through hydrogen
gas $\sigma_{2}$. In the above dynamic equations, the compressibility factors
of the gasses are assumed to be constants so that the equations of states are
ideally given by $p^{(m)}_{k}=\sigma^{2}_{m}\rho_{k}^{(m)}$, where
$p_{k}^{(m)}$ is the partial pressure. The summation of partial pressures
results in the equation of state
$p_{k}=(\sigma^{2}_{1}\eta_{k}^{(1)}+\sigma^{2}_{2}\eta_{k}^{(2)})\rho_{k}$,
where $p_{k}=(p^{(1)}_{k}+p^{(2)}_{k})$ is the total pressure,
$\rho_{k}=(\rho^{(1)}_{k}+\rho^{(2)}_{k})$ is the total density,
$\eta_{k}^{(1)}=\rho_{k}^{(1)}/\rho_{k}$ is the concentration of natural gas,
and $\eta_{k}^{(2)}=\rho_{k}^{(2)}/\rho_{k}$ is the concentration of hydrogen.
The concentration as defined here refers to mass fraction, so that the
volumetric fraction of hydrogen in the mixture is substantially greater.

Friction forces between the interior wall of a pipe and gas flowing through it
cause pressure to decrease in the direction of flow, as reflected in the
momentum equation (2). Compressor stations receive gas at low pressure and
reduce its volume to increase its pressure to levels required for
transportation and customer contracts. In addition to compressors, regulators
are installed to reduce the pressure of the received gas to within limits that
are compatible with lower pressure distribution systems. For convenience, we
assume that a compressor is located at the inlet and a regulator is located at
the outlet of each pipeline with respect to the prescribed positive flow
direction. For each pipeline $k\in\mathcal{E}$, compression and regulation are
modeled with time-varying multiplicative compressor ratio
$\underline{\mu}_{k}(t)\geq 1$ and regulator ratio $\overline{\mu}_{k}(t)\geq
1$, with orientations illustrated in Figure 1.

Natural gas and hydrogen are injected into the network at each supply node
$i\in\mathcal{V}_{s}$ with specified time-varying profiles of natural gas
density $\bm{s}_{i}^{(1)}(t)$ and hydrogen density $\bm{s}_{i}^{(2)}(t)$.
Alternatively, pressure
$(\bm{p}_{s})_{i}=(\sigma_{1}^{2}\bm{s}_{i}^{(1)}+\sigma_{2}^{2}\bm{s}_{i}^{(2)})$
and concentration
$\bm{\alpha}_{i}^{(m)}=\bm{s}_{i}^{(m)}/(\bm{s}_{i}^{(1)}+\bm{s}_{i}^{(2)})$
may be specified at slack nodes $i\in\mathcal{V}_{s}$ instead of the
constituent densities. Gas is withdrawn downstream at each withdrawal node
$j\in\mathcal{V}_{w}$ with specified time-varying mass flux $\bm{w}_{j}(t)$.
For $m=1,\,2$ and all $j\in\mathcal{V}_{w}$, define nodal density variables
$\bm{\rho}_{j}^{(m)}(t)$ and nodal concentration variables
$\bm{\eta}_{j}^{(m)}(t)$. All of the nodal quantities in this study are
identified with bold symbols. Inlet and outlet edge variables are defined by
attaching underlines below and overlines above the associated edge variables,
respectively. For example, $\underline{\varphi}_{k}(t)=\varphi_{k}(t,0)$ and
$\overline{\varphi}_{k}(t)=\varphi_{k}(t,\ell_{k})$. The boundary conditions
for the flow of the mixture are defined for $m=1,\,2$ by

$\displaystyle\underline{\rho}_{k}^{(m)}$ $\displaystyle=$
$\displaystyle\underline{\mu}_{k}\bm{s}_{i}^{(m)},\qquad\overline{\rho}_{k}^{(m)}=\overline{\mu}_{k}\bm{\rho}_{j}^{(m)},$
(3) $\displaystyle\underline{\rho}_{k}^{(m)}$ $\displaystyle=$
$\displaystyle\underline{\mu}_{k}\bm{\rho}_{i}^{(m)},\qquad\overline{\rho}_{k}^{(m)}=\overline{\mu}_{k}\bm{\rho}_{j}^{(m)},$
(4) $\displaystyle\bm{\eta}_{j}^{(m)}\bm{w}_{j}$ $\displaystyle=$
$\displaystyle\sum_{k\in_{\mapsto}j}\overline{\eta}_{k}^{(m)}\overline{\varphi}_{k}-\sum_{k\in
j_{\mapsto}}\underline{\eta}_{k}^{(m)}\underline{\varphi}_{k},$ (5)

where (3) is defined for $k:i\mapsto j$ with $i\in\mathcal{V}_{s}$, (4) is
defined for $k:i\mapsto j$ with $i\in\mathcal{V}_{w}$, and (5) is defined for
$j\in\mathcal{V}_{w}$. The configuration of the boundary conditions in a
pipeline segment is depicted in Fig. 1 (see Fig. 2 for a network example). The
conditions in (3)-(4) represent the effects of compression and regulation, and
the conditions in (5) represent the conservation of mass flow of each
constituent through withdrawal nodes. It is assumed that the final operating
state returns to its initial state, resulting in periodic temporal constraints

$\rho_{k}^{(m)}(0,x)=\rho_{k}^{(m)}(T,x).$ (6)

for all $k\in\mathcal{E}$ and $x\in[0,\ell_{k}]$. Periodicity in time requires
the parameters $\bm{s}_{i}^{(m)}(t)$, $\bm{w}_{j}(t)$,
$\underline{\mu}_{k}(t)$, and $\overline{\mu}_{k}(t)$ to be periodic with
period $T$. We assume that the boundary conditions are smooth, slowly-varying,
and bounded in their respective domains to ensure the existence of a smooth,
slowly-varying, bounded solution. The flow of the mixture of gases in the
network is defined by the initial-boundary value system of PDEs (1)-(6).

Figure 1: Configuration of the pipeline segment $k:i\mapsto j$ with
$i\in\mathcal{V}_{s}$ and $j\in\mathcal{V}_{w}$.

Gas network operators require pressure, compression, and regulation to be
within satisfactory limitations to ensure the safety of transportation and the
quality of gas delivered to customers. These limitations are modeled for all
$k\in\mathcal{E}$ with inequality constraints of the form

$p^{\min}_{k}\leq\sigma_{1}^{2}\rho_{k}^{(1)}+\sigma_{2}^{2}\rho_{k}^{(2)}\leq
p^{\max}_{k},\quad 1\leq\underline{\mu}_{k},\,\overline{\mu}_{k}\leq 2,$ (7)

where $p^{\min}_{k}$ and $p^{\max}_{k}$ are specified bounds on pressure for
each pipeline $k\in\mathcal{E}$. Compression $\underline{\mu}_{k}$ and
regulation $\overline{\mu}_{k}$ are the control actuators in the network that
are designed to minimize total consumption. Because pressure down-regulation
does not consume considerable energy, its proportion of total energy used for
operating a pipeline may be neglected. The total energy required for
compression is given by

$J=\sum_{k\in\mathcal{E}}\int_{0}^{T}c_{k}|\underline{\varphi}_{k}(t)|\left((\underline{\mu}_{k}(t))^{(\nu-1)/\nu}-1\right)dt,$
(8)

where $c_{k}$ is related to the efficiency of the compressor
$\underline{\mu}_{k}$ and $\nu$ is the isentropic exponent [33] (which is
assumed to be a weighted average of those of natural gas and hydrogen with
weights equal to their respective mean concentration injections). The
continuous optimal control problem is

$\begin{array}[]{ll}\text{min}&J\triangleq\text{compressor energy in
\eqref{eq:J}},\\\ \text{s.t.}&\text{dynamic constraints:
\eqref{eq:flow1}-\eqref{eq:flow2}},\\\ &\text{boundary conditions:
\eqref{eq:bc1}-\eqref{eq:bc}},\\\ &\text{temporal constraints:
\eqref{eq:ic}},\\\ &\text{inequality constraints:
\eqref{eq:ineq}}.\end{array}$ (9)

The decision variables are partial densities, mass fluxes, compressor ratios,
and regulator ratios throughout the network.

##  III Network Flow Control Discretization

The intial-boundary value system of PDEs described in the previous section
will be discretized in space to obtain an initial-value system of ODEs.
Spatial discretization is formalized by refining the graph of the gas network.
A graph refinement $(\hat{\mathcal{E}},\hat{\mathcal{V}})$ of the graph
$(\mathcal{E},\mathcal{V})$ is made by adding auxiliary nodes to $\mathcal{V}$
that subdivide the edges of $\mathcal{E}$ so that $\ell_{k}\leq\ell$ for all
$k\in\hat{\mathcal{E}}$, where $\ell\leq 10$ (km) is sufficiently small [12].
The refined graph inherits the prescribed direction of the parent graph. For
sufficiently fine network refinement, the relative difference of the flow
variables between adjacent nodes is small in magnitude by continuity of the
flow variables. We assume that the graph has been sufficiently refined and
that the hats may be omitted moving forward.

The system of ODEs is obtained by integrating the dynamic equations in (1)-(2)
along the length of each refined pipeline segment so that

$\displaystyle\int_{0}^{\ell}\partial_{t}\rho^{(m)}dx=-\int_{0}^{\ell}\partial_{x}\left(\frac{\rho^{(m)}}{\rho^{(1)}+\rho^{(2)}}\varphi\right)dx,$
$\displaystyle\int_{0}^{\ell}\partial_{x}\left(\sigma^{2}_{1}\rho^{(1)}+\sigma^{2}_{2}\rho^{(2)}\right)dx=-\frac{\lambda}{2D}\int_{0}^{\ell}\frac{\varphi|\varphi|}{\rho^{(1)}+\rho^{(2)}}dx,$

where edge subscripts have been removed for readability. The above integrals
of space derivatives are evaluated using the fundamental theorem of calculus.
The remaining integrals are evaluated by approximating pipeline density with
outlet density and pipeline flux with inlet flux. These approximations are
independent of $x$ and may be factored out of the integrals. The above
equations become

$\displaystyle\ell\dot{\overline{\rho}}^{(m)}$ $\displaystyle=$
$\displaystyle\underline{\eta}^{(m)}\underline{\varphi}-\overline{\eta}^{(m)}\overline{\varphi},$
(10)
$\displaystyle\sum_{n=1}^{2}\sigma_{n}^{2}\left(\overline{\rho}^{(n)}-\underline{\rho}^{(n)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda\ell}{2D}\frac{\underline{\varphi}\left|\underline{\varphi}\right|}{\overline{\rho}^{(1)}+\overline{\rho}^{(2)}},$
(11)

where a dot above a variable represents the time-derivative of the variable.
We now write the discretized system in matrix form. Define the $E\times E$
diagonal matrices $L$ and $K$ with diagonal entries $L_{kk}=\ell_{k}$ and
$K_{kk}=\lambda_{k}/(2D_{k})$. Define the time-varying (transposed) incidence
matrix $M$ of size $E\times V$ componentwise by

$\displaystyle M_{ki}=\begin{cases}\overline{\mu}_{k}(t),&\text{ edge
$k\in_{\mapsto}i$ enters node $i$},\\\ -\underline{\mu}_{k}(t),&\text{ edge
$k\in i_{\mapsto}$ leaves node $i$,}\\\ 0,&\text{else.}\end{cases}$ (12)

Define the $E\times r$ submatrix $M_{s}$ of $M$ by the removal of columns
$i\in\mathcal{V}_{w}$, the $E\times(V-r)$ submatrix $M_{w}$ of $M$ by the
removal of columns $i\in\mathcal{V}_{s}$, and the positive and negative parts
of $M_{w}$ by $\overline{M}_{w}$ and $\underline{M}_{w}$ so that
$M_{w}=(\overline{M}_{w}+\underline{M}_{w})/2$ and
$|M_{w}|=(\overline{M}_{w}-\underline{M}_{w})/2$, where $r$ denotes the number
of supply nodes and $|A|$ denotes the componentwise absolute value of a matrix
$A$. Define the signed matrices $Q_{w}=$ sign$(M_{w})$, $\overline{Q}_{w}=$
sign$(\overline{M}_{w})$, $\underline{Q}_{w}=$ sign$(\underline{M}_{w})$, and
similarly for $M_{s}$. These signed matrices are well-defined by the
inequalities in (7).

Define inlet and outlet edge mass flux vectors by
$\underline{\varphi}=(\underline{\varphi}_{1},\dots,\underline{\varphi}_{E})^{T}$
and
$\overline{\varphi}=(\overline{\varphi}_{1},\dots,\overline{\varphi}_{E})^{T}$,
and similarly for inlet and outlet edge concentrations. Moreover, define the
vectors
$\bm{\rho}^{(m)}=(\bm{\rho}_{r+1}^{(m)},\dots,\bm{\rho}_{V}^{(m)})^{T}$,
$\bm{\eta}^{(m)}=(\bm{\eta}_{r+1}^{(m)},\dots,\bm{\eta}_{V}^{(m)})^{T}$, and
$\bm{\alpha}^{(m)}=(\bm{\alpha}_{1}^{(m)},\dots,\bm{\alpha}_{r}^{(m)})^{T}$,
where the subscripts of the entries are indexed according to the node labels
in $\mathcal{V}$. Applying the above matrix definitions, the discretized
equations in (10)-(11) together with the boundary conditions in (3)-(5) become

$\displaystyle
L\overline{M}_{w}\dot{\bm{\rho}}^{(m)}=\underline{\eta}^{(m)}\odot\underline{\varphi}-\overline{\eta}^{(m)}\odot\overline{\varphi},$
(13)
$\displaystyle\sum_{m=1}^{2}\sigma_{m}^{2}\left(M_{w}\bm{\rho}^{(m)}+M_{s}\bm{s}^{(m)}\right)=-\frac{LK(\underline{\varphi}\odot|\underline{\varphi}|)}{\overline{M}_{w}(\bm{\rho}^{(1)}+\bm{\rho}^{(2)})},$
(14)
$\displaystyle\bm{\eta}^{(m)}\odot\bm{w}=\overline{Q}_{w}^{T}\left(\overline{\eta}^{(m)}\odot\overline{\varphi}\right)+\underline{Q}_{w}^{T}\left(\underline{\eta}^{(m)}\odot\underline{\varphi}\right),$
(15)

where $\odot$ is the Hadamard product, and the ratio of vectors on the right-
hand-side of (14) is understood to be componentwise. It is assumed that
regulators vary slowly so that the time derivative of $\overline{M}_{w}$ is
insignificant, justifying its removal from (13). Multiplying both sides of
(13) on the left by $\overline{Q}_{w}^{T}$ and using (15), we may combine (13)
and (15) to form the equation
$\overline{Q}_{w}^{T}L\overline{M}_{w}\dot{\bm{\rho}}^{(m)}=[Q_{w}^{T}(\underline{\eta}^{(m)}\odot\underline{\varphi})-\bm{\eta}^{(m)}\odot\bm{w}]$,
where we have used $Q_{w}=(\underline{Q}_{w}+\overline{Q}_{w})$. By the
definitions of supply and withdrawal concentrations, the above equations
become

$\displaystyle\overline{Q}_{w}^{T}L\overline{M}_{w}\dot{\bm{\rho}}^{(m)}=Q_{w}^{T}[(|\underline{Q}_{w}|\bm{\eta}^{(m)}+|\underline{Q}_{s}|\bm{\alpha}^{(m)})\odot\underline{\varphi}]$
$\displaystyle-\bm{\eta}^{(m)}\odot\bm{w},$ (16)
$\displaystyle\sum_{m=1}^{2}\sigma_{m}^{2}\left(M_{w}\bm{\rho}^{(m)}+M_{s}\bm{s}^{(m)}\right)=-\frac{LK(\underline{\varphi}\odot|\underline{\varphi}|)}{\overline{M}_{w}(\bm{\rho}^{(1)}+\bm{\rho}^{(2)})}.$
(17)

Periodic temporal constraints in (6) reduce to

$\displaystyle\bm{\rho}^{(m)}(0)=\bm{\rho}^{(m)}(T).$ (18)

Pressure, compression, and regulation inequality constraints in (7) reduce to

$\bm{p}^{\min}_{j}\leq\sigma_{1}^{2}\bm{\rho}_{j}^{(1)}+\sigma_{2}^{2}\bm{\rho}_{j}^{(2)}\leq\bm{p}^{\max}_{j},\quad
1\leq\underline{\mu}_{k},\,\overline{\mu}_{k}\leq 2,$ (19)

where $\bm{p}^{\min}_{j}$ and $\bm{p}^{\max}_{j}$ are specified bounds for
each node $j\in\mathcal{V}_{w}$. The reduced-model optimal control problem is
formulated as

$\begin{array}[]{ll}\text{min}&J\triangleq\text{compressor energy in
\eqref{eq:J}},\\\ \text{s.t.}&\text{dynamic constraints:
\eqref{partial1}-\eqref{partial2}},\\\ &\text{temporal constraints:
\eqref{ic2}},\\\ &\text{inequality constraints: \eqref{Disineq}}.\end{array}$
(20)

We now present a few results on the discretization method. Proposition 1 below
shows that the discretized system in (16)-(17) approaches the continuous
system in (1)-(2) in a single pipeline as the distance between adjacent nodes
of the refined pipeline approaches zero. Proposition 2 shows that the number
of density variables in (16)-(17) reduces to half this number for homogeneous
mixtures. Moreover, there we show that the discretized system reduces to the
steady-state equations in the time-invariant setting.

Proposition 1. Consider a single pipeline of length $\ell$, and refine its
graph as a chain connection of $E$ segments of uniform length
$\Delta\ell=\ell/E$, diameter $D$, and friction factor $\lambda$. Suppose the
gas mixture is supplied to the pipeline at the inlet $\mathcal{V}_{s}=\\{1\\}$
with boundary conditions as in (3) and withdrawn from only the outlet so that
$\bm{w}_{j}=0$ for $j\not=E+1$. Suppose for simplicity that there are no
compressors or regulators. Then the resulting system in (16)-(17) is a
consistent spatial discretization of (1)-(5).

Proof. The matrix $\overline{M}_{w}$ is the $E\times E$ identity matrix,
$\underline{M}_{w}$ is the $E\times E$ lower off-diagonal matrix with nonzero
entries $(\underline{M}_{w})_{n+1,n}=-1$, and $\underline{Q}_{s}$ is an
$E\times 1$ unit vector with one nonzero entry given by
$(\underline{Q}_{s})_{1}=-1$. For the intermediate segment $n:n\mapsto(n+1)$
with $2\leq n\leq E-1$, the associated dynamics in (16)-(17) are given by

$\displaystyle\dot{\bm{\rho}}^{(m)}_{n+1}+\frac{1}{\Delta\ell}\left(\bm{\eta}_{n+1}^{(m)}\underline{\varphi}_{n+1}-\bm{\eta}_{n}^{(m)}\underline{\varphi}_{n}\right)$
$\displaystyle=$ $\displaystyle\;0$
$\displaystyle\frac{1}{\Delta\ell}\sum_{m=1}^{2}\sigma_{m}^{2}\left(\bm{\rho}_{n+1}^{(m)}-\bm{\rho}_{n}^{(m)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda}{2D}\frac{\underline{\varphi}_{n}|\underline{\varphi}_{n}|}{\bm{\rho}_{n+1}^{(1)}+\bm{\rho}_{n+1}^{(2)}}.$

Taking the limit $\Delta\ell\rightarrow 0$, the above equations approach the
dynamics in (1)-(2). Similarly, as $\Delta\ell\rightarrow 0$, it can be shown
that the first and last segments of the pipe reduce to the dynamics (1)-(2)
with boundary conditions (3)-(5). $\square$

Proposition 2. The network system in (16)-(17) reduces to a system with total
density as the only density variable if the concentration is homogeneous. If
the concentration of hydrogen is zero, then the system reduces further to the
pure natural gas equations. Furthermore, if the concentration, supply
pressure, and withdrawal flux are time-invariant, yielding a time-invariant
solution, then the system in (16)-(17) reduces to the steady-state balance
laws.

Proof. Suppose that $\bm{\alpha}^{(2)}:=\bm{\alpha}$ and
$\bm{\eta}^{(2)}:=\bm{\eta}$ are constant vectors that are known. The equation
of state may be written as
$(\sigma_{1}^{2}\bm{\rho}^{(1)}+\sigma_{2}^{2}\bm{\rho}^{(2)})=\bm{a}^{2}\odot\bm{\rho}$,
where $\bm{\rho}=(\bm{\rho}^{(1)}+\bm{\rho}^{(2)})$ is the total density
vector and
$\bm{a}=(\sigma_{1}^{2}\bm{(}1-\bm{\eta})+\sigma_{2}^{2}\bm{\eta})^{1/2}$ is a
generalized sound speed vector (the square-root operation is componentwise).
Moreover,
$(\sigma_{1}^{2}\bm{s}^{(1)}+\sigma_{2}^{2}\bm{s}^{(2)})=\bm{b}^{2}\odot\bm{s}$,
where $\bm{s}=(\bm{s}^{(1)}+\bm{s}^{(2)})$ and
$\bm{b}=(\sigma_{1}^{2}\bm{(}1-\bm{\alpha})+\sigma_{2}^{2}\bm{\alpha})^{1/2}$.
Because $\bm{a}$ and $\bm{b}$ are constant vectors, the superposition of (16),
for $m=1,\,2$, results in

$\displaystyle\overline{Q}_{w}^{T}L\overline{M}_{w}\dot{\bm{\rho}}$
$\displaystyle=Q_{w}^{T}\underline{\varphi}-\bm{w},$ (21) $\displaystyle
M_{w}\left(\bm{a}^{2}\odot\bm{\rho}\right)+M_{s}\left(\bm{b}^{2}\odot\bm{s}\right)$
$\displaystyle=-\frac{LK(\underline{\varphi}\odot|\underline{\varphi}|)}{\overline{M}_{w}\bm{\rho}}.\quad$
(22)

Partial densities may be determined using
$\bm{\rho}^{(m)}=\bm{\eta}^{(m)}\odot\bm{\rho}$. This shows that, under the
assumption of constant concentration, the reduced system in (21)-(22) may be
used to determine the solution to (16)-(17). If $\bm{\alpha}=0$ and
$\bm{\eta}=0$, then the above equations reduce to the single gas endpoint
discretization method with $\bm{a}=\bm{b}=\sigma_{1}$ [34]. If $\bm{w}$,
$\bm{\eta}$, $\bm{s}$, and $\bm{\rho}$ are constant, then the system (21)-(22)
reduces to the Weymouth equations for a mixture of gases [5]. $\square$

##  IV Implementation

The optimal control problem in (20) may be expressed as

$\displaystyle\min\qquad$ $\displaystyle\int_{0}^{T}\mathcal{F}(x(t),u(t))dt$
(23a) s.t. $\displaystyle F\frac{d}{dt}[R^{(m)}x(t)]=f^{(m)}(x(t),u(t),p(t)),$
(23b) $\displaystyle e(x(T),x(0),u(T),u(0))=0,$ (23c) $\displaystyle
q(x(t),u(t))\geq 0,$ (23d)

where $x=(\bm{\rho}^{(1)},\bm{\rho}^{(2)},\underline{\varphi})^{T}$ is the
state, $u=(\\{\underline{\mu}_{k},\overline{\mu}_{k}\\})^{T}$ is the control,
and $p=(\bm{s}^{(1)},\bm{s}^{(2)},\bm{w})^{T}$ is a fixed (potentially time-
varying) vector of parameters. The function $e$ represents temporal periodic
constraints, $q$ represents inequality constraints, $R^{(m)}$ is a row
selector matrix that maps $x$ into $\bm{\rho}^{(m)}$, and
$F=\overline{Q}_{w}^{T}L\overline{M}_{w}$. A nonlinear program is obtained by
discretizing the time interval $[0,T)$ into $N$ subintervals with equally-
spaced collocation points $t_{n}=(n-1)T/N$ for $n=1,\dots,N$.

The vector-valued functions $x(t)$, $u(t)$, and $p(t)$ are interpolated with
piecewise-linear vector-valued functions $\hat{x}(t)$, $\hat{u}(t)$, and
$\hat{p}(t)$, respectively. For $n=1,\dots,N$, the function $\hat{x}(t)$ is
defined for $t\in[t_{n},t_{n+1})$ by

$\hat{x}(t)=x(t_{n})+\frac{x(t_{n+1})-x(t_{n})}{T/N}(t-t_{n}),$

where $x(t_{N+1})=x(t_{1})$ by the assumption of periodicity. The functions
$\hat{u}(t)$ and $\hat{p}(t)$ are defined similarly. The integral in the
objective function is approximated using the trapezoidal rule, resulting in

$\int_{0}^{T}\mathcal{F}(x(t),u(t))dt\approx\sum_{n=1}^{N}\frac{T}{N}\mathcal{F}(\hat{x}(t_{n}),\hat{u}(t_{n})).$

The time derivative of $R^{(m)}x(t)$ is approximated using the finite
difference $d/dt[R^{(m)}x(t)]\approx N/T(R^{(m)}x(t_{n+1})-R^{(m)}x(t_{n}))$.
Define the $N\times N$ differentiation matrix $D$ componentwise by
$D_{n,n}=-T/N$, $D_{n,n+1}=T/N$, and $D_{N,1}=T/N$. The differentiation matrix
includes a periodic constraint to the vector on which it operates. Define the
stacked state vector $X=(\hat{x}(t_{1}),\dots,\hat{x}(t_{N}))^{T}$, input
$U=(\hat{u}(t_{1}),\dots,\hat{u}(t_{N}))^{T}$, parameter
$P=(\hat{p}(t_{1}),\dots,\hat{p}(t_{N}))^{T}$, and extend $D$, $F$, $R^{(m)}$,
and all of the other matrices using Kronecker products with identity matrices
of appropriate dimensions. Then the nonlinear program may be written as

$\displaystyle\min\qquad$
$\displaystyle\sum_{n=1}^{N}\frac{T}{N}\mathcal{F}(X_{n}X,U_{n}U)$ (24a) s.t.
$\displaystyle FDR^{(m)}X_{n}X=f^{(m)}(X_{n}X,U_{n}U,P_{n}P),$ $\displaystyle
m=1,2,\;n=1,\dots,N,$ (24b) $\displaystyle q(X_{n}X,U_{n}U)\geq
0,\;n=1,\dots,N,$ (24c)

where $X_{n}$, $U_{n}$, and $P_{n}$ are row selector matrices that satisfy
$X_{n}X=\hat{x}(t_{n})$, $U_{n}U=\hat{u}(t_{n})$, and $P_{n}P=\hat{p}(t_{n})$
for $n=1,\dots,N$.

Solution of the NLP in (24) is implemented in Matlab with the interior-point
algorithm using the function `fmincon`, and is evaluated on a MacBook Air
8‑core CPU with 8GB of unified memory. The gradient of the objective and
Jacobian of the constraints are supplied to the function for improved
performance. The Hessian of the Lagrangian function is set to the default
finite-difference approximation. Optimal control of compression and regulation
is obtained from the optimal solution $U=U^{*}$. The optimized time-series for
compressors and regulators, and any specified parameters, are linearly
interpolated to provide control functions to the ODEs (16)-(17). The system is
simulated in Matlab using the function `ode15s` for validation of the solution
and an improved prediction of pressure and mass flux. The steady-state
solution is used as the starting point for optimization, and the initial state
of the optimal solution is used as the initial condition for simulation. In
the following, we distinguish between the solution of the optimization problem
(24) and the solution of the ODEs (16)-(17) that are driven by optimal
compression and regulation. We consider this comparison an important
validation of the presented optimal control scheme, because feasibility of the
coarsly discretized physical system in the optimal solution of the NLP (24)
does not necessarily guarantee that the control solution obtained by solving
(24) results in that same physical solution in a simulation with controlled
error. For the case studies in the next section, the two solutions are
compared using the average $L^{2}$ norm of the relative difference given by

$\frac{1}{\hat{E}}\sum_{k\in\hat{\mathcal{E}}}\left(\frac{1}{T}\int_{0}^{T}\left(2\frac{\underline{\varphi}_{k}(t)-\underline{\phi}_{k}(t)}{\underline{\varphi}_{k}(t)+\underline{\phi}_{k}(t)}\right)^{2}dt\right)^{1/2}\times
100,$ (25)

where $\underline{\varphi}_{k}$ is the optimized flux in edge
$k\in\hat{\mathcal{E}}$, $\underline{\phi}_{k}$ is the simulated flux, and
$\hat{E}$ is the cardinality of $\hat{\mathcal{E}}$. In addition, the maximum
absolute relative difference is also documented as

$\max_{k\in\hat{\mathcal{E}}}\left(\max_{t\in[0,T]}\left|2\frac{\underline{\varphi}_{k}(t)-\underline{\phi}_{k}(t)}{\underline{\varphi}_{k}(t)+\underline{\phi}_{k}(t)}\right|\right)\times
100.$ (26)

Similar metrics are used for the difference between optimized and simulated
pressure trajectories.

Figure 2: Top: Configuration of the network. Pipeline dimensions: blue to
purple = 50 km, purple to cyan = 30 km, purple to green = 20 km, and green to
cyan is 30 km. Bottom: Withdrawal flow profiles, color-coded to correspond to
associated nodes. The network has two compressors, indicated by triangles,
located at the start and end of the pipe directed from blue to purple.

##  V Case Study

The optimal control algorithm is demonstrated on a cyclic network whose
configuration and dimensions are shown on the top of Fig. 2. We create a
refined network with a uniform discretization length of $\ell_{k}=10$ km for
all $k\in\hat{\mathcal{E}}$. The diameters and friction factors of the refined
pipelines are uniform and equal to $D_{k}=0.5$ m and $\lambda_{k}=0.011$ for
all $k\in\hat{\mathcal{E}}$. The speeds of sounds of the gases are
$\sigma_{1}=338.38$ m/s and $\sigma_{2}=4\sigma_{1}$. We use $N=20$ time steps
with $\nu=1.28$ and compressor efficiency values $c_{1}=c_{5}=D_{k}^{4}/T$ in
(8). The minimum and maximum pressures in (19) are $\bm{p}^{\min}_{j}=5$ MPa
and $\bm{p}^{\max}_{j}=12$ MPa for all $j\in\hat{\mathcal{V}}_{w}$.
Discretization results in 780 optimization variables, 740 equality
constraints, and 520 inequality constraints in the NLP (24).

Figure 3: Left column: Hydrogen concentration profiles at the supply node.
Right column: Corresponding optimal compressor responses. Figure 4: Left
column: Optimized pressure at refined withdrawal nodes driven by the
respective concentration profiles from Fig. 3. Right column: Simulated
pressure at refined withdrawal nodes driven by inlet concentration and optimal
compression profiles from Fig. 3. The average $L^{2}$ norm of the relative
difference in pressures, using the metric (25), are approximately 0.769%,
0.770%, and 0.769% for the top, middle, and bottom rows. The maximum relative
difference in pressures, using the metric (26), are approximately 2.154%,
2.038%, and 1.971% for the top, middle, and bottom rows. The combined
computational times for optimization and simulation are approximately 10.0 s,
8.0 s, and 7.7 s for the top, middle, and bottom rows, respectively.

The purple, green, and cyan nodes in the network graph in Figure 2 represent
stations where gas is withdrawn with color-coordinated flow profiles depicted
on the bottom of Fig. 2. The red and yellow triangles represent two compressor
stations whose time-dependent operations are optimized in a model-predictive
manner. The blue node is the supply station for a mixture of natural gas and
hydrogen with a fixed pressure of $(\bm{p}_{s})_{\text{blue}}=5$ MPa that is
immediately boosted by the red compressor station. We demonstrate three
solutions for this network, each of which are subject to the same above
boundary conditions but differ in the injected concentration of hydrogen at
the supply node. The left column of Fig. 3 depicts the specified hydrogen
concentration profiles at the supply node for the three solutions and the
right column shows the associated results for optimal compression ratios of
the two color-coordinated compressor stations. The total compressor energy
values in as defined in equation (8) are $J=0.787$, $J=0.824$, and $J=0.860$
(non-dimensionalized units) corresponding to the optimal compression ratios
given in Fig. 3 from top to bottom, respectively. Figures 4 and 5 depict
pressure and mass flux solutions, respectively, where the left-hand-side
columns show the optimized solutions and the right-hand-side columns show the
validating simulation. The results demonstrate that minor variations in
hydrogen concentration may have substantial effects on pressure and compressor
activity. In particular, the pressure trajectories in Figure 4 increase by
approximately 2 MPa from $t=0$ to $t=8$ hours for a fixed 10% hydrogen
injection, while the change is approximately 3.5 MPa when the hydrogen
injection concentration slowly varies from 8% to 12% over the 24 hour time
horizon. Observe that because this concentration is defined as the mass
fraction, the volume fraction here is actually over 30%, and exhibits
variations of up to 10%. A detailed analysis of how hydrogen blending impacts
energy transport capacity of gas pipeline networks is outside the scope of
this study, in which we focus on modeling flow dynamics and demonstrating our
optimal control approach.

Figure 5: Left column: Optimized mass flux in refined edges driven by the
respective concentration profiles from Fig. 3. Right column: Simulated mass
flux in refined edges driven by inlet concentration and optimal compression
profiles from Fig. 3. The $L^{2}$ relative difference metric values using (25)
are approximately 3.994%, 4.608%, and 5.258% for the top, middle, and bottom
rows. The maximum relative difference values using metric (26) are
approximately 12.967%, 16.713%, and 21.509%, respectively.

##  VI Conclusions

We synthesized a control system model for the distributed flow of mixtures of
two gases with different physical properties through large-scale pipeline
networks with time-varying injections, withdrawals, and control actions of
compressors. The motivation is to develop analysis methods to evaluate recent
proposals for blending of clean hydrogen into natural gas pipelines as an
interim means for carbon emissions reduction that allows utilization of
existing infrastructure for its planned lifetime [35]. The partial
differential equations for gas dynamics on pipelines and balance conditions at
junctions are approximated using lumped elements to a sparse nonlinear
differential algebraic equation system. In contrast to previously published
results, we model the mixing of constituents in time throughout the network.
The optimal control technique we develop is able to rapidly produce validated
solutions, even though representing dynamics of the gas mixture requires
doubling the state space with respect to models for a single gas, and worsens
numerical conditioning. We show that the reduced model is a consistent
approximation of the original system, use it as the dynamic constraints in a
model-predictive optimal control method for minimizing the energy expended by
applying time-varying compressor operating profiles to guarantee time-varying
delivery profiles subject to system pressure limits. The optimal control
problem is implemented after time discretization using a nonlinear program,
with validation of the results done using a transient simulation.

The developed control system model and computational optimal control scheme
can be used to solve a variety of problem formulations for gas transport
networks. The objective function could be modified to reflect the economic
value of pipeline transport, in terms of natural gas and hydrogen flow
provided by suppliers, and energy received by consumers. Including a price of
carbon emissions mitigation due to replacement of natural gas with hydrogen
could indicate optimal locations for integrating hydrogen supplies.

## References

  * [1] P. L. Wong and R. Larson, “Optimization of natural-gas pipeline systems via dynamic programming,” _IEEE Transactions on Automatic Control_ , vol. 13, no. 5, pp. 475–481, 1968.
  * [2] P. B. Percell and M. J. Ryan, “Steady state optimization of gas pipeline network operation,” in _PSIG annual meeting_. OnePetro, 1987.
  * [3] S. Wu, R. Z. Rios-Mercado, E. A. Boyd, and L. R. Scott, “Model relaxations for the fuel cost minimization of steady-state gas pipeline networks,” _Mathematical and Computer Modelling_ , vol. 31, no. 2-3, pp. 197–220, 2000\.
  * [4] D. De Wolf and Y. Smeers, “The gas transmission problem solved by an extension of the simplex algorithm,” _Management Science_ , vol. 46, no. 11, pp. 1454–1465, 2000.
  * [5] S. Misra, M. W. Fisher, S. Backhaus, R. Bent, M. Chertkov, and F. Pan, “Optimal compression in natural gas networks: A geometric programming approach,” _IEEE Transactions on Control of Network Systems_ , vol. 2, no. 1, pp. 47–56, 2014.
  * [6] M. Abbaspour and K. S. Chapman, “Nonisothermal Transient Flow in Natural Gas Pipeline,” _Journal of Applied Mechanics_ , vol. 75, no. 3, 05 2008, 031018\. [Online]. Available: https://doi.org/10.1115/1.2840046
  * [7] A. Osiadacz, “Simulation of transient gas flows in networks,” _International Journal for Numerical Methods in Fluids_ , vol. 4, no. 1, pp. 13–24, 1984.
  * [8] A. Zlotnik, M. Chertkov, and S. Backhaus, “Optimal control of transient flow in natural gas networks,” in _54th IEEE Conference on Decision and Control_. IEEE, 2015, pp. 4563–4570.
  * [9] A. J. Osiadacz and M. Chaczykowski, “Comparison of isothermal and non-isothermal pipeline gas flow models,” _Chemical Engineering Journal_ , vol. 81, no. 1-3, pp. 41–51, 2001.
  * [10] A. R. D. Thorley and C. H. Tiley, “Unsteady and transient flow of compressible fluids in pipelines—a review of theoretical and some experimental studies,” _International journal of heat and fluid flow_ , vol. 8, no. 1, pp. 3–15, 1987.
  * [11] R. Z. Ríos-Mercado _et al._ , “A reduction technique for natural gas transmission network optimization problems,” _Annals of Operations Research_ , vol. 117, no. 1, pp. 217–234, 2002.
  * [12] S. Grundel, N. Hornung, B. Klaassen, P. Benner, and T. Clees, “Computing surrogates for gas network simulation using model order reduction,” in _Surrogate-Based Modeling and Optimization_. Springer, 2013, pp. 189–212.
  * [13] S. Nalley and A. LaRose, “Annual energy outlook 2022,” _U.S. Energy Information Agency_ , 2022.
  * [14] B. L. Salvi and K. A. Subramanian, “Sustainable development of road transportation sector using hydrogen energy system,” _Renewable and Sustainable Energy Reviews_ , vol. 51, pp. 1132–1155, 2015.
  * [15] M. Götz, J. Lefebvre, F. Mörs, A. M. Koch, F. Graf, S. Bajohr, R. Reimert, and T. Kolb, “Renewable power-to-gas: A technological and economic review,” _Renewable energy_ , vol. 85, pp. 1371–1390, 2016.
  * [16] M. Ozturk and I. Dincer, “A comprehensive review on power-to-gas with hydrogen options for cleaner applications,” _International Journal of Hydrogen Energy_ , vol. 46, no. 62, pp. 31 511–31 522, 2021.
  * [17] T. van der Hoeven, _Math in gas and the art of linearization_. Energy Delta Institute Groningen, The Netherlands, 2004.
  * [18] M. W. Melaina, O. Antonia, and M. Penev, “Blending hydrogen into natural gas pipeline networks: A review of key issues,” _National Renewable Energy Laboratory_ , 2013.
  * [19] A. A. Abd, S. Z. Naji, T. C. Thian, and M. R. Othman, “Evaluation of hydrogen concentration effect on the natural gas properties and flow performance,” _International Journal of Hydrogen Energy_ , vol. 46, no. 1, pp. 974–983, 2021.
  * [20] T. Blacharski _et al._ , “The effect of hydrogen transported through gas pipelines on the performance of natural gas grid,” _AGH Drilling, Oil, Gas_ , vol. 33, no. 2, 2016.
  * [21] Ł. Zabrzeski, P. Janusz, K. Liszka, M. Łaciak, and A. Szurlej, “Hydrogen-natural gas mixture compression in case of transporting through high-pressure gas pipelines,” _IOP Conference Series: Earth and Environmental Science_ , vol. 214, p. 012137, 2019. [Online]. Available: https://doi.org/10.1088/1755-1315/214/1/012137
  * [22] F. E. Uilhoorn, “Dynamic behaviour of non-isothermal compressible natural gases mixed with hydrogen in pipelines,” _International journal of hydrogen energy_ , vol. 34, no. 16, pp. 6722–6729, 2009.
  * [23] M. Chaczykowski, F. Sund, P. Zarodkiewicz, and S. M. Hope, “Gas composition tracking in transient pipeline flow,” _Journal of Natural Gas Science and Engineering_ , vol. 55, pp. 321–330, 2018.
  * [24] G. Guandalini _et al._ , “Dynamic modeling of natural gas quality within transport pipelines in presence of hydrogen injections,” _Applied Energy_ , vol. 185, pp. 1712–1723, 2017.
  * [25] S. Elaoud, Z. Hafsi, and L. Hadj-Taieb, “Numerical modelling of hydrogen-natural gas mixtures flows in looped networks,” _Journal of Petroleum Science and Engineering_ , vol. 159, pp. 532–541, 2017.
  * [26] Z. Hafsi, S. Elaoud, and M. Mishra, “A computational modelling of natural gas flow in looped network: Effect of upstream hydrogen injection on the structural integrity of gas pipelines,” _Journal of Natural Gas Science and Engineering_ , vol. 64, pp. 107–117, 2019.
  * [27] D. Fan, J. Gong, S. Zhang, G. Shi, Q. Kang, Y. Xiao, and C. Wu, “A transient composition tracking method for natural gas pipe networks,” _Energy_ , vol. 215, p. 119131, 2021.
  * [28] B. G. Agaie _et al._ , “Reduced-order modelling for high-pressure transient flow of hydrogen-natural gas mixture,” _The European Physical Journal Plus_ , vol. 132, no. 5, pp. 1–16, 2017.
  * [29] H. H. Rachford and R. G. Carter, “Optimizing pipeline control in transient gas flow,” in _PSIG annual meeting_. OnePetro, 2000.
  * [30] A. Zlotnik, K. Sundar, A. M. Rudkevich, A. Beylin, and X. Li, “Optimal control for scheduling and pricing intra-day natural gas transport on pipeline networks,” in _58th Conference on Decision and Control_. IEEE, 2019, pp. 4887–4884.
  * [31] K. Sundar and A. Zlotnik, “State and parameter estimation for natural gas pipeline networks using transient state data,” _IEEE Transactions on Control Systems Technology_ , vol. 27, no. 5, pp. 2110–2124, 2019.
  * [32] S. Misra, M. Vuffray, and A. Zlotnik, “Monotonicity properties of physical network flows and application to robust optimal allocation,” _Proceedings of the IEEE_ , vol. 108, no. 9, pp. 1558–1579, 2020.
  * [33] I. Marić, A. Galovic, and T. Smuc, “Calculation of natural gas isentropic exponent,” _Flow Measurement and Instrumentation_ , vol. 16, pp. 13–20, 03 2005.
  * [34] C. Himpe, S. Grundel, and P. Benner, “Model order reduction for gas and energy networks,” _Journal of Mathematics in Industry_ , vol. 11, no. 1, pp. 1–46, 2021.
  * [35] A. S. K. Raju, A. Martinez-Moralez, O. Lever, and L. Asiedu-Akrofi, “Hydrogen blending impacts study,” _California Public Utilities Commission_ , 2022\. [Online]. Available: http://docs.cpuc.ca.gov/PublishedDocs/Efile/G000/M493/K760/ 493760600.PDF


# Measurement of the Free Neutron Lifetime using the Neutron Spectrometer on
NASA’s Lunar Prospector Mission

Jack T. Wilson<EMAIL_ADDRESS>David J. Lawrence  Patrick N. Peplowski
The Johns Hopkins Applied Physics Laboratory,
11101 Johns Hopkins Road,
Laurel, Md. 20723, USA.  Vincent R. Eke  Jacob A. Kegerreis  Institute for
Computational Cosmology,
Durham University, South Road,
Durham DH1 3LE, UK.

###### Abstract

We use data from the Lunar Prospector Neutron Spectrometer to make the second
space-based measurement of the free neutron lifetime finding $\tau_{n}=887\pm
14_{\text{stat}}{\>{}^{+7}_{-3\>\text{syst}}}$ s, which is within 1$\sigma$ of
the accepted value. This measurement expands the range of planetary bodies
where the neutron lifetime has been quantified from space, and by extending
the modeling to account for non-uniform elemental composition, we mitigated a
significant source of systematic uncertainty on the previous space-based
lifetime measurement. This modeling moves space-based neutron lifetime
measurement towards the ultimate goal of reducing the magnitude of the
systematics on a future space-measurement to the level of those seen in
laboratory-based experiments.

##  I Introduction

The neutron mean lifetime $\tau_{n}$ is an important parameter for the weak
interaction, and more precise knowledge of its value is needed in particle
physics, nuclear physics, and cosmology [1, 2]. Specifically, $\tau_{n}$ is a
key input to calculations of primordial helium abundance, and the
uncertainties in these predictions are presently dominated by those on
$\tau_{n}$ [3]. Primordial nucleosynthesis is one of the major lines of
evidence for the big bang, along with the observed expansion of the universe
and cosmic microwave background. Additionally, in combination with beta-decay
correlations, $\tau_{n}$ is used to test the unitarity of the
Cabibbo–Kobayashi–Maskawa (CKM) matrix, which is one of the most important
low-energy tests of the standard model [1].

Currently, there are two competing values for $\tau_{n}$ based on the results
of two different classes of laboratory experiments. The ‘bottle’ experiments
involve counting the number of neutrons that survive within a material,
magnetic, and/or gravitational trap as a function of time. The ‘beam’
experiments involve measuring the rate of production of $\beta$-decay products
in a neutron beam passing through a trapping region. The average beam
measurement ${\tau_{n}^{\textrm{beam}}=888\pm 2}$ s differs by about 4$\sigma$
from the more precise ultra-cold trapped neutron average
${\tau_{n}^{\textrm{bottle}}=879.4\pm 0.6}$ s. This discrepancy, which has
persisted for 15 years, has become known as the ‘neutron lifetime puzzle’. The
most likely explanation for the discrepancy is the presence of an unaccounted
for systematic error in one, or both, classes of experiment. However, given
the direction of the disagreement, a physical explanation is possible where
the neutron decays to unobserved particles outside of the standard model with
a branching fraction of approximately 1% [e.g., 4, 5].

The feasibility of a third technique, measuring $\tau_{n}$ from space, was
recently demonstrated using data taken by NASA’s MESSENGER spacecraft during
its flybys of Venus and Mercury [6]. The opportunity to make a space-based
measurement of $\tau_{n}$ is made possible by the fact that planetary surfaces
are constantly bombarded by galactic cosmic rays (GCRs) that collide with
atomic nuclei, which leads to the liberation of large numbers of high-energy
neutrons. These neutrons have their energy moderated downwards in subsequent
collisions with near-surface nuclei. Some fraction of the neutrons undergo a
sufficiently large number of collisions that they approach thermal equilibrium
with the atmosphere or solid surface. These thermal neutrons typically have
velocities on the order of a few km s-1. Thus, their time of flight between
emission and detection by a spacecraft hundreds to thousands of kilometers
above the surface is of order $\tau_{n}$. The neutron flux that ultimately
escapes from a planetary body into space is characteristic of that planet’s
elemental composition on depth scales of order the neutron mean free path,
which is typically $\sim$ 10 cm in silicate rocks [7]. As measuring planets’
surface compositions often forms a major goal of planetary missions, several
neutron spectrometers have been included on such missions to learn about
planetary surface composition [8, 9, 10, 11].

Although MESSENGER spent several years in orbit around Mercury, uncertainty in
that planet’s surface composition greatly constrained the accuracy of the
resulting measurement, which had to rely on a small amount of data taken
during the spacecraft’s flyby of Venus [6]. In this paper we consider data
taken by NASA’s Lunar Prospector (LP) mission during its initial highly
elliptical orbits of the Moon [12]. Unlike Mercury, the Moon’s surface
composition is well understood and can be included in the modeling of the
detected neutron count rate. Thus, our results demonstrate the practicality of
reducing this major systematic uncertainty in the space-based approach to
measuring $\tau_{n}$.

##  II  LP elliptical orbital data

LP arrived in orbit at the Moon on 11 January 1998 and spent 18 months
gathering data in mapping orbits that were close to circular. The LP neutron
spectrometer (NS) was designed to explore the Moon’s surface composition with
a focus on searching for water ice at the poles and mapping the Moon’s
elemental composition [13]. The NS consisted of two cylindrical gas
proportional counters 5.7 cm in diameter and 20 cm in length filled with 10
atmospheres of 3He. One of the detectors was covered in 0.63 mm of Cd to
shield it from thermal neutrons and the other with 0.63 mm of Sn, which has a
similar atomic number to Cd but a much smaller neutron capture cross section
and a negligible effect on the number of neutrons detected [13]. The
spacecraft body was cylindrical with the detectors mounted on the end of a 2.5
m boom perpendicular to the cylinder’s long axis to reduce the spacecraft’s
contribution to the measured signal. The spacecraft was spinning about its
long axis at 12 revolutions per minute and the integration time for a single
observation was 32 s. Thus, the response of the detector was an angle-averaged
one with respect to rotation about LP’s long axis.

Before LP entered its circular mapping orbits it carried out a set of nine
highly elliptical orbits [12]. We used the variation in thermal neutron flux
measured during these elliptical orbits, as determined by the difference in
the two 3He gas proportional counters, to infer $\tau_{n}$ by minimizing a
chi-squared comparison to Monto Carlo simulated fluxes. The spacecraft
altitude, along with the measured and modeled neutron count rates for these
orbits, is shown in Fig. 1.

Figure 1: (a) LP altitude during the initial elliptical orbits. (b) The
measured count rates in the Cd- and Sn-covered NS detectors. (c) The modeled
count rates in the Cd- and Sn- covered detectors as described in section III.
(d) The residuals of the thermal neutron counts (i.e., the difference between
the Sn- and Cd- covered detectors) normalized by the uncertainty due to
counting statistics. The grey regions show the data used to produce the final
result based on the cuts described in section IV. Time on the x-axis is
measured from January 1st, 1998. Figure 2: Low-altitude 32-s (solid) and
mission averaged (dotted) LPNS spectra from both the Cd- and Sn- covered 3He
gas proportional counters. The grey shaded region shows those channels
involved in determining the measured neutron counts. The dashed curves show
the high-altitude background, measured during LP’s approach to the Moon when
the spacecraft altitude was $>1.5\times 10^{4}$ km.

The observations of each 3He detector were recorded as 32-channel pulse height
spectra every 32 s (Fig. 2). The measured spectra include the 764 keV energy
deposition peak from the
${n+{}^{3}\textrm{He}\rightarrow{}^{3}\textrm{H}+{}^{1}\textrm{H}}$ reaction,
which is the signal of interest. Due to the spacecraft design the continuum
background originating from the interaction of GCRs with the spacecraft and
detector is small (Fig. 2). The conversion of recorded spectra to count rates
involved removing this background and then summing over the channels that
measure the neutron absorption peak. Due to the short duration of the
elliptical phase of the mission, GCR and gain variation corrections were not
necessary, thus the data reduction is much simplified compared with analyses
of the full LPNS mission dataset [e.g., 12].

The mean LPNS background spectrum was estimated by averaging the spectra taken
during the spacecraft’s initial approach to the Moon when its altitude was
greater than $1.5\times 10^{4}$ km. Background-subtracted spectra were then
produced by subtracting a solid-angle-scaled high-altitude background
$b^{\prime}$ from each individual spectrum in the time series, i.e.,

$\displaystyle b^{\prime}$ $\displaystyle=fb,$ $\displaystyle f$
$\displaystyle=\frac{1+\cos(\theta)}{2},$ $\displaystyle\theta$
$\displaystyle=\sin^{-1}\left(\frac{R_{\leftmoon}}{R_{\leftmoon}+h}\right),$

where $b$ is the high-altitude background spectrum, $R_{\leftmoon}$ is the
mean radius of the Moon and $h$ the spacecraft’s altitude above the Moon’s
surface. This solid angle scaling is a consequence of the fact that the
background is largely due to the interaction of GCRs with the spacecraft,
which are obscured by the Moon when the spacecraft is at low altitudes. To
calculate the neutron count rate, the background-subtracted spectra were
summed over the channels containing the neutron-absorption peak (shown with
the grey shaded region in Fig. 2) before being divided by the observation
period. There is a small shift in the location of the peak in the pulse height
spectra associated with temperature changes in the detectors, however the
summation region is sufficiently large that these variations have no effect on
the results. Uncertainties are those resulting from the Poisson statistics of
the observed spectra. Data from both of the 3He detectors are shown in Fig.
1(b).

##  III Neutron production and transport modeling

As in Wilson  _et al._ [6] we made use of a comparison between the measured
data and the result of models of neutron production, transport, and detection
to estimate $\tau_{n}$. The modeled count rates were determined using a
combination of three separate calculations. First, we used the particle
transport code MCNPX [14] to model the neutron production following GCR
interaction with the Moon’s surface and the flux due to the escape of these
neutrons for the set of compositions shown in Table 1, with the Moon taken to
be a uniform sphere of radius 1738 km. The GCR flux was assumed to be
isotropic, with an energy distribution described by the parameterization of
Masarik and Reedy [15] with solar modulation parameter equal to 550 MV. Due to
the low energy of the incident protons (E $\lesssim$ 10 MeV) the contribution
of the solar wind to neutron production is negligible outside of brief solar
energetic particle events [16], none of which occurred during LP’s orbital
insertion.

In the second step of the modeling, an extended version of the formalism of
Feldman  _et al._ [17] was used to calculate the flux at the spacecraft
altitude by analytically extending the surface fluxes calculated in the
previous step. Here, we have modified the approach to include the effect of
variations in the planet’s surface composition on the neutron flux. The Moon’s
surface is conventionally taken to consist of three major compositional
terranes: the Procellarum KREEP terrane (PKT), a high-Th, magnesium and iron
rich province on the lunar nearside; the South Pole Aitken (SPA) terrane, an
impact feature with intermediate-Th abundances on the lunar farside that may
include material from the upper mantle and lower crust; and the Feldspathic
Highland terrane (FHT), an anorthositic region exposing the result of early
lunar crustal differentiation [18]. To capture the range of thermal neutron
variation across the lunar surface, we add two supplemental regions: the
nearly Pure Anorthosite (PAN) terrane[19]; and the non-PKT maria (nPKT), a
mafic but less Th-rich companion to the PKT. The PAN region has a thermal
neutron flux significantly higher than the surrounding regions due to its
deficit in rare earth elements (REE) with large neutron absorption cross
sections. We used the LP low-altitude (30 km average) thermal neutron data to
define its extent. The nPKT has a thermal neutron flux intermediate between
the PKT and highlands due to its high Fe but low REE composition and its
boundaries were drawn based on LP low-altitude Th gamma-ray line and fast
neutron data, using Th as a proxy for REE [20] and fast neutrons as
representative of Fe. The extent of each of the regions is shown in Fig. 3.

The compositions of each of the regions to be used in the models were based on
our current best understanding of lunar near-surface bulk composition as
described in several studies of LP NS and GRS data [20, 21, 22, 23]. The major
and radioactive element compositions of each region were determined by
averaging, over each region, the elemental abundances derived in Prettyman
_et al._ [21], which analyzed LP GRS spectra taken during the initial circular
mapping orbit when the spacecraft had an average altitude of 100 km. For the
rare earth elements Sm and Gd we used the abundances derived in Elphic  _et
al._ [20], which are based on thermal neutron measurements taken during the
low altitude phase of the LP mission when the spacecraft had an average
altitude around 30 km. The H abundances were derived from Lawrence  _et al._
[23]. Although these compositions provide our best estimate of the elemental
abundances of the lunar surface the elemental abundances determined using this
method do not yield macroscopic neutron absorption cross sections $\Sigma_{a}$
that match the average values of $\Sigma_{a}$ measured for each region [22].
As our measurement is determined by the thermal neutron flux, having lunar
soils with the correct $\Sigma_{a}$ in our models is important. Therefore, the
final step in setting the soil compositions was to alter the elemental
abundances so that $\Sigma_{a}$ in the simulated soils matched the measured
values. We chose to modify the relative abundances of the most abundant
element, oxygen, while maintaining the relative abundances between each of the
other elments. The size of these final modifications to the oxygen abundances
is small, ranging from a 0.3% decrease for the FAN composition to a 5%
decrease for PKT soil. The compositions used in this study are shown in Table
1.

Figure 3: (a) The five regions with distinct compositions included in the
neutron count-rate models: Procellarum KREEP terrane (PKT), South Pole Aitken
Terrane (SPAT), Feldspathic Highland terrane (FHT), Pure Anorthosite terrane
(PAN), and non-PKT maria (nPKT). (b) The model thermal neutron counts, i.e.,
the difference between the Sn-covered and Cd-covered detectors. (c) The
measured thermal detector counts. The grey contours in (b) and (c) show the
regions defined in (a).

The non-uniformity of the Moon’s surface composition was included in the
modeling of neutron transport by modifying the equations developed in Feldman
_et al._ [17] such that the flux at the detectors at an altitude $R$, latitude
$\lambda_{R}$, and longitude $\varphi_{R}$ is related to the surface flux at
latitude $\lambda$, and longitude $\varphi$ by

$\Phi_{R}\left[\lambda_{R},\varphi_{R}\right](K_{R},\mu_{R},\phi_{R})=\left(\frac{K_{R}}{K}\right)\Phi[\lambda,\varphi](K,\mu,\phi)e^{-\frac{\Delta
t_{R}}{\tau_{n}}},$

where $K$ is the kinetic energy of the neutron at emission,
$K_{R}=K-V(R-R_{\leftmoon})$ is the kinetic energy at $R$, with
$V=\frac{GMm}{R_{\leftmoon}}$, $M$ the mass of the Moon, $m$ the neutron mass,
$\mu$ the cosine of the angle of emission with respect to the local zenith
$\theta$, $\mu_{R}=\sqrt{1-(R_{\leftmoon}/R)^{2}(K/K_{R})(1-\mu^{2})}$, and
$\Delta t_{R}$ is the time for transit for a neutron travelling from the
surface to an altitude $R$ [17], which is given by the expression

$\Delta
t_{R}=\frac{R_{\leftmoon}(m/2V)^{\frac{1}{2}}}{2(1-K/V)^{\frac{3}{2}}}\left(2\mu\left(1-\frac{K}{V}\right)^{\frac{1}{2}}\left(\frac{K}{V}\right)^{\frac{1}{2}}\times\right.\\\
\left(1-\left|\frac{\tan{\theta}}{\tan{\theta_{R}}}\right|\right)+\sin^{-1}\left(\frac{B}{(A^{2}+B^{2})^{\frac{1}{2}}}\right)\\\
\left.+\sin^{-1}\left(\frac{1-2K_{R}/V_{R}}{(A^{2}+B^{2})^{\frac{1}{2}}}\right)\right)$

with

$A=\sqrt{4\left(\frac{K}{V}\right)\left(1-\frac{K}{R}\right)\mu^{2}},\quad
B=\left(\frac{2K}{V}-1\right)\\\ \text{for}\quad\frac{K}{V}<1.$

A similar expression exists for $K/V>1$, see Feldman  _et al._ [17] for more
details and a derivation of the above equation.
$\Phi[\lambda,\varphi](K,\mu,\phi)$ is the flux at the surface of the Moon,
which is found by looking up the appropriate MCNPX output flux for the model
composition at latitude $\lambda$, and longitude $\varphi$. The location of
emission of a neutron is uniquely determined by the location and velocity at
detection by Kepler’s laws and spherical trigonometry such that

$\lambda=\sin^{-1}\left(\sin(\lambda_{R})\cos(\delta)+\cos(\lambda_{R})\sin(\delta)cos(\pi-\phi_{R})\right)$

and

$\varphi=\varphi_{R}+\tan^{-1}\left(\frac{\sin(\pi-\phi_{R})\sin(\delta)\cos(\lambda_{R})}{\cos(\delta)-\sin(\lambda_{R})\sin(\lambda)}\right),$

where $\delta$ is the difference between the true anomaly at neutron detection
$\nu_{0}$ and emission $\nu_{1}$. The true anomaly can be calculated from the
semi-major axis

$a=-\frac{VR}{2(K_{V})}$

and ellipticity

$e=\sqrt{1+4\frac{K}{V}\left(\frac{K}{V}-1\right)(1-\mu^{2})}$

via

$\nu_{0}=\cos^{-1}\left(\frac{a(1-e^{2})-R}{eR}\right).$

The third step in the generation of model count rates was the creation of an
MCNPX [14] model of the NS detectors and the calculation of their response to
neutrons that are incident with different momenta. Due to the small size of
the spacecraft and the location of the detectors on the end of a 2.5 m boom
(Fig. 5), a full simulation of the spacecraft was not necessary. Instead the
small spacecraft-originating background was removed from the data as described
in section II. The NS detectors were modeled in MCNPX as cylinders 20 cm long
and 5.7 cm in diameter filled with 3He with a density of $1.2\times 10^{-3}$ g
cm-3. Two detectors were simulated, one covered in each of 0.63 mm of Cd and
Sn with respective densities of 8.67 g cm-3 and 7.29 g cm-3. These detectors
were illuminated with plane wave neutrons with energies spanning the range to
which the detectors are sensitive, between $10^{-9}$ MeV and $10^{-2}$ MeV,
and a range of incidence angles covering the sphere. To calculate the
effective area $A$ of the detectors the total number of neutrons absorbed in
the active region of the detectors was tallied, where the active region was
defined as between 4.975 cm from one end and 3.705 cm from the other, and
multiplied by the number of primary events per unit area of the simulated
neutron source. This estimate of the active region is based on recent
calibration measurements of similar 3He sensors being used for NASA’s Psyche
mission [24]. The result of the simulation of the Sn-covered tube following
illumination with $10^{-9}$ MeV neutrons is shown in Fig. 4(a). As LP rotates
approximately six times during a single observation period of the NS, the
effective detector response was taken to be the instantaneous detector
response, calculated in MCNPX, averaged over rotation about the spacecraft’s
rotation axis. This is shown in Fig. 4(b). A cartoon of the spacecraft and
detectors is shown in Fig 5.

Figure 4: (a) The result of the MCNPX-based efficiency calculation for a Sn-
covered 3He tube to $10^{-9}$ MeV neutrons. (b) An example of angle-averaged
detector response used in the count rate simulation, based on the
instantaneous response shown in (a). The coordinates are conventional
spherical polars where the spacecraft’s rotation axis coincides with the
y-axis and the NS detector’s long axis coincides with the z-axis.

The final step in the generation of model count rates is to combine the models
of surface neutron flux, neutron transport and detector effective area by
integrating over neutron detection angle and energy such that

$J(\lambda_{R},\varphi_{R})=\int_{0}^{\infty}dK_{R}\int_{-1}^{1}d\mu_{R}\int_{0}^{2\pi}d\phi_{R}\left(\frac{-v\cdot\hat{n}}{v_{n}}\right)\times\\\
A(K_{R}^{\prime},\mu_{R}^{\prime},\phi_{R}^{\prime})\Phi_{R}\left[\lambda_{R},\varphi_{R}\right](K_{R},\mu_{R},\phi_{R}),$

where $v$ is the spacecraft velocity, $v_{n}$ the neutron velocity, and
$\hat{n}$ the detector normal. The velocities are in the Moon’s reference
frame and primed quantities are those in the frame of the detector. The
integration over neutron flux at the detector includes the contribution of
both upward going neutrons coming directly from the lunar surface and downward
going neutrons on elliptical trajectories that are returning to the Moon’s
surface. This effectively integrates over neutrons originating from across the
lunar surface.

The results of these calculations for both the Cd-covered and Sn-covered 3He
detectors is shown in Fig. 1(c). Fig. 3(b) shows the modeled thermal neutron
counts, which is the difference between modeled Sn-covered and Cd-covered
detector counts. It is the thermal neutron counts that will be used to
determine $\tau_{n}$ due to their greater sensitivity to that parameter, which
results from the lower velocity of these neutrons yielding a greater time of
flight.

During LP’s mapping phase, the spacecraft was spin stabilized about an axis
normal to the plane of the ecliptic. However, in the early elliptical orbits
this axis varied and was typically around 60∘ to the elliptic. The
spacecraft’s orientation was recorded before and after each maneuver and this
measured orientation was included in the modeling.

Figure 5:  A cartoon of the LP spacecraft and NS looking along the spacecraft
rotation axis. The locations of the Sn-covered and Cd-covered tubes are shown
in red and blue, respectively. The axis of rotation, pointing out of the page,
is indicated with the circled dot.

As the early elliptical orbits used in this study were not part of the science
phase of the LP mission they have been little used and not extensively
verified. We observed that there was a temporal offset between certain
spacecraft ephemerides (importantly altitude and velocity) and the detected
counts in the NS and the Gamma-Ray Spectrometer (GRS). That is, the times of
closest approach were not coincident with maxima in the rate of gamma ray
detection. To find the size of this temporal offset we tried shifting the
spacecraft ephemerides in time by increments of the NS integration period of
32 s. A $\chi^{2}$ comparison of solid-angle subtended by the Moon, following
various time offsets, with GRS counts integrated over all energies is shown in
Fig. 6, where the solid-angle has been scaled linearly to fit the GRS counts.
The magnitude of this offset was found to be $900.0\pm 0.5$ s, which was
estimated by fitting a quadratic to the points shown in Fig. 6. A comparison
with the GRS data is preferable to one with the NS measurements due to the
Doppler effects in the neutron data resulting from the particles’ velocities
being comparable with that of the spacecraft. We accommodated the offset in
ephemerides by translating the altitude and velocity parameters in time with
respect to the NS-measured counts. It will be shown in Section IV that the
final agreement between the data and model is good, which confirms the success
of this correction. Additionally, the spacecraft ephemerides for the first
half of the first elliptical orbit were missing so these data are not included
in the modeling.

Figure 6: $\chi^{2}$ statistic based on a comparison between the GRS counts
and solid-angle subtended by the Moon, with a linear scaling applied to minize
the misfit, for different temporal offsets applied to the spacecraft
ephemerides. The number of degrees of freedom $\nu$ is shown in the figure.

Similar to Wilson  _et al._ [6], the basis of the measurement technique is a
comparison of the NS-measured data with models of the constructed count rate
assuming different values of $\tau_{n}$. However, unlike that earlier study,
here we will use only the change in detected neutron counts with altitude to
measure the $\tau_{n}$, not the absolute count rate. In this approach each
model is normalized directly to the data with a single multiplicative factor
chosen to minimize the misfit between the model and data. This change
mitigates the largest source of systematic uncertainty in the previous
measurement, which was associated with errors in the normalization. Direct
normalization of the models to the data is made practical by the larger amount
of data taken by LP during its elliptical orbits compared with MESSENGER’s
flyby of Venus.

##  IV Results

We used the thermal neutron data (i.e., the difference between the Sn- and Cd-
covered NS detectors) to make the measurement of $\tau_{n}$, due to the larger
changes in the thermal neutron count rate compared with those in either the
Cd-covered or Sn-covered detectors, for given changes in $\tau_{n}$. Focussing
on the region of measurement-space that is most sensitive to the signal of
interest is important when attempting to measure $\tau_{n}$ at the Moon, due
to its relatively low mass and consequently small binding energy (0.029 eV
compared with 0.7 eV at Venus and a typical thermal neutron velocity of 0.1
eV) leading to a less steep change in the thermal neutron flux with altitude
around the Moon than at Venus [25].

A $\chi^{2}$ comparison of modeled thermal counts, derived assuming various
lifetimes, with the measured data using the uncertainties due to counting
statistics is shown in Fig. 7(a). This comparison includes all of the
observations shown in Fig. 1, including those taken after the initial nine
highly elliptical orbits when the spacecraft was in a near circular orbit.
Although qualitative agreement between the model and data can be seen in Figs.
1 & 3 it is clear from the value of $\chi^{2}$ shown in Fig. 7(a) and the
distribution of residuals in Fig. 7(b) that even the best-fitting model does
not adequately describe the data. Therefore, the value of $\tau_{n}$ obtained
from the analysis in Fig. 7(a) is subject to significant systematic errors.

Figure 7: (a) $\chi^{2}$ comparison of the thermal count model and data with
differing $\tau_{n}$ values using all of the observations shown in Fig. 1. The
number of degrees of freedom $\nu$ is shown in the figure. (b) Histogram of
the residuals of the best-fit model normalized by the uncertainty due to
counting statistics. (c) and (d) are similar to (a) and (b), respectively, but
make use only of the observations taken at longitudes greater than 180∘during
the nine initial highly elliptical orbits.

To avoid some of these systematic uncertainties we confine our analysis to
those observations taken where LP’s longitude was greater than 180∘ and before
the near-circular mapping orbit was reached (the grey regions in Fig 1 show
the extent of this restricted data set). With the first of these constraints
we removed the eastern lunar maria from the study as the residuals in this
region were anomalously large, which suggests that the complex composition of
this region has not been fully incorporated into the model. The effect of the
second constraint was to increase the average altitude of observation, which
reduces the need to precisely and accurately describe the lunar surface on
small scales in the model. The results of these data cuts on the analysis are
shown in Fig. 7(c,d). Although the statistical uncertainty in the result was
increased by the reduction in data volume, the new residuals are close to the
expected standard normal distribution. However, the minimum value of
$\chi^{2}$ is still not entirely consistent with a completely well-fitting
model given the expected statistical uncertainties on the measured data. This
implies that there exist additional small systematic uncertainties present in
the analysis.

We expect the largest source of systematic error to result from uncertainties
in the compositions used to generate the modeled counts. Based on the
elemental covariances reported in Prettyman  _et al._ [21] we explore the
magnitude of this systematic by modifying the compositions used for each of
the regions described in section III. Specifically, we perturbed the Fe
abundance by one standard deviation, in both the positive and negative
directions, and changed the abundances for each of the other elements by the
amount expected given these changes in Fe concentrations and the reported
covariances. Fe was chosen as variations in that element are the dominant
cause of changes in $\Sigma_{a}$ across the lunar surface [20]. We consider
this to represent a ‘$\pm$1$\sigma$’ perturbation of the reference
compositions. The elemental abundances of these perturbed compositions are
given in Table 2. As expected both of these models provide a less good fit to
the data with the minimum value of $\chi^{2}$ being greater than that
calculated for the reference composition (Fig. 8). A comparison of the
$+1\sigma$ models with the restricted data reveals that the best-fit lifetime
is approximately 7 s greater than the reference model, for the negative
perturbation the best-fit lifetime is around 3 s less than the reference case.
This analysis provides an estimate of the systematic uncertainty in our
analysis that results from uncertainties in the lunar composition that is less
than 10 s. Additional systematic uncertainties related to the definition of
the regions used in the model exist. Exploration of such systematics is a
critical outstanding task for understanding whether a spaced-based measurement
of $\tau_{n}$ could ever be competitive with laboratory measurements. However,
given the physical motivation for our selections of these regions and the
large parameter space from which they might be chosen, developing this task is
left for a detailed study of the potential systematics on a dedicated mission
to measure $\tau_{n}$.

Figure 8: $\chi^{2}$ comparison of the modelled restricted thermal neutron
count and data with differing $\tau_{n}$ values. The red curve shows a
comparison with the reference model as in Fig. 7c. The grey curves are based
on models using perturbations from the reference composition. The black curve
compares the data with a model including a latitude-dependent perturbation to
the count rates intended to simulate the effect of temperature variation. The
dotted lines indicate the locations of the minima of each of the curves.

An additional source of systematic error is the effect of variation in lunar
temperature with time of day and latitude, which was not included in the
modeling but will affect the thermal neutron flux above the Moon. Substantial
diurnal variation in thermal neutron flux due to temperature changes is not
expected as the diurnal temperature variation at a depth of 20 cm, which is
approximately the mean depth of origin of neutrons, is less than 30 K [26,
27]. However, latitudinal variations in temperature are more significant. The
mean surface temperature on the Moon, which is similar to the temperature
beneath the diurnal thermal wave envelope, varies between approximately 70 K
at the poles and 210 K at the equator [28] as $\cos^{\frac{1}{2}}\lambda$.
Lawrence  _et al._ [29] showed that similar variations in temperature gives
rise to a 2% variation in thermal neutron count rate when considering LP above
a surface of uniform composition and temperature. We estimated the magnitude
of the effect on $\tau_{n}$ of latitudinal variation in temperature by
superimposing a 2% equator-to-pole variation, of $\cos^{\frac{1}{2}}\lambda$
functional form, on the modeled thermal neutron count rate. The effect of this
modification on $\chi^{2}$ is illustrated in Figure 8. Fitting a parabola to
the reference and surface-temperature adjusted curves in the figure implies a
difference in $\tau_{n}$ of 0.2 s, which is negligible compared with both the
statistical errors and the composition-derived systematic described above. The
increase in $\chi^{2}$ over the reference model suggests that this simplistic
inclusion of the effect of temperature variation fails to accurately capture
the true variation in count rates due to surface temperature. For a precise
measurement of $\tau_{n}$ using data from a dedicated mission these diurnal
and latitudinal temperature variations could be included in the model in a
similar manner to the current implementation of variable composition. The
remaining systematic uncertainty would then be the result of the limits of our
knowledge of the lunar subsurface temperature.

As noted in Section III a temporal correction to the ephemerides was required.
The 1-$\sigma$ uncertainty on this correction, determined by the temporal
offset with $\chi^{2}$ equal to one greater than the minimum value, is 0.5 s,
much less than one observation period of 32 s as shown in Fig. 6. This
ephemerides offset is unusually large and its presence is a consequence of the
data not being part of the science phase of the mission, so not previously
being well validated. We would therefore not expect similar ephemerides-
derived systematics to be present on any future mission. Ordinarily,
uncertainties in spacecraft position are on the order of meters [30] rather
than the hundreds of meters travelled by LP during a 0.5-s interval.

For a space-based method similar to that described here to become competitive
with current and near-future laboratory-based measurements, the above sources
of systematic uncertainty along with potentially significant contributions due
to uncertainties in the instrument response function and the Monte Carlo
modeling, including nuclear cross-section uncertainties, will need to be
corrected for or mitigated [30].

We find $\tau_{n}=887\pm 14_{\text{stat}}{\>{}^{+7}_{-3\>\text{syst}}}$ s,
which is a 0.5$\sigma$ difference from the Particle Data Group (PDG) value of
$879\pm 0.4$ s [31]. This result is shown in comparison with laboratory-based
results and the previous space based result from Wilson  _et al._ [6] in Fig.
9. The measurement presented here further demonstrates the feasibility of
measuring $\tau_{n}$ using a space-based experiment. Additionally, the
extension of the modeling to include non-uniform surface compositions extends
the range of bodies where the space-based $\tau_{n}$ measurement might be
performed and reduces what was a major systematic uncertainty in the earlier
space-based $\tau_{n}$ measurement. The combined space-based measurement with
standard error (shown as a grey region in Fig. 9) of
$\tau_{n}^{\textrm{space}}=883\pm 17$ s is in good agreement with the PDG
value.

Figure 9:  Particle Data Group and other recent measurements of $\tau_{n}$
[31, 32, 33, 34], along with those of this paper and Wilson  _et al._ [6]. The
shaded regions represent the standard error on the uncertainty-weighted mean
lifetime in each class of measurement.

##  V Conclusions

Using data taken by the neutron detectors on NASA’s Lunar Prospector
spacecraft during the elliptical orbits completed immediately following the
spacecraft’s capture into a lunar orbit we found the neutron lifetime to be
$\tau_{n}=887\pm 14_{\text{stat}}{\>{}^{+7}_{-3\>\text{syst}}}$ s. In
combination with the previous measurement of $\tau_{n}$ from space using
MESSENGER data [6], this result firmly establishes the feasibility of making a
measurement of $\tau_{n}$ from space. The statistical uncertainties are large
compared with laboratory measurements due to the relatively short duration of
the elliptical phase of the LP mission. Consequently, small amounts of data
are taken, which is a result of the LP mission not being designed for this
analysis. However, the modeling developed in this study, which enables non-
uniform surfaces to be simulated, provides a key component to carrying out
more thorough measurements of $\tau_{n}$ around airless bodies. Since the
space-based method of constraining $\tau_{n}$ has entirely separate systematic
uncertainties to the two existing classes of laboratory experiments, future
space-based measurements with higher statistical precision and well-
constrained systematics may provide a route to make progress beyond the
current disagreement between the bottle and beam results. Achieving this goal
would require a reduction of the systematics to the 1 s level, which
necessitates further development of the space-based method that builds on the
result of this paper.

###### Acknowledgements.

We wish to thank the anonymous reviewer whose comments spurred substantial
improvement in this paper. We gratefully acknowledge the support of the U.S.
Department of Energy Office of Nuclear Physics (Grant No. DE-SC0019343). VRE
was supported by the Science and Technology Facilities Council (STFC) grant
ST/P000541/1. JAK acknowledges support from STFC grants ST/N001494/1 and
ST/T002565/1.

##  Appendix A Model Compositions

Table 1: The elemental mass fractions of the major crustal terranes used in the reference model: the Procellarum KREEP terrane (PKT), non-PKT maria (nPKT), South Pole Aitken terrane (SPAT), Feldspathic highland terrane (FHT), and pure anorthosite terrane (PAN) [18, 19]. Element | FHT | PAN | SPA | nPKT | PKT
---|---|---|---|---|---
H | 5.13$\times 10^{-5}$ | 2.86$\times 10^{-5}$ | 5.37$\times 10^{-5}$ | 4.35$\times 10^{-5}$ | 4.91$\times 10^{-5}$
O | 0.455 | 0.463 | 0.444 | 0.418 | 0.403
Na | 0.00228 | 0.00236 | 0.00225 | 0.00235 | 0.00227
Mg | 0.0368 | 0.0271 | 0.0489 | 0.0578 | 0.0562
Al | 0.135 | 0.146 | 0.114 | 0.0695 | 0.075
Si | 0.214 | 0.221 | 0.209 | 0.198 | 0.186
K | 0.000441 | 0.000306 | 0.00101 | 0.0012 | 0.00254
Ca | 0.12 | 0.125 | 0.115 | 0.104 | 0.114
Ti | 0.00278 | 0.00216 | 0.00539 | 0.0208 | 0.0204
Fe | 0.0344 | 0.0125 | 0.0596 | 0.128 | 0.141
Sm | 8.43$\times 10^{-7}$ | 0 | 3.32$\times 10^{-6}$ | 1.8$\times 10^{-6}$ | 7.86$\times 10^{-6}$
Gd | 9.87$\times 10^{-7}$ | 0 | 3.89$\times 10^{-6}$ | 2.11$\times 10^{-6}$ | 9.19$\times 10^{-6}$
Th | 8.69$\times 10^{-8}$ | 8.83$\times 10^{-8}$ | 2.23$\times 10^{-6}$ | 2.72$\times 10^{-6}$ | 6.61$\times 10^{-6}$
U | 2.36$\times 10^{-8}$ | 2.40$\times 10^{-8}$ | 6.06$\times 10^{-7}$ | 7.41$\times 10^{-7}$ | 1.80$\times 10^{-6}$
Table 2: The elemental mass fractions of the major crustal terranes used in
the perturbed, ‘1$\sigma$’, models: the Procellarum KREEP terrane (PKT), non-
PKT maria (nPKT), South Pole Aitken Terrane (SPAT), Feldspathic Highland
Terrane (FHT), and Pure Anorthosite terrane (PAN) [18, 19]. $\bf{+1\sigma}$
---
Element | FHT | PAN | SPA | nPKT | PKT
H | 5.21$\times 10^{-5}$ | 2.92$\times 10^{-5}$ | 5.44$\times 10^{-5}$ | 4.39$\times 10^{-5}$ | 4.93$\times 10^{-5}$
O | 0.46 | 0.466 | 0.452 | 0.43 | 0.417
Na | 0.00232 | 0.00241 | 0.00228 | 0.00237 | 0.00229
Mg | 0.0373 | 0.0277 | 0.0495 | 0.0584 | 0.0566
Al | 0.122 | 0.132 | 0.101 | 0.0551 | 0.0613
Si | 0.218 | 0.228 | 0.212 | 0.199 | 0.186
K | 0.000454 | 0.00032 | 0.00103 | 0.00122 | 0.00257
Ca | 0.122 | 0.128 | 0.117 | 0.105 | 0.114
Ti | 0.0031 | 0.00241 | 0.00604 | 0.0222 | 0.022
Fe | 0.0348 | 0.013 | 0.0592 | 0.127 | 0.139
Sm | 8.56$\times 10^{-7}$ | 0 | 3.37$\times 10^{-6}$ | 1.82$\times 10^{-6}$ | 7.9$\times 10^{-6}$
Gd | 1.00$\times 10^{-6}$ | 0 | 3.94$\times 10^{-6}$ | 2.13$\times 10^{-6}$ | 9.24$\times 10^{-6}$
Th | 8.95$\times 10^{-8}$ | 9.13$\times 10^{-8}$ | 2.26$\times 10^{-6}$ | 2.74$\times 10^{-6}$ | 6.64$\times 10^{-6}$
U | 2.78$\times 10^{-8}$ | 2.83$\times 10^{-8}$ | 6.51$\times 10^{-7}$ | 7.88$\times 10^{-7}$ | 1.85$\times 10^{-6}$
$\bf{-1\sigma}$
Element | FHT | PAN | SPA | nPKT | PKT
H | 5.06$\times 10^{-5}$ | 2.8$\times 10^{-5}$ | 5.31$\times 10^{-5}$ | 4.31$\times 10^{-5}$ | 4.88$\times 10^{-5}$
O | 0.45 | 0.46 | 0.437 | 0.406 | 0.39
Na | 0.00225 | 0.0023 | 0.00222 | 0.00233 | 0.00226
Mg | 0.0362 | 0.0265 | 0.0483 | 0.0573 | 0.0559
Al | 0.147 | 0.16 | 0.126 | 0.0836 | 0.0885
Si | 0.209 | 0.214 | 0.206 | 0.198 | 0.186
K | 0.000427 | 0.000293 | 0.000988 | 0.00118 | 0.00252
Ca | 0.118 | 0.122 | 0.114 | 0.103 | 0.113
Ti | 0.00246 | 0.00192 | 0.00474 | 0.0194 | 0.0188
Fe | 0.034 | 0.0121 | 0.0599 | 0.13 | 0.143
Sm | 8.31$\times 10^{-7}$ | 0 | 3.28$\times 10^{-6}$ | 1.79$\times 10^{-6}$ | 7.81$\times 10^{-6}$
Gd | 9.72$\times 10^{-7}$ | 0 | 3.84$\times 10^{-6}$ | 2.09$\times 10^{-6}$ | 9.14$\times 10^{-6}$
Th | 8.43$\times 10^{-8}$ | 8.53$\times 10^{-8}$ | 2.19$\times 10^{-6}$ | 2.7$\times 10^{-6}$ | 6.58$\times 10^{-6}$
U | 1.96$\times 10^{-8}$ | 1.98$\times 10^{-8}$ | 5.63$\times 10^{-7}$ | 6.95$\times 10^{-7}$ | 1.76$\times 10^{-6}$

## References

  * Wietfeldt and Greene [2011] F. E. Wietfeldt and G. L. Greene, Rev. Mod. Phys. 83, 1173 (2011).
  * Dubbers and Schmidt [2011] D. Dubbers and M. G. Schmidt, Rev. Mod. Phys. 83, 1111 (2011).
  * Cyburt  _et al._ [2016] R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, Reviews of Modern Physics 88, 015004 (2016), arXiv:1505.01076 [astro-ph.CO] .
  * Fornal and Grinstein [2018] B. Fornal and B. Grinstein, Phys. Rev. Lett. 120, 191801 (2018).
  * Berezhiani [2019] Z. Berezhiani, European Physical Journal C 79, 484 (2019), arXiv:1807.07906 [hep-ph] .
  * Wilson  _et al._ [2020] J. T. Wilson, D. J. Lawrence, P. N. Peplowski, V. R. Eke, and J. A. Kegerreis, Physical Review Research 2, 023316 (2020), arXiv:2006.10008 [nucl-ex] .
  * Lingenfelter  _et al._ [1961] R. E. Lingenfelter, E. H. Canfield, and W. N. Hess, Journal of Geophysical Research (1896-1977) 66, 2665 (1961).
  * Feldman  _et al._ [2004] W. C. Feldman, K. Ahola, B. L. Barraclough, R. D. Belian, R. K. Black, R. C. Elphic, D. T. Everett, K. R. Fuller, J. Kroesche, D. J. Lawrence, S. L. Lawson, J. L. Longmire, S. Maurice, M. C. Miller, T. H. Prettyman, S. A. Storms, and G. W. Thornton, Journal of Geophysical Research (Planets) 109, E07S06 (2004).
  * Boynton  _et al._ [2004] W. Boynton, W. Feldman, I. Mitrofanov, L. Evans, R. Reedy, S. Squyres, R. Starr, J. Trombka, C. d’Uston, J. Arnold, P. Englert, A. Metzger, H. Wänke, J. Brückner, D. Drake, C. Shinohara, C. Fellows, D. Hamara, K. Harshman, K. Kerry, C. Turner, M. Ward, H. Barthe, K. Fuller, S. Storms, G. Thornton, J. Longmire, M. Litvak, and A. Ton’chev, Space Science Reviews 110, 37 (2004).
  * Goldsten  _et al._ [2007] J. O. Goldsten, E. A. Rhodes, W. V. Boynton, W. C. Feldman, D. J. Lawrence, J. I. Trombka, D. M. Smith, L. G. Evans, J. White, N. W. Madden, P. C. Berg, G. A. Murphy, R. S. Gurnee, K. Strohbehn, B. D. Williams, E. D. Schaefer, C. A. Monaco, C. P. Cork, J. Del Eckels, W. O. Miller, M. T. Burks, L. B. Hagler, S. J. Deteresa, and M. C. Witte, Space Sci. Rev. 131, 339 (2007).
  * Prettyman  _et al._ [2011] T. H. Prettyman, W. C. Feldman, H. Y. McSween, R. D. Dingler, D. C. Enemark, D. E. Patrick, S. A. Storms, J. S. Hendricks, J. P. Morgenthaler, K. M. Pitman, and R. C. Reedy, Space Science Reviews 163, 371 (2011).
  * Maurice  _et al._ [2004] S. Maurice, D. J. Lawrence, W. C. Feldman, R. C. Elphic, and O. Gasnault, Journal of Geophysical Research (Planets) 109, E07S04 (2004).
  * Feldman  _et al._ [1999] W. Feldman, B. Barraclough, K. Fuller, D. Lawrence, S. Maurice, M. Miller, T. Prettyman, and A. Binder, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 422, 562 (1999).
  * Pelowitz [2005] D. B. Pelowitz,  _MCNPX User’s Manual, version 2.5.0._ , Los Alamos National Laboratory, Los Alamos, N. Mex. (2005).
  * Masarik and Reedy [1996] J. Masarik and R. C. Reedy, J. Geophys. Res. 101, 18891 (1996).
  * McKinney  _et al._ [2006] G. W. McKinney, D. J. Lawrence, T. H. Prettyman, R. C. Elphic, W. C. Feldman, and J. J. Hagerty, Journal of Geophysical Research (Planets) 111, E06004 (2006).
  * Feldman  _et al._ [1989] W. C. Feldman, D. M. Drake, R. D. O’dell, F. W. Brinkley, and R. C. Anderson, J. Geophys. Res. 94, 513 (1989).
  * Jolliff  _et al._ [2000] B. L. Jolliff, J. J. Gillis, L. A. Haskin, R. L. Korotev, and M. A. Wieczorek, J. Geophys. Res. 105, 4197 (2000).
  * Peplowski  _et al._ [2016] P. N. Peplowski, A. W. Beck, and D. J. Lawrence, Journal of Geophysical Research (Planets) 121, 388 (2016).
  * Elphic  _et al._ [2000] R. C. Elphic, D. J. Lawrence, W. C. Feldman, B. L. Barraclough, S. Maurice, A. B. Binder, and P. G. Lucey, J. Geophys. Res. 105, 20333 (2000).
  * Prettyman  _et al._ [2006] T. H. Prettyman, J. J. Hagerty, R. C. Elphic, W. C. Feldman, D. J. Lawrence, G. W. McKinney, and D. T. Vaniman, Journal of Geophysical Research (Planets) 111, E12007 (2006).
  * Peplowski  _et al._ [2015] P. N. Peplowski, D. J. Lawrence, W. C. Feldman, J. O. Goldsten, D. Bazell, L. G. Evans, J. W. Head, L. R. Nittler, S. C. Solomon, and S. Z. Weider, Icarus 253, 346 (2015).
  * Lawrence  _et al._ [2015] D. J. Lawrence, P. N. Peplowski, J. B. Plescia, B. T. Greenhagen, S. Maurice, and T. H. Prettyman, Icarus 255, 127 (2015).
  * Peplowski  _et al._ [2020] P. N. Peplowski, Z. W. Yokley, M. Liebel, S. Cheng, R. C. Elphic, S. F. Hoogerheide, D. J. Lawrence, and J. S. Nico, Nuclear Instruments and Methods in Physics Research A 982, 164574 (2020).
  * Hess  _et al._ [1961] W. N. Hess, E. H. Canfield, and R. E. Lingenfelter, J. Geophys. Res. 66, 665 (1961).
  * Vasavada  _et al._ [1999] A. R. Vasavada, D. A. Paige, and S. E. Wood, Icarus 141, 179 (1999).
  * Little  _et al._ [2003] R. C. Little, W. C. Feldman, S. Maurice, I. Genetay, D. J. Lawrence, S. L. Lawson, O. Gasnault, B. L. Barraclough, R. C. Elphic, T. H. Prettyman, and A. B. Binder, Journal of Geophysical Research (Planets) 108, 5046 (2003).
  * Williams  _et al._ [2017] J.-P. Williams, D. Paige, B. Greenhagen, and E. Sefton-Nash, Icarus 283, 300 (2017).
  * Lawrence  _et al._ [2006] D. J. Lawrence, W. C. Feldman, R. C. Elphic, J. J. Hagerty, S. Maurice, G. W. McKinney, and T. H. Prettyman, Journal of Geophysical Research (Planets) 111, E08001 (2006).
  * Lawrence  _et al._ [2021] D. J. Lawrence, J. T. Wilson, and P. N. Peplowski, Nuclear Instruments and Methods in Physics Research A 988, 164919 (2021), arXiv:2011.06095 [nucl-ex] .
  * Zyla [2020] A. Zyla, P. (Particle Data Group), Prog. Thror. Exp. Phys. 2020, 083C01 (2020).
  * Yue  _et al._ [2013] A. T. Yue, M. S. Dewey, D. M. Gilliam, G. L. Greene, A. B. Laptev, J. S. Nico, W. M. Snow, and F. E. Wietfeldt, Phys. Rev. Lett. 111, 222501 (2013), arXiv:1309.2623 [nucl-ex] .
  * Byrne  _et al._ [1996] J. Byrne, P. G. Dawber, C. G. Habeck, S. J. Smidt, J. A. Spain, and A. P. Williams, EPL (Europhysics Letters) 33, 187 (1996).
  * Spivak [1988] P. E. Spivak, Soviet Physics - JETP (English Translation) 97, 1735 (1988).