# Lensing of Gravitational Waves as a Novel Probe of Graviton Mass

Adrian Ka-Wai Chung<EMAIL_ADDRESS>Theoretical Particle Physics and
Cosmology Group, Department of Physics, King’s College London, University of
London, Strand, London, WC2R 2LS, United Kingdom  Tjonnie G.F. Li
<EMAIL_ADDRESS>Department of Physics, The Chinese University of Hong Kong,
Shatin, N.T., Hong Kong  Institute for Theoretical Physics, KU Leuven,
Celestijnenlaan 200D, B-3001 Leuven, Belgium  Department of Electrical
Engineering (ESAT), KU Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium

###### Abstract

The diffraction patterns of lensed gravitational waves encode information
about their propagation speeds. If gravitons have mass, the dispersion
relation and speed of gravitational waves will be affected in a frequency-
dependent manner, which would leave traces in the diffraction pattern if the
waves are lensed. In this paper, we study how the alternative dispersion
relation induced by massive gravitons affects gravitational waves lensed by
point-mass lenses, such as intermediate-mass black holes. We find that the
waveform morphology of lensed dispersive gravitational waves depends on the
graviton mass more sensitively than their unlensed counterpart. Together with
lensing amplification, the waveform-morphology modifications due to lensing
can improve the measurement accuracy of the graviton mass. A single lensed
gravitational-wave signal enables us to measure the graviton mass with an
accuracy comparable with the combined measurement across $\mathcal{O}(10^{3})$
unlensed signals. Our method allows us to incorporate lensed gravitational-
wave signals into existing graviton-mass measurements. Our method can also be
generalized to other lens types, gravitational-wave sources, and detector
networks, preparing ourselves for thoroughly understanding the nature of
gravitational waves in the era of lensed gravitational-wave astronomy.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

When gravitational waves (GWs) propagate near a massive compact object, the
propagation direction will be changed due to the gravity of the compact object
Lensing_GW_01 ; Lensing_GW_02 ; Lensing_GW_03 ; Lensing_IMBH_02 ;
Takahashi_2003 ; Takahashi_2017 . This phenomenon is known as lensing. The
lensed waves interfere among themselves to form new wave patterns. The lensing
pattern depends on the nature of the lens, lensing geometry, and interference
between different (possible) rays. Lensing is a crucial astrophysical probe.
For example, lensing of GWs might provide us with the information of the
existence of intermediate-mass black holes (IMBHs) (of mass $\sim
10^{2}M_{\odot}-10^{3}M_{\odot}$) Lensing_IMBH_01 , cosmological expansion
Hubble_constant ; Cremonese2021  and testing general relativity Suvodip01 ;
Suvodip02 ; Collett_Bacon_2017 ; Fan_2017 ; Baker . Since the direct detection
of GWs, searches for lensed GWs have been popular Lensing_O1_O2 ;
O3a_lensing_paper ; Lensing_LIGO_01 . Nonetheless, no conclusive evidence of
lensed GW signals has been found as of the time of writing O3a_lensing_paper .

Another popular aspect of studies is the dispersion relation of GWs LIGO_01 ;
LIGO_02 ; LIGO_03 ; LIGO_04 ; LIGO_05 ; LIGO_06 ; LIGO_07 ; LIGO_08 ; LIGO_09
; LIGO_10 ; LIGO_11 ; LIGO_12 ; LIGO_13 ; LIGO_14 ; LIGO_15 . According to
general relativity, in vacuum, GWs propagate at the speed of light and obey
the dispersion relation of $\omega=ck$, where $\omega$ is the angular
frequency, $c$ is the speed of light, and $k$ is the wave number. This
dispersion relation implies that gravitons are massless. On the other hand,
some alternative theories explore the possibility of massive gravitons, see,
e.g., MG_review_01  for a review. If the gravitons have mass, GWs obey an
alternative dispersion relation, leading to different propagation speeds for
various frequencies. As GWs propagate, dephasing will be developed across
different frequencies. By (not) measuring the dephasing, we can constrain the
graviton mass m_g_bound_review_01 ; m_g_bound_review_02 ; Will_2012 .
Alternatively, the near-field behavior of black holes has also been suggested
as a probe of the graviton mass Chung_2019 . As of the time of writing, no
evidence of massive gravitons has been found using this method LIGO_03 ;
LIGO_06 ; LIGO_07 ; LIGO_11 .

These methods are capable of extracting the graviton mass from unlensed GW
signals. On the other hand, lensed GW signals are expected to be detected in
the future Ken_2016 ; Lensing_LIGO_01 . To thoroughly understand the nature of
GWs in different astrophysical scenarios, developing a test of the dispersion
of lensed GWs becomes increasingly pressing. Moreover, since the amplitude of
lensed GWs shows more variations across the frequency than the amplitude and
phase of unlensed waves Takahashi_2003 ; lensingGW , the waveform morphology
of lensed dispersive GWs may depend on the graviton mass more sensitively than
the unlensed GWs. Besides, the amplification introduced by lensing may
contribute to an improved measurement accuracy of the graviton mass compared
to the unlensed case. Furthermore, the dispersion relation of GWs
corresponding to the massive graviton also changes the time delay of waves of
different frequencies in different directions, leading to additional features
of the resultant lensing pattern. These considerations prompt us to explore
measuring the graviton mass from lensed GW signals. Measuring the graviton
mass by lensing also makes relevant tests more complete in at least two ways.
First, lensing involves the strength of gravity intermediately between the
near and far fields. Second, our work enables us to incorporate lensed GW
signals into the measurement of the graviton mass, which has thus far focused
only on unlensed signals, better preparing ourselves for the era of lensed GW
astronomy.

Throughout the paper, $m_{g}$ is in the unit of $c=\hbar=1$ (so $h=2\pi$),
while the mass of compact objects (such as black hole and lens) are in the
unit of $c=G=1$.

##  II Lensing Pattern of Gravitational Waves with Dispersion

###  II.1 Assumptions and approximations

This work makes a few assumptions:

  1. A.1

Following Will_2012 ; Ajith_2010 ; Yagi_2015 , we assume perfect screening of
gravity due to the mass of graviton Yunes_2018 . In other words, general-
relativistic limits are recovered at a length scale shorter than the Compton
wavelength of the graviton. This assumption implies that we will ignore the
effects on the dynamics of binary black hole mergers due to the graviton mass.
In the context of lensing, this assumption implies that at a sufficiently far
distance $r$ the Newtonian gravitational potential due to a black hole (point-
mass lens) of mass $M$ is given by $-\frac{M}{r}$.

  2. A.2

We focus only on the effects on GW lensing due to the graviton mass. Other
consequences of lensing, such as modifications on polarization
Lensing_polarization  and phase shift Lensing_phase_shift_01 ;
Lensing_birefringence , will be omitted. These are acceptable approximations
because including these effects will include more contrasting features to
graviton-mass measurement. In this work, we focus on improvements on the
constraint or measurement accuracy of the graviton mass by considering the
combined effects of lensing and dephasing induced by the graviton mass.

###  II.2 Method

If gravitons have mass, phenomenologically, the dispersion relation of GWs
will be altered to Will_2012

$\omega^{2}=k^{2}+m_{g}^{2},$ (1)

where $m_{g}$ is the mass of graviton 111Alternatively, this equation can be
interpreted as a definition of the massive graviton which leads to the
dispersion of gravitational perturbations. In this work, we refer “the mass of
graviton” to $m_{g}$ defined by Eq. 1. If $m_{g}\ll k$, the propagation speed
of dispersive GWs that obey this dispersion relation can be approximated by
the following equation:

$v_{g}(f)\approx 1-\frac{1}{8\pi^{2}}\frac{m_{g}^{2}}{f^{2}}.$ (2)

When propagating in a flat space-time, the dispersive GWs obeying Eq. 1 will
acquire a dephasing due to the difference in propagation speeds among
different frequencies Will_2012 ,

$\begin{split}&\Psi_{\rm disp}(f;m_{g})=-\frac{\pi
D_{0}}{\lambda_{g}^{2}}\frac{1}{(1+z)f},\end{split}$ (3)

where $\lambda_{g}=1/m_{g}$ is Compton’s wavelength of the graviton, $D_{0}$
is the propagation distance from the source to the detector and $z$ is the
redshift of the source binary. Thus, in the frequency domain, the waveform of
unlensed dispersive GWs is

$\tilde{h}_{\text{\rm disp}}(f)=\tilde{h}(f)e^{i\Psi_{\text{disp}}(f)},$ (4)

where $\tilde{h}(f)$ is the original (unlensed) GR waveform (see, e.g.,
(Sathyaprakash_2009, ; IMRPhenom_2011, ; IMRPhenom_2014, ) for GR waveform
approximants).

When encountering a massive compact object, such as an intermediate-mass black
hole, GWs will be lensed. The lensing effect is characterized by the
amplification function (or transmission factor) (Schneider_01, ; Schneider_02,
), $F$, which is the ratio of lensed-wave amplitude to unlensed-wave
amplitude,

$\tilde{h}_{L}(f)=F(f)\tilde{h}(f),$ (5)

where $\tilde{h}_{L}(f)$ is the lensed waveform and $\tilde{h}(f)$ is the
unlensed waveform. Given a lensing geometry, $F(f)$ can be computed by
Takahashi_2003 ; Takahashi_2017 ; Nakamura_02

$\begin{split}F(f;\vec{\theta}_{s})=&\frac{D_{L}D_{S}}{D_{LS}}\xi_{0}^{2}\frac{(1+z_{L})}{i}\frac{f}{v_{g}}\\\
&\times\int d^{2}\vec{\theta}_{L}\exp\left[2\pi
ift_{d}(\vec{\theta}_{L},\vec{\theta}_{s})\right],\end{split}$ (6)

where $v_{g}$ is GW propagation speed; $D_{L},D_{S}$, and $D_{LS}$ are,
respectively, the lens-to-observer distance, the source-to-observer distance,
and the source-to-lens distance; $z_{L}$ is the redshift of lens;
$\vec{\theta}_{s}$ is the displacement from optical axis to the source on
source plane; $\vec{\theta}_{L}$ is the displacement from optical axis to lens
on lens plane; and $t_{d}$ is the time delay between the lensed ray and
unlensed ray,

$t_{d}(\vec{\theta},\vec{\theta}_{s})=\frac{(1+z_{L})}{v_{g}}\Bigg{[}\frac{D_{L}D_{S}}{2D_{LS}}|\vec{\theta}_{s}-\vec{\theta}|^{2}-\psi(\vec{\theta}_{s})\Bigg{]},$
(7)

where $\psi(\vec{\theta}_{\rm s})$ is the lensing potential. Overall, $t_{d}$
also depends on $v_{g}$, $\vec{\theta}_{s}$, and lens $\vec{\theta}_{L}$, and
$\xi_{0}$ is a length scale.

We note that the amplification function Eq. 6 depends on $\frac{f}{v_{g}}$ as
a whole. Thus, the amplification function of GWs of the massive graviton is
just that of GWs without dispersion with the following replacement:

$f\rightarrow\beta(f)f,$ (8)

where

$\beta(f)=\frac{c}{v_{g}(f)}\approx 1+\frac{1}{2}\frac{m_{g}^{2}}{f^{2}}.$ (9)

From Eq. 9, we expect that the modifications to the lensing pattern due to the
dispersion relation Eq. 1 are manifest for $m_{g}\geq 10^{-14}\rm eV$,
corresponding to the energy scale of $hf$ at $f=10\rm Hz$.

As a proof of principle, in this work we focus on the case of a point-mass
lens, such as a black hole. For a point-mass lens, the amplification function
can be analytically evaluated as Takahashi_2003 ; Takahashi_2017

$\begin{split}&F(f;M_{\rm len},y,m_{g})\\\
&=\exp\Big{(}\frac{\pi}{4}w\beta\Big{)}\Big{(}\frac{w}{2}\beta\Big{)}^{i\frac{w}{2}\beta}\Gamma\Big{(}1-i\frac{w}{2}\beta\Big{)}{}_{1}F_{1}\Big{(}i\frac{w}{2}\beta,1;i\frac{w}{2}\beta
y^{2}\Big{)},\end{split}$ (10)

where $M_{\rm len}$ is the redshifted mass of the lens, $y$ is the impact
parameter of lensing, $\Gamma$ is the (complex) Gamma function, ${}_{1}F_{1}$
is confluent hypergeometric function, and $w=8\pi M_{\rm len}f$ is the
dimensionless frequency. The resulting lensed waveform of GWs corresponding to
the massive gravitons can be written as

$\tilde{h}_{L}(f;m_{g})=F(f;M_{\rm len},y,m_{g})\tilde{h}(f)e^{i\Psi_{\rm
disp}(f)}.$ (11)

Note that, according to Eq. 2, GWs of different frequencies travel at
different speeds. The only constant achromatic speed is the speed of light.
Therefore, the effects described by Eq. 11 are not degenerate with a constant
change of propagation speed of GWs. Thus, the effects of the massive gravitons
can be distinguished upon gravitational-wave detection.

Figure 1: The lensing amplification function due to an intermediate mass black
hole of redshifted mass $M_{\rm len}$ of 400 $M_{\odot}$ at $y=0.9$
corresponding to $m_{g}=0$ (solid blue line), $m_{g}=10^{-14}\rm eV$ (dashed
red line), and $m_{g}=10^{-22}\rm eV$ (dashed dotted green). At
$m_{g}=10^{-22}\rm eV$, which is close to the existing constraints by
gravitational-wave detection, the change of the amplification due to $m_{g}$
is not very manifest as the solid blue line and dotted green line overlap
almost entirely.

Fig. 1 plots the $F(f)$ corresponding to the lensing by an intermediate-mass
black hole of redshifted lens mass of $400M_{\odot}$ and $y=0.9$ for $m_{g}=0$
(solid blue), $m_{g}=10^{-14}\rm eV$ (dashed red) and $m_{g}=10^{-22}\rm eV$
(dotted green) as a function of $f$. For $m_{g}=10^{-14}\rm eV$, we find that
the modifications of the amplification function is manifest for the low-
frequency regime ($f\leq 10^{2}\rm Hz$), in which $\beta$ changes
significantly with $m_{g}$. As GW frequency increases, the changes of the
amplification function due to the alternative dispersion become increasingly
less manifest because the ultrarelativistic limit $E\approx p$ has been
attained. For $m_{g}=10^{-22}\rm eV$, over $f\in[10,10^{3}]\rm Hz$, the
modifications due to the graviton mass are not visible, as expected because
$|\beta(f)-1|\sim 10^{-15}\ll 1$ for $m_{g}=10^{-22}\rm eV$.

Eq. 11 suggests that GW lensing may help to improve the measurement of $m_{g}$
in at least three ways.

  1. R.1

Lensing changes the waveform morphology of the signal. Specifically, because
of the modulation by the amplification function (Fig. 1), the amplitude and
phase of lensed GWs show more variations across the frequencies than the
unlensed GWs. This beating pattern may make the waveform morphology of lensed
dispersive GWs depend on the graviton mass more sensitively than unlensed
waves.

  2. R.2

Lensing increases the signal-to-noise ratio (SNR).

  3. R.3

The graviton mass modifies the amplification function, making the waveform
morphology of lensed dispersion GWs depend on the graviton mass even more
sensitively.

However, judging from Fig. 1, for $m_{g}$ close to the existing constraints of
the graviton mass ($\sim 10^{-22}\rm eV$ LIGO_07 ; LIGO_11 ; LIGO_15 ), the
changes of the amplification function due to $m_{g}$ are not significant.
Therefore, R.3 is unlikely to contribute to any significant improvement. In
what follows, we focus on investigating the roles of R.1 and R.2.

As a first step, we compare the similarity of the waveform of both lensed and
unlensed dispersive GWs of a given $m_{g}$ to the dispersive waves of other
$m_{g}$. In general, the similarity between two waveforms, $\tilde{h}_{1}(f)$
and $\tilde{h}_{2}(f)$, can be gauged by the match between $\tilde{h}_{1}$ and
$\tilde{h}_{2}$, defined as

$\mathcal{M}=\frac{\braket{\tilde{h}_{1}}{\tilde{h}_{2}}}{\sqrt{\braket{\tilde{h}_{1}}{\tilde{h}_{1}}\braket{\tilde{h}_{2}}{\tilde{h}_{2}}}},$
(12)

where the braket notation denotes the noise-weighted inner product
Inner_product ,

$\braket{\tilde{h}_{1}}{\tilde{h}_{2}}=4\text{Re}\int_{0}^{+\infty}df\frac{\tilde{h}_{1}(f)\tilde{h}_{2}^{\dagger}(f)}{S_{n}(f)},$
(13)

and $S_{n}(f)$ is the one-sided power-spectral density of the detector.
Throughout this work, we assume GW signals are detected by the Advanced LIGO
and Virgo detectors operating at their design sensitivity LIGO_noise_curve ;
Virgo_noise_curve . To investigate how sensitive lensed dispersive GWs depends
on $m_{g}$, we chose

$\begin{split}&\tilde{h}_{1}(f)=\tilde{h}_{L}(f;m_{g}=m_{g}^{\rm inj}),\\\
&\tilde{h}_{2}(f)=\tilde{h}_{L}(f;m_{g}),\end{split}$ (14)

where $m_{g}^{\rm inj}$ is a given value of $m_{g}$ and $\tilde{h}_{L}$ is the
lensed waveform defined by Eq. 11 and LABEL:eq:Amp_fn_mg. Using this waveform,
we have defined a match as a function of $m_{g}$ for lensed dispersive GWs.
Similarly, we can define a match for unlensed dispersive GWs by replacing
$\tilde{h}_{L}(f;m_{g}=m_{g}^{\rm inj})\rightarrow\tilde{h}(f;m_{g}=m_{g}^{\rm
inj})$ and $\tilde{h}_{L}(f;m_{g})\rightarrow\tilde{h}(f;m_{g})$, where
$\tilde{h}(f;m_{g})$ is defined by Eq. 4. As
$\sqrt{\braket{\tilde{h}_{1}}{\tilde{h}_{1}}}$ and
$\sqrt{\braket{\tilde{h}_{2}}{\tilde{h}_{2}}}$ are, respectively, the SNRs of
$\tilde{h}_{1}$ and $\tilde{h}_{2}$222Note that, throughout this work, SNR is
defined with respect to the Advanced LIGO and Virgo detectors at the design
sensitivity. , $\mathcal{M}$ does not depend on the SNR of the waveform
considered. Alternatively, $\mathcal{M}$ can be viewed as a normalized inner
product between the two waveforms, and its magnitude is always smaller than 1.
If $\tilde{h}_{1}$ and $\tilde{h}_{2}$ have more similarity, $\mathcal{M}$ is
closer to unity. In particular, if $\tilde{h}_{1}(f)\propto\tilde{h}_{2}(f)$,
meaning that $\tilde{h}_{1}$ and $\tilde{h}_{2}$ have the same morphology,
$\mathcal{M}=1$.

For the explicit calculations of $\mathcal{M}$, we consider:

  1. U.1

an unlensed waveform due to a GW150914-like source binary black hole LIGO_14
at a luminosity distance of 400 Mpc, whose SNR is 46,

  1. L.1

a lensed waveform of the unlensed waveform by an IMBH of reshifted mass of
$400M_{\odot}$ and impact parameter $y=0.9$, whose SNR is 57.

This mass of lens is chosen because IMBHs of similar masses are hoped to be
discovered by GW lensing Lensing_IMBH_01 . This value of $y$ is chosen because
IMBH lensing is more likely to occur at a larger $y$ (see the subsequent
discussion of the prior of $y$). The existing constraints on $m_{g}$ by GW
detection LIGO_12 ; LIGO_15 ; LIGO_16  suggest that we can probe the existence
of massive gravitons of $\sim 10^{-22}$ eV via GW detection. Thus, we consider
$m_{g}^{\rm inj}=10^{-22}\rm eV$.

Figure 2: The match, a function gauging the similarity between waveforms, of
lensed GWs and unlensed GWs as a function of the graviton mass. In particular,
we compare the similarity between lensed dispersive GWs of $m_{g}=10^{-22}\rm
eV$ to lensed dispersive GWs of other $m_{g}$ (solid blue) and the similarity
between unlensed dispersive GWs of $m_{g}=10^{-22}\rm eV$ to unlensed
dispersive GWs of other $m_{g}$. The match of lensed dispersive GWs shows a
narrower peak, suggesting that the waveform morphology of lensed dispersive
GWs vary more sensitively with $m_{g}$. As we shall see, this character of
lensed GWs can lead to a better measurement accuracy of $m_{g}$ over unlensed
GWs.

Fig. 2 shows the match $\mathcal{M}$ of the lensed waveform L.1 (solid blue
line) and the unlensed lensed waveform U.1 (dashed red line) as a function of
$m_{g}$. For both waveforms, $\mathcal{M}$ peaks at unity at the injected
$m_{g}$; this is reasonable because at the injected $m_{g}$,
$\tilde{h}_{1}(f)=\tilde{h}_{2}(f)$. Nonetheless, the match of the lensed
waveform peaks more sharply than that of the unlensed signal, suggesting that,
at the same SNR, the morphology of lensed dispersive GWs vary more sensitively
than the unlensed waves at $m_{g}\approx 10^{-22}\rm eV$. This feature hints
that lensing may have advantages in measuring $m_{g}$ over unlensed waves. For
$m_{g}$ far away from the injected $m_{g}$, $\mathcal{M}(m_{g})$ drops well
below unity for both the lensed and the unlensed waveform, because the
waveform at those values of $m_{g}$ is significantly different from that at
$m_{g}=10^{-22}\rm eV$. At last, we note that $\mathcal{M}$ for both lensed
and unlensed waves is oscillatory for $m_{g}$ far away from $m_{g}^{\rm inj}$.
This is because $\tilde{h}_{2}(f)\sim\tilde{h}_{1}(f)e^{i\left[\Psi_{\rm
disp}(f;m_{g}^{\rm inj})-\Psi_{\rm disp}(f;m_{g})\right]}$, where $\Psi_{\rm
disp}(f;m_{g})$ is defined by Eq. 3, leading to

$\braket{\tilde{h}_{1}}{\tilde{h}_{2}}\propto\int_{0}^{+\infty}df\frac{|\tilde{h}_{1}(f)|^{2}}{S_{n}(f)}\cos\left\\{\frac{\pi
D_{0}\left[m_{g}^{2}-\left(m_{g}^{\rm
inj}\right)^{2}\right]}{(1+z)f}\right\\}.$ (15)

The cosine term in the integrand makes $\mathcal{M}$ oscillatory over $m_{g}$.

##  III Parameter Estimation

###  III.1 Mock signals

To further investigate how lensing modifications of waveform morphology (R.1)
and SNR (R.2) may help to improve the measurement of $m_{g}$, we analyze a
mock signal of L.1 and U.1 that is injected into simulated Gaussian noises
assuming the design sensitivity of the Advanced LIGO and Virgo detectors
LIGO_noise_curve ; Virgo_noise_curve . We also inject

  1. L.2

a lensed signal which is identical to L.1 except the source binary is at 500
Mpc.

The luminosity distance of the source binary of L.2 is increased so that the
SNR of L.2 is the same as that of U.1. On the other hand, it is estimated that
the Advanced LIGO and Virgo detectors will detect 0.05 IMBH-lensed events per
year (or 1 IMBH-lensed event per $\sim$ 20 years) Lensing_IMBH_01 . At its
design sensitivity, the Advanced LIGO and Virgo are expected to detect $\leq
360$ unlensed events per year LIGO_detection_rate . Thus, a more fair
comparison will be with the posterior of $m_{g}$ combined across $20\times
360\sim$ 7000 unlensed signals. In practice, the combined measurement accuracy
of $m_{g}$ will be dominated by the signal with the best measurement accuracy,
which depends on the SNR of the signal m_g_bound_review_02 . Thus, we first
simulated a population of $\sim 7000$ binary black-hole mergers according to
LIGO_17 , each of which has an SNR of $\geq 10$, approximately the minimum SNR
for an event to be detectable by the Advanced LIGO and Virgo detectors LIGO_16
; LIGO_15 ; LIGO_sensitivity . Then, we inject the fourth signal, which is

  1. P.1

the unlensed signal that has the largest SNR (130) among the simulated 7000
unlensed events.

We represent the measurement of $m_{g}$ combined across these $7000$ simulated
signals by the posterior of $m_{g}$ of P.1.

###  III.2 Bayesian inference

We denote parameters describing the source binary by $\vec{\theta}_{\rm BBH}$
and parameters describing lensing by $\vec{\theta}_{\rm lens}=(M_{\rm
len},y)$. By Bayes’ theorem, the posterior of $m_{g}$, $\vec{\theta}_{\rm
len}$ and $\vec{\theta}_{\rm BBH}$ is given by

$\begin{split}&p(\vec{\theta}_{\rm BBH},\vec{\theta}_{\rm
lens},m_{g}|\tilde{d},H,I)\\\ &\propto p_{\rm BBH}(\vec{\theta}_{\rm
BBH}|H,I)p_{\rm lens}(\vec{\theta}_{\rm lens}|H,I)p_{m}(m_{g}|H,I)\\\
&\quad\times p(\tilde{d}|\vec{\theta}_{\rm BBH},\vec{\theta}_{\rm
lens},m_{g},H,I),\end{split}$ (16)

where $p_{\rm BBH}(\vec{\theta}_{\rm BBH}|H,I)$, $p_{\rm
lens}(\vec{\theta}_{\rm lens}|H,I)$ and $p_{m}(m_{g}|H,I)$ are, respectively,
the prior of $\vec{\theta}_{\rm BBH}$, $\vec{\theta}_{\rm lens}$ and $m_{g}$,
given the hypothesis $H$ that GWs may exhibit dispersion relation due to the
massive gravitons and background information $I$, such as that the signal is
lensed, the amplification function LABEL:eq:Amp_fn_mg and lensing geometry
etc. Since $\vec{\theta}_{\rm BBH}$, $\vec{\theta}_{\rm lens}$ and $m_{g}$
should be independent, we have assumed that their priors are factorized.
$p(\tilde{d}|\vec{\theta}_{\rm BBH},\vec{\theta}_{\rm lens},m_{g},H,I)$ is the
likelihood that a binary black hole of $\vec{\theta}_{\rm BBH}$ and lens of
$\vec{\theta}_{\rm len}$ will generate detected strain data $\tilde{d}_{D}$,

$\begin{split}&p(\tilde{d}|m_{g},\vec{\theta}_{\rm len},\vec{\theta}_{\rm
BBH},\vec{\theta},H,I)\propto\exp\left(-\frac{1}{2}\braket{\tilde{n}(f)}{\tilde{n}(f)}\right),\\\
&\tilde{n}(f)=\tilde{h}_{\rm D}(f;m_{g},\vec{\theta}_{\rm
len},\vec{\theta}_{\rm BBH})-\tilde{d}_{\rm D},\end{split}$ (17)

where $\tilde{h}_{\rm D}(m_{g},\vec{\theta}_{\rm len},and\vec{\theta}_{\rm
BBH})$ is the frequency-domain responses corresponding to detector $D$ by the
waveform equation Eq. 11.

Following Lensing_IMBH_01 , we place a uniform prior for $M_{\rm len}$. For
$y$, we place a prior which is uniform for $y^{2}\in[0,1]$ instead of $y$. For
$m_{g}$, we place a prior which is uniform for $\log_{10}m_{g}\in[-26,-20]$,
covering the magnitude of the most updated constraints on $m_{g}$ LIGO_15  by
GWs and for us to explore tighter constraints. At last, the marginalized
posterior of $m_{g}$ can be obtained by marginalizing Eq. 16 over
$\vec{\theta}_{\rm BBH}$ and $\vec{\theta}_{\rm lens}$,

$\begin{split}&p(m_{g}|\tilde{d},H,I)\\\ &=\int d\vec{\theta}_{\rm BBH}\int
d\vec{\theta}_{\rm lens}p(\vec{\theta}_{\rm BBH},\vec{\theta}_{\rm
lens},m_{g}|\tilde{d},H,I).\end{split}$ (18)

###  III.3 Mock signals of $m_{g}=0$

We first analyze U.1, L.1, L.2 and P.1 that are generated by assuming
$m_{g}=0$. The frequency-domain strains of U.1 and P.1 are generated using the
IMRPhenoPv2 template IMRPhenom_2011 ; IMRPhenom_2014 , a phenomenological
waveform template calibrated against numerical-relativity simulations, using
the LALSimulation library lalsuite . The simulated unlensed signals contain
the inspiral, merger, and ringdown phase. We then map U.1 into L.1 by
multiplying the frequency-domain waveform of U.1 by the amplification function
Eq. 6. L.2 is also generated according to these procedures. When inferring L.1
and L.2, we use the waveform model of Eq. 11 with the dephasing due to the
massive gravitons included and infer $m_{g}$ along with $\vec{\theta}_{\rm
BBH}$ and $\vec{\theta}_{\rm len}$. For U.1 and P.1, we infer with the
waveform model with $F(f;m_{g})$ in Eq. 11 set to be $1$ for all frequencies
and $M_{\rm len}$ and $y$ are removed from inference.

Figure 3: The corner plot shows the marginalized posterior of the redshifted
lens mass $M_{\rm len}$, $y$ and $\log m_{g}$ and their correlations, although
we also infer the parameters of the source binary together as free parameters.
The posteriors are estimated from a mock lensed signal due to a GW150914-like
binary lensed by a black hole of redshifted mass $M_{\rm len}$ of
$400M_{\odot}$ at $y=0.9$. The green lines denote the injected values for
$M_{\rm len}$ and $y$. The red line on the marginalized posterior of
$\log_{10}m_{g}$ denotes the $3\sigma$ confidence interval (CI) from the lower
limit of the prior of $\log_{10}m_{g}$. We conclude that we can bound the
graviton mass while accurately measuring the lensing-related parameters.

The diagonal of Fig. 3 shows the posterior of redshifted lens mass $M_{\rm
len}$, $y$ and $\log m_{g}$ obtained from L.1. The off-diagonal plots show the
two-dimensional posterior distributions among the variables. The green
vertical lines mark the injected values. The red vertical line marks the
$3\sigma$ interval of the marginalized posterior of $\log_{10}m_{g}$ from
$m_{g}=10^{-26}\rm eV$. From Fig. 3, we find that the posterior of
$\log_{10}m_{g}$ has no support for $\log_{10}m_{g}>-23.2$ because our
measurement of GWs rules out the possibility of an excessive large $m_{g}$.
From the posterior of $M_{\rm len}$ and $y$, we conclude that we can
accurately estimate the lensing-related parameters while testing the graviton
mass with lensing. Moreover, judging from Fig. 3 there are no strong
correlations between the lensing-related parameters and $m_{g}$.

We now compare the constraints on $m_{g}$ obtained from different
nondispersive GW signals. Fig. 4 shows the posterior of $\log_{10}m_{g}$ of
L.1 (solid blue line), its unlensed counterpart U.1 (dashed red line), L.2
(dashed dotted dotted black line) and P.1 (dashed dotted green). We notice
that all posteriors are in step-function shape because the measurement rules
out large values for $m_{g}$, which would produce discernible effects on the
waveform. All posteriors correspond to a similar $3\sigma$ confidence interval
(CI), ranging from $3.3\times 10^{-24}\rm eV$ to $1.3\times 10^{-23}\rm eV$.
In particular, L.1 yields a constraint (3$\sigma$ CI) on $m_{g}$ of $5.5\times
10^{-24}\rm eV$, slightly better than the constraint on $m_{g}$ by U.1,
corresponding to $1.3\times 10^{-23}\rm eV$. At the same SNR, we find that the
$3\sigma$ CI of L.2 is $1.3\times 10^{-23}\rm eV$, almost the same as that of
U.1. Even with large SNR, P.1 yields a constraint on $m_{g}$ of $\sim
3.3\times 10^{-24}\rm eV$, slightly better than the constraint by all the
other signals. These results conclude that lensing and increasing the SNR do
not significantly improve the constraints.

Figure 4: The marginalized posterior of $m_{g}$ in solid blue line is obtained
from a lensed signal (L.1) due to a GW150914-like binary lensed by a black
hole of redshifted mass of $400M_{\odot}$ at $y=0.9$. To gauge the improvement
of the constraint, we compare the result of a lensed signal with its unlensed
counterpart (U.1, dashed red line). To examine the role of the signal-to-noise
ratio (SNR), we inject another lensed signal (L.2, in dot-dot-dashed line)
that is identical to L.1 except that the source binary is more farther away so
that it has the same SNR of U.1. At last, lensing rate calculations suggest
that the Advanced LIGO and Virgo may detect an IMBH-lensed signal per $\sim$
7000 unlensed signals. For a fair comparison, we also compare the lensed
result with the best constraint likely to be obtained from 7000 simulated
unlensed signals (P.1, dotted green line, see its generation in the main
text). The solid vertical line denotes the 3$\sigma$ CI of the posterior of
U.1 ($1.3\times 10^{-23}\rm eV$) and the dashed vertical line denotes the
3$\sigma$-CI of L.1 ($5.5\times 10^{-24}\rm eV$). All posteriors correspond to
a similar 3-$\sigma$ interval, ranging from ranging from $3.3\times
10^{-24}\rm eV$ to $1.3\times 10^{-23}\rm eV$.

###  III.4 Mock signals of $m_{g}=10^{-22}\rm eV$

On the other hand, we find that lensing can help to improve the measurement of
$m_{g}$ from dispersive GWs. Fig. 5 shows the posterior of $\log_{10}m_{g}$
obtained from U.1 (dashed red line), L.1 (solid blue line), L.2 (dot-dot-
dashed line) and P.1 (dot-dashed green line) that are generated by assuming
$m_{g}=10^{-22}\rm eV$ (solid vertical black line). The shaded region
illustrates the $3\sigma$ CI of the posterior of L.1. The embed figure shows
the zoomed-in comparison of the posterior of L.1 and P.1. We notice that, at
the same SNR, lensing still improves the measurement accuracy of $m_{g}$ over
its unlensed counterpart. The posterior of L.2 shows more support for $m_{g}$
close to the injected $m_{g}$ than U.1, indicating that the the posterior of
L.2 is more accurate than that of U.1. This is because lensing modulates the
amplitude and phase of GWs, so that lensed GWs depend on $m_{g}$ more
sensitively, increasing the detectability of dispersive GWs, as indicated by
Fig. 2. The posterior of L.1 and P.1 peaks at a $m_{g}$ closer to the injected
$m_{g}$ because of larger SNR. Nevertheless, L.1 still leads to significant
improvement of the measurement accuracy of $m_{g}$ to an extent comparable to
P.1. From the results of Fig. 5, we find that both the lensing modifications
of waveform morphology and amplification can contribute to the improved
measurement accuracy of the graviton mass.

Figure 5: The marginalized posterior of $m_{g}$ obtained from L.1 (solid
blue), U.1 (dashed red), L.2 (dashed-dotted-dotted black line) and P.1
(dashed-dotted green). The embed figure shows the zoomed-in comparison of the
posterior of L.1 and P.1. The shaped region denotes the 3$\sigma$ confidence
interval of the posterior obtained from the lensed signal. For all signals, we
assume $m_{g}=10^{-22}\mathrm{eV}$ (vertical line in black). The lensing
geometry is the same as that considered in Fig. 4 and Section III.3. By
comparing the posterior of L.2 and U.1, we find that, even at the same SNR
lensing modifications of waveform morphology contribute to improving
measurement accuracy of $m_{g}$.

##  IV Conclusions

This paper studies the lensing pattern by a point-mass lens of GWs with an
isotropic dispersion relation due to the massive gravitons. Although the
graviton mass close to the existing constraints leads to no significant
effects on the lensing amplification, we find that lensing modifies the
waveform morphology of dispersive GWs, making the morphology changes more
sensitively with the graviton mass, which helps to improve the measurement of
the graviton mass. The improvement can also be further enhanced by the
increase of signal-to-noise ratio due to lensing. By detecting a lensed
gravitational-wave signal, we can measure the graviton mass with an accuracy
comparable with the combined measurement across $\mathcal{O}(10^{3})$ unlensed
signals. Our work lays the foundation for measuring the graviton mass in the
era of detectable lensed GWs, making the existing analyses that focus
primarily on unlensed signals more complete.

Other than the improvement of measurement accuracy of the graviton mass, our
method enjoys several advantages. First, compared to other proposed methods of
testing the speed of GWs by observing lensing Collett_Bacon_2017 ; Fan_2017 ,
our approach requires no observation of the electromagnetic counterpart(s) of
a given event. Therefore, our method is more stand-alone and is easier to be
performed. Second, our method is independent of the nature of the source
binaries. Although in this paper, we focused on GWs generated by binary black
holes, our method can be straightforwardly applied to other types of
coalescence, such as binary neutron star coalescence Lensing_GW_BNS . This
flexibility greatly extends the scope of graviton-mass measurement. Lastly,
our method makes the test of graviton mass more complete. While the far-field
propagation of GWs Will_2012 ; Ajith_2010  and near-field behavior of black
holes Chung_2019  have been proposed to constrain the mass of graviton, our
test bridges the intermediate region between these two tests. Along with other
tests of general relativity via observing the lensing of GWs (such as
Lensing_TGR ), our test demonstrates the strong potential to understand the
nature of space-time via observing gravitational-wave lensing.

In this work, we ignore the effects of (i) the change of polarization of GWs
due to lensing Lensing_polarization ; (ii) the change of the behavior of the
source compact binary due to massive graviton, as is the case in LIGO_07 ;
LIGO_11 ; and (iii) the change of the gravitational field around the lens by
the graviton mass. Also, our study focusing on the case of point-mass lens.
These ignored effects and the lensing of dispersive GWs of other lens types
remain fully explored. If we include these effects in our measurement, the
accuracy can be further enhanced.

In the future, we plan to extend our studies to other types of lenses, which
may help further improve the measurement accuracy. Our study has thus far
focused on the point-mass lens, such as intermediate-mass black holes, which
leads to microlensing. In reality, it may be very rare for gravitational waves
to be lensed by an intermediate-mass black hole of $\sim 400M_{\odot}$.
Moreover, the population properties and lensing rates of intermediate-mass
black holes are uncertain. On the other hand, strong lensing due to different
types of lenses, such as galaxies or galaxy clusters galaxy_strong_lensing_01
, are expected to be more common, roughly 1 per $\sim 600$ unlensed events at
the design sensitivity of LIGO and Virgo Ken_2016 . Upon strong lensing, a GW
signal may split into multiple images whose properties, such as image position
and the arrival time differences, may depend on the graviton mass even more
sensitively than the diffraction pattern considered in this work Baker ;
private_communication . To extend our test to strong lensing, we need to study
the strong lensing of dispersive GWs by lenses with structures, such as
galactic lenses, singular isothermal sphere, and other possible extended mass
distribution lensing_mass_distribution . We would also like to investigate the
performance of our test for the detection by proposed space-based detectors,
such as the Laser Interferometer Space Antenna LISA , which are capable of
exquisite phase measurement and much better constraints. Therefore, in the
future, we can measure the graviton mass with unparalleled accuracy by
observing lensed gravitational-wave signals.

## Acknowledgements

The authors are indebted to valuable discussion among the lensing working
group of LIGO. A.K.-W.C. would like to acknowledge Patrick C.K. Cheong,
Srashti Goyal and Shasvath Kapadia for stimulating discussions, Jose Maria
Ezquiaga Bravo, Mark H.Y. Cheung, Otto A. Hannuksela, Alvin K.Y. Li and
Ignacio Magana for their comments on the manuscript and relevant presentations
and Robin S.H. Yuen for his advice about computer programming. A.K.-W.C. was
supported by the Hong Kong Scholarship for Excellence Scheme (HKSES). The work
described in this paper was partially supported by grants from the Research
Grants Council of the Hong Kong (Project No. CUHK 24304317 and CUHK 14306218),
The Croucher Foundation of Hong Kong, and the Research Committee of the
Chinese University of Hong Kong. This manuscript carries a report number of
KCL-PH-TH 2021/41 and LIGO Document number of P2100192-v2.

This research has made use of data, software and/or web tools obtained from
the GW Open Science Center (https://www.gw-openscience.org), a service of LIGO
Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration.
LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the
French Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with
contributions by Polish and Hungarian institutes.

## References

  * (1) https://dcc.ligo.org/LIGO-P1200087-v42/public.
  * (2) https://dcc.ligo.org/LIGO-T2000012/public.
  * (3) P. Ajith, M. Hannam, S. Husa, Y. Chen, B. Brügmann, N. Dorband, D. Müller, F. Ohme, D. Pollney, C. Reisswig, L. Santamaría, and J. Seiler.  Inspiral-merger-ringdown waveforms for black-hole binaries with nonprecessing spins.  Phys. Rev. Lett., 106:241101, Jun 2011.
  * (4) B. P. Abbott et al. .  Binary black hole mergers in the first advanced ligo observing run.  Phys. Rev. X, 6:041015, Oct 2016.
  * (5) B. P. Abbott et al. .  Gw151226: Observation of gravitational waves from a 22-solar-mass binary black hole coalescence.  Phys. Rev. Lett., 116:241103, Jun 2016.
  * (6) B. P. Abbott et al. .  Observation of gravitational waves from a binary black hole merger.  Phys. Rev. Lett., 116:061102, Feb 2016.
  * (7) B. P. Abbott et al. .  Properties of the Binary Black Hole Merger GW150914.  prl, 116(24):241102, June 2016.
  * (8) B. P. Abbott et al. .  Tests of general relativity with gw150914.  Phys. Rev. Lett., 116:221101, May 2016.
  * (9) B. P. Abbott et al. .  Gw170104: Observation of a 50-solar-mass binary black hole coalescence at redshift 0.2.  Phys. Rev. Lett., 118:221101, Jun 2017.
  * (10) B. P. Abbott et al. .  Gw170608: Observation of a 19 solar-mass binary black hole coalescence.  The Astrophysical Journal Letters, 851(2):L35, 2017.
  * (11) B. P. Abbott et al. .  Gw170814: A three-detector observation of gravitational waves from a binary black hole coalescence.  Phys. Rev. Lett., 119:141101, Oct 2017.
  * (12) B. P. Abbott et al. .  Gw170817: Observation of gravitational waves from a binary neutron star inspiral.  Phys. Rev. Lett., 119:161101, Oct 2017.
  * (13) B. P. Abbott et al. .  Prospects for Observing and Localizing Gravitational-Wave Transients with Advanced LIGO, Advanced Virgo and KAGRA.  Living Rev. Rel., 21:3, 2018.  [Living Rev. Rel.19,1(2016)].
  * (14) B. P. Abbott et al. .  Gwtc-1: A gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs.  Phys. Rev. X, 9:031040, Sep 2019.
  * (15) B. P. Abbott et al. .  Tests of general relativity with the binary black hole signals from the ligo-virgo catalog gwtc-1.  Phys. Rev. D, 100:104036, Nov 2019.
  * (16) B. P. Abbott et al. .  GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses.  arXiv e-prints, page arXiv:2004.08342, April 2020.
  * (17) B. P. Abbott et al. .  GW190814: Gravitational waves from the coalescence of a 23 solar mass black hole with a 2.6 solar mass compact object.  The Astrophysical Journal, 896(2):L44, jun 2020.
  * (18) B. P. Abbott et al. .  GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run.  arXiv e-prints, page arXiv:2010.14527, October 2020.
  * (19) B. P. Abbott et al. .  Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog.  arXiv e-prints, page arXiv:2010.14533, October 2020.
  * (20) B. P. Abbott et al. .  Tests of General Relativity with Binary Black Holes from the second LIGO-Virgo Gravitational-Wave Transient Catalog.  arXiv e-prints, page arXiv:2010.14529, October 2020.
  * (21) B. P. Abbott et al. .  Search for lensing signatures in the gravitational-wave observations from the first half of LIGO-Virgo’s third observing run.  arXiv e-prints, page arXiv:2105.06384, May 2021.
  * (22) Vishal Baibhav, Emanuele Berti, Davide Gerosa, Michela Mapelli, Nicola Giacobbo, Yann Bouffanais, and Ugo N. Di Carlo.  Gravitational-wave detection rates for compact binaries formed in isolation: Ligo/virgo o3 and beyond.  Phys. Rev. D, 100:064060, Sep 2019.
  * (23) Tessa Baker and Mark Trodden.  Multimessenger time delays from lensed gravitational waves.  Phys. Rev. D, 95:063512, Mar 2017.
  * (24) Haugan Mark P Bontz, Robert J.  A diffraction limit on the gravitational lens effect.  Astrophysics and Space Science, August 1981.
  * (25) Alessandra Buonanno, Bala R. Iyer, Evan Ochsner, Yi Pan, and B. S. Sathyaprakash.  Comparison of post-newtonian templates for compact binary inspiral signals in gravitational-wave detectors.  Phys. Rev. D, 80:084043, Oct 2009.
  * (26) Adrian Ka-Wai Chung and Tjonnie Guang Feng Li.  Phenomenological inclusion of alternative dispersion relations to the teukolsky equation and its application to bounding the graviton mass with gravitational-wave measurements.  Phys. Rev. D, 99:124023, Jun 2019.
  * (27) T. E. Collett and D. Bacon.  Testing the Speed of Gravitational Waves over Cosmological Distances with Strong Gravitational Lensing.  Physical Review Letters, 118(9):091101, March 2017.
  * (28) Paolo Cremonese, Jose María Ezquiaga, and Vincenzo Salzano.  Breaking the mass-sheet degeneracy with gravitational wave interference in lensed events.  arXiv e-prints, page arXiv:2104.07055, April 2021.
  * (29) Claudia de Rham.  Massive Gravity.  Living Reviews in Relativity, 17(1):7, December 2014.
  * (30) Claudia de Rham, J. Tate Deskins, Andrew J. Tolley, and Shuang-Yong Zhou.  Graviton mass bounds.  Rev. Mod. Phys., 89:025004, May 2017.
  * (31) S. Deguchi and W. D. Watson.  Diffraction in Gravitational Lensing for Compact Objects of Low Mass.  Astrophys. J. , 307:30, August 1986.
  * (32) P. Schneider J. Ehlers E.E. Falco.  Gravitational Lenses.  Springer-Verlag Berlin Heidelberg, 1992.
  * (33) X.-L. Fan, K. Liao, M. Biesiada, A. Piórkowska-Kurpas, and Z.-H. Zhu.  Speed of Gravitational Waves from Strongly Lensed Gravitational Waves and Electromagnetic Signals.  Physical Review Letters, 118(9):091102, March 2017.
  * (34) Éanna É. Flanagan and Scott A. Hughes.  Measuring gravitational waves from binary black hole coalescences. i. signal to noise for inspiral, merger, and ringdown.  Phys. Rev. D, 57:4535–4565, Apr 1998.
  * (35) Srashti Goyal, K. Haris, Ajit Kumar Mehta, and Parameswaran Ajith.  Testing the nature of gravitational-wave polarizations using strongly lensed signals.  arXiv e-prints, page arXiv:2008.07060, August 2020.
  * (36) Srashti Goyal and Shasvath Kapadia.  private communication, May 2021.
  * (37) O. A. Hannuksela, K. Haris, K. K. Y. Ng, S. Kumar, A. K. Mehta, D. Keitel, T. G. F. Li, and P. Ajith.  Search for Gravitational Lensing Signatures in LIGO-Virgo Binary Black Hole Events.  apjl, 874(1):L2, March 2019.
  * (38) Otto A. Hannuksela, Thomas E. Collett, Mesut Çalıskan, and Tjonnie G. F. Li.  Localizing merging black holes with sub-arcsecond precision using gravitational-wave lensing.  mnras, 498(3):3395–3402, August 2020.
  * (39) D. Hansen, N. Yunes, and K. Yagi.  Projected constraints on Lorentz-violating gravity with gravitational waves.  Phys. Rev. D, 91(8):082003, April 2015.
  * (40) Shaoqi Hou, Xi-Long Fan, and Zong-Hong Zhu.  Gravitational lensing of gravitational waves: Rotation of polarization plane.  Phys. Rev. D, 100:064028, Sep 2019.
  * (41) Tjonnie G. F. Li Ken K. Y. Ng, Kaze W. K. Wong and Tom Broadhurst.  Precise ligo lensing rate predictions for binary black holes.  Phys. Rev. D, 2017.
  * (42) D. Keppel and P. Ajith.  Constraining the mass of the graviton using coalescing black-hole binaries.  Phys. Rev. D, 82:122001, Dec 2010.
  * (43) Kwun-Hang Lai, Otto A. Hannuksela, Antonio Herrera-Martín, Jose M. Diego, Tom Broadhurst, and Tjonnie G. F. Li.  Discovering intermediate-mass black hole lenses through gravitational wave lensing.  Phys. Rev. D, 98(8):083005, October 2018.
  * (44) LIGO Scientific Collaboration.  LIGO Algorithm Library - LALSuite.  free software (GPL), 2018.
  * (45) Jose María Ezquiaga, Daniel E. Holz, Wayne Hu, Macarena Lagos, and Robert M. Wald.  Phase effects from strong gravitational lensing of gravitational waves.  arXiv e-prints, page arXiv:2008.12814, August 2020.
  * (46) Jose María Ezquiaga and Miguel Zumalacárregui.  Gravitational wave lensing beyond general relativity: birefringence, echoes and shadows.  arXiv e-prints, page arXiv:2009.12187, September 2020.
  * (47) S. Mirshekari, N. Yunes, and C. M. Will.  Constraining Lorentz-violating, modified dispersion relations with gravitational waves.  Phys. Rev. D, 85(2):024041, January 2012.
  * (48) C. J. Moore, R. H. Cole, and C. P. L. Berry.  Gravitational-wave sensitivity curves.  Classical and Quantum Gravity, 32:015014, Jan 2015.
  * (49) Suvodip Mukherjee, Tom Broadhurst, Jose M. Diego, Joseph Silk, and George F. Smoot.  Inferring the lensing rate of LIGO-Virgo sources from the stochastic gravitational wave background.  mnras, 501(2):2451–2466, February 2021.
  * (50) Suvodip Mukherjee, Benjamin D. Wandelt, and Joseph Silk.  Multimessenger tests of gravity with weakly lensed gravitational waves.  Phys. Rev. D, 101(10):103509, May 2020.
  * (51) Suvodip Mukherjee, Benjamin D. Wandelt, and Joseph Silk.  Probing the theory of gravity with gravitational lensing of gravitational waves and galaxy surveys.  mnras, 494(2):1956–1970, May 2020.
  * (52) Takahiro T. Nakamura and Shuji Deguchi.  Wave optics in gravitational lensing.  Progress of Theoretical Physics Supplement, 133:137–153, 1999.
  * (53) Hans C. Ohanian.  On the focusing of gravitational radiation.  International Journal of Theoretical Physics, June 1974.
  * (54) G. Pagano, O. A. Hannuksela, and T. G. F. Li.  LENSINGGW: a PYTHON package for lensing of gravitational waves.  aap, 643:A167, November 2020.
  * (55) Yi Pan, Alessandra Buonanno, Andrea Taracchini, Lawrence E. Kidder, Abdul H. Mroué, Harald P. Pfeiffer, Mark A. Scheel, and Béla Szilágyi.  Inspiral-merger-ringdown waveforms of spinning, precessing black-hole binaries in the effective-one-body formalism.  Phys. Rev. D, 89:084006, Apr 2014.
  * (56) Peter T. H. Pang, Otto A. Hannuksela, Tim Dietrich, Giulia Pagano, and Ian W. Harry.  Lensed or not lensed: determining lensing magnifications for binary neutron star mergers from a single detection.  mnras, 495(4):3740–3750, May 2020.
  * (57) Scott Perkins and Nicolás Yunes.  Probing screening and the graviton mass with gravitational waves.  Classical and Quantum Gravity, 36(5):055013, March 2019.
  * (58) Travis Robson, Neil J Cornish, and Chang Liug.  The construction and use of lisa sensitivity curves.  Classical and Quantum Gravity, 36(10):105011, 2019.
  * (59) P. Schneider.  Introduction to Gravitational Lensing and Cosmology, pages 1–89.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2006.
  * (60) Graham P. Smith, Christopher Berry, Matteo Bianconi, Will M. Farr, Mathilde Jauzac, Richard Massey, Johan Richard, Andrew Robertson, Keren Sharon, Alberto Vecchio, and et al.  Strong-lensing of gravitational waves by galaxy clusters.  Proceedings of the International Astronomical Union, 13(S338):98–102, 2017.
  * (61) R. Takahashi.  Arrival Time Differences between Gravitational Waves and Electromagnetic Signals due to Gravitational Lensing.  Astrophys. J. , 835:103, January 2017.
  * (62) Ryuichi Takahashi and Takashi Nakamura.  Wave Effects in the Gravitational Lensing of Gravitational Waves from Chirping Binaries.  Astrophys. J. , 595(2):1039–1051, October 2003.
  * (63) Ryuichi Takahashi and Takashi Nakamura.  Wave effects in the gravitational lensing of gravitational waves from chirping binaries.  The Astrophysical Journal, 595(2):1039, 2003.
  * (64) Slava G. Turyshev and Viktor T. Toth.  Gravitational lensing by an extended mass distribution.  arXiv e-prints, page arXiv:2106.06696, June 2021.
  * (65) Clifford M. Will.  Bounding the mass of the graviton using gravitational-wave observations of inspiralling compact binaries.  Phys. Rev. D, 57:2061–2068, Feb 1998.


# Monte Carlo calculation of the average neutron depolarization for the
NPDGamma experiment

K. B. Grammer<EMAIL_ADDRESS>J. D. Bowman  University of Tennessee,
Knoxville, TN, USA  Oak Ridge National Laboratory, Oak Ridge, TN, USA

###### Abstract

The NPDGamma experiment measures the asymmetry in $\gamma$-ray emission in the
capture of polarized neutrons on liquid parahydrogen. The beam polarization is
measured using 3He spin analysis, but this measurement does not account for
the contribution of depolarization from spin-flip scattering primarily due to
orthohydrogen in the bulk liquid. This is a systematic effect that dilutes the
experimental asymmetry and is modeled using Monte Carlo. Methods for tracking
neutron spin in MCNPX were developed in order to calculate the average neutron
polarization upon capture for use as a multiplicative correction to the
measured beam polarization for the NPDGamma
experiment.mytitlenotemytitlenotefootnotetext: This manuscript has been
authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the
U.S. Department of Energy. The United States Government retains and the
publisher, by accepting the article for publication, acknowledges that the
United States Government retains a non-exclusive, paid-up, irrevocable,
worldwide license to publish or reproduce the published form of this
manuscript, or allow others to do so, for United States Government purposes.
The Department of Energy will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
(http://energy.gov/downloads/doe-public-access-
plan).titlenote2titlenote2footnotetext: © 2019\. This manuscript version is
made available under the CC-BY-NC-ND 4.0 license
http://creativecommons.org/licenses/by-nc-nd/4.0/

###### keywords:

Monte Carlo, MCNPX, NPDGamma experiment

##  1 Introduction

The NPDGamma experiment [1] uses a polarized neutron beam in order to measure
the parity violating asymmetry in the angular distribution of the emitted 2.2
MeV $\gamma$-rays from neutron capture on liquid parahydrogen [2]

$\vec{\textrm{n}}+\textrm{p}\to\textrm{d}+\gamma.$ (1)

The physics asymmetry, $A_{\gamma}$, is defined by

$\frac{d\sigma}{d\Omega}\propto\frac{1}{4\pi}(1+A_{\gamma}\hat{\sigma}_{\textrm{n}}\cdot\hat{k}_{\gamma}).$
(2)

where $\hat{\sigma}_{\textrm{n}}$ is the neutron spin, and $\hat{k}_{\gamma}$
is the $\gamma$-ray momentum. Neutrons exit the $10\times 12$ cm neutron guide
on the Fundamental Neutron Physics Beamline (FNPB) [3] and enter a supermirror
polarizer [4], after which is the NPDGamma apparatus [5, 6]. Polarized
neutrons then enter the resonant frequency spin rotator (RFSR) [7] and are
rotated according to a pattern
($\uparrow\downarrow\downarrow\uparrow\downarrow\uparrow\uparrow\downarrow$)
that cancels beam power fluctuations to second order and the signals from
opposite spin states are used to isolate the asymmetry signal. The neutrons
are incident on the 16-liter liquid hydrogen target [8] that is surrounded by
48 CsI(Tl) detectors arranged in 4 rings of 12 detectors each [9] that detect
the $\gamma$-rays from neutron capture. The average neutron beam polarization
is measured with an empty liquid hydrogen target vessel, and a 3He spin filter
[10] followed by a 3He transmission monitor [11] downstream of the vessel are
used to measure the beam polarization. A uniform 9.5 gauss magnetic field
aligned within 3 mrad to the vertical axis ($\hat{z}$) transports polarized
neutrons to the target vessel [1].

The raw asymmetry for a given spin sequence is determined using the super
ratio method for each pair of detectors, $i$, where a pair of opposed
detectors, $l=0\ldots 5$ and $m=l+6$, is separated by $180^{\circ}$,

$A_{\textrm{raw}}^{i}=\frac{\sqrt{\alpha_{i}}-1}{\sqrt{\alpha_{i}}+1},$ (3)

where $\alpha_{i}$ is the super ratio given by

$\alpha_{i}=\frac{N_{\uparrow}^{l}}{N_{\downarrow}^{l}}\frac{N_{\downarrow}^{m}}{N_{\uparrow}^{m}},$
(4)

and $N_{\uparrow}$ and $N_{\downarrow}$ are the spin up and spin down signals,
respectively. The physics asymmetry [2] for a given detector is given by

$A_{\gamma}=\frac{A_{\textrm{raw}}^{i}-A_{\textrm{app}}^{i}}{P_{\textrm{n}}\Delta^{i}_{\textrm{dep}}(\lambda)\Delta_{\textrm{sf}}},$
(5)

where $A_{\textrm{app}}^{i}$ is a term encompassing apparatus asymmetries
including noise and beam fluctuations, $P_{\textrm{n}}$ is the measured
neutron beam polarization, $\Delta_{\textrm{sf}}$ is the spin flip efficiency,
and $\Delta^{i}_{\textrm{dep}}(\lambda)$ is the detector-dependent neutron
depolarization correction. The depolarization correction depends on the
geometry of the target, detector array, and neutron beam and on the neutron
wavelength. Since the detector rings are azimuthally symmetric about the
target, the polarization correction can be expressed as a function of detector
rings ($\Delta^{r}_{\textrm{dep}}(\lambda)$) rather than a function of
detectors or detector pairs.

There are two spin isomers of the hydrogen molecule: parahydrogen
$(J=0,2,4\ldots)$ and orthohydrogen $(J=1,3,5\ldots)$. Room temperature
hydrogen gas has approximately a 3:1 ortho-to-para ratio, with the equilibrium
parahydrogen concentration increasing as the temperature decreases [12]. The
energy separation between the $J=0$ ground state of parahydrogen and the $J=1$
lowest energy orthohydrogen state is 14.7 meV and higher order rotational
energy levels in the hydrogen molecule are given by

$E_{\textrm{rot}}(J)=CJ(J+1),$ (6)

where $C=7.36\textrm{~{}meV}$ is the rotational constant for the hydrogen
molecule rigid rotor and $J$ is the angular momentum of the hydrogen molecule.
The energy separation between the lowest rotational energy levels in the
hydrogen molecule is used as the criterion for discerning scattering event
types in MCNPX. The hydrogen target vessel is initially filled with normal
hydrogen gas. The hydrogen gas is then liquefied and operated at 15.6 K and
the slow natural conversion process is accelerated by circulating the liquid
through an ortho-para converter (OPC) [13] containing 150 ml of hydrous iron
(III) oxide 30–50 mesh powder[14] as a catalyst. Production data is taken
after the hydrogen in the vessel has been condensed and circulated through the
OPC over the course of weeks such that the ortho-para ratio has reached a
steady state and is assumed to be constant throughout the bulk liquid.

The orthohydrogen scattering cross section for neutrons with energy 1-10 meV
(figure 1) is approximately 2 orders of magnitude larger than both the
parahydrogen scattering cross section [15] and the hydrogen absorption cross
section [16]. The transmission of neutrons through the target as well as the
average polarization of neutrons as they propagate through the target is
strongly sensitive to the orthohydrogen concentration of the bulk hydrogen.
Neutrons can undergo spin flip scattering in the apparatus components and from
scattering on orthohydrogen in the bulk liquid hydrogen. Since there is no
reliable means of measuring the spin flip scattering rate, the depolarization
effect is modeled using Monte Carlo.

Figure 1: Neutron scattering cross sections on parahydrogen and orthohydrogen
at 20 K from ENDF-VII [15] and the neutron absorption cross section [16].

There was no in situ method to determine the relative orthohydrogen to
parahydrogen concentration for the NPDGamma experiment. The methods to monitor
the orthohydrogen concentration available at the time of the construction of
the apparatus required feeding electrical signals out of the cryogenic region
and hydrogen safety requirements limited the number of wires in the vessel.
However, conversion towards a steady state condition with a time constant on
the order of days is observed via measurement of the relative neutron
transmission and produced a measurement of the parahydrogen scattering cross
section [17]. The cross section measurement also yielded a determination of an
upper bound (0.15%) on the orthohydrogen concentration in the bulk liquid
while the lower limit (0.015%) is given by thermodynamic equilibrium at 15.6
K. The upper and lower bounds on the orthohydrogen concentration are used to
as inputs for the calculation of $\Delta^{r}_{\textrm{dep}}$.

The polarization correction, $\Delta^{r}_{\textrm{dep}}$ in equation 5, is
calculated with two methods: the first uses a modified version of MCNPX 2.7
[18] and the second uses the Young and Koppel free gas model [19]. These two
models will be described next and the two calculation methods are combined to
determine the geometry and wavelength dependent polarization for the NPDGamma
analysis.

##  2 Spin-flip scattering formalism

We follow the spin-flip scattering formalism outlined in Moon et al. [20]. The
double differential cross section for a scattering process with momentum
transfer $\vec{K}=\vec{k}-\vec{k}^{\prime}$ that changes the neutron spin from
$s$ to $s^{\prime}$ and the scattering system state from $q$ to $q^{\prime}$
is given by [20]

$\begin{array}[]{rcl}\frac{d^{2}\sigma^{ss^{\prime}}}{d\Omega^{\prime}dE^{\prime}}&=&\sum_{q}P_{q}\sum_{q^{\prime}}\frac{k^{\prime}}{k}\|\langle
q^{\prime}|\sum_{i}e^{\imath\vec{K}\cdot\vec{r_{i}}}U_{i}^{ss^{\prime}}|q\rangle\|^{2}\\\
&&\times\delta(\Delta E_{\textrm{n}}+\Delta E_{q}),\\\ \end{array}$ (7)

where $P_{q}$ is the probability the system is in state $q$, with a sum over
the atomic sites $i$, and the $\delta$ function ensures energy conservation
for the change in energy of the neutron ($\Delta E_{\textrm{n}}$) and the
scattering system ($\Delta E_{q}$). The scattering amplitude is given by [20],

$U_{i}^{ss^{\prime}}=\langle s^{\prime}|a_{\textrm{coh},i}-p_{i}\vec{S_{\perp
i}}\cdot\vec{\sigma}+a_{\textrm{inc},i}\vec{I_{i}}\cdot{\sigma}|s\rangle,$ (8)

which has a coherent nuclear scattering term ($a_{\textrm{coh},i}$), a
magnetic scattering term ($p_{i}\vec{S_{\perp i}}\cdot\vec{\sigma}$), and an
incoherent nuclear scattering term
($a_{\textrm{inc},i}\vec{I_{i}}\cdot{\sigma})$. There are then four
combinations of $s$ and $s^{\prime}$ for neutrons with spin in the $\hat{z}$
direction [20],

$\begin{array}[]{rcl}U^{++}&=&a_{\textrm{coh},i}-p_{i}S_{\perp
zi}+a_{\textrm{inc},i}I_{zi},\\\ U^{--}&=&a_{\textrm{coh},i}+p_{i}S_{\perp
zi}-a_{\textrm{inc},i}I_{zi},\\\ U^{+-}&=&-p_{i}(S_{\perp xi}+\imath S_{\perp
yi})+a_{\textrm{inc},i}(I_{\perp xi}+\imath I_{\perp yi}),\\\
U^{-+}&=&-p_{i}(S_{\perp xi}-\imath S_{\perp yi})+a_{\textrm{inc},i}(I_{\perp
xi}-\imath I_{\perp yi}),\\\ \end{array}$ (9)

where $a_{\textrm{coh}}$ is the bound coherent nuclear scattering length,
$p_{i}$ is the magnetic scattering amplitude, $\vec{S}_{\perp}$ is the
projection the neutron spin onto the plane perpendicular to the scattering
vector, $\vec{I}$ is the nuclear spin, and $a_{\textrm{inc}}$ is the bound
incoherent scattering length.

As described by Moon[20], coherent scattering and random isotopic ordering
does not contribute to spin-flip scattering. While magnetic scattering
contributes to spin-flip scattering in general, it is not relevant for the
nonmagnetic hydrogen target in the NPDGamma experiment. Thus, only the
incoherent scattering term contributes to spin-flip scattering. The
differential scattering cross section per atom becomes simple for the spin-
flip and non-spin-flip cases,

$\frac{d\sigma^{+-}}{d\Omega}=\frac{d\sigma^{-+}}{d\Omega}=\frac{2}{3}a_{\textrm{inc}}^{2}I(I+1),$
(10)
$\frac{d\sigma^{++}}{d\Omega}=\frac{d\sigma^{--}}{d\Omega}=\frac{1}{3}a_{\textrm{inc}}^{2}I(I+1),$
(11)

and which should be compared to equation 15 below. The spin flip terms for a
neutron with spin along the $\hat{z}$ axis, the vertical polarization axis in
the NPDGamma experiment, contain only $S_{\perp xi}$ and $S_{\perp yi}$ terms,
such that only nuclear spins perpendicular to the neutron spin can cause a
spin flip and this gives rise to the $\frac{2}{3}$ factor in equation 10. The
probability that any given particle interaction leads to a spin flip is given
by

$P_{\textrm{flip}}=\frac{2}{3}\frac{\sigma_{\textrm{incoh}}}{\sigma_{\textrm{incoh}}+\sigma_{\textrm{coh}}+\sigma_{\textrm{abs}}},$
(12)

where $\sigma_{\textrm{incoh}}$ is the incoherent scattering cross section,
$\sigma_{\textrm{coh}}$ is the coherent scattering cross section, and
$\sigma_{\textrm{abs}}$ is the absorption cross section.

###  2.1 Neutron-hydrogen spin-flip scattering

Transition | $P_{\textrm{flip}}$
---|---
H${}_{\textrm{ortho}}\to$H${}_{\textrm{ortho}}$ | 0.659(1)
H${}_{\textrm{ortho}}\to$H${}_{\textrm{para}}$ | 1.0
H${}_{\textrm{para}}\to$H${}_{\textrm{para}}$ | 0.0
H${}_{\textrm{para}}\to$H${}_{\textrm{ortho}}$ | 1.0
Table 1: Molecular hydrogen scattering interactions with spin-flip
probabilities as implemented in our MCNPX model using ENDF/B-VII cross
sections [15]. The uncertainty for
H${}_{\textrm{ortho}}\to$H${}_{\textrm{ortho}}$ is due to the neutron energy
spectrum.

There are four possibilities of initial and final states for neutron
scattering from a hydrogen molecule (see table 1).
H${}_{\textrm{para}}\to$H${}_{\textrm{para}}$ scattering events do not change
the spin of the neutron for the neutron energies relevant to the NPDGamma
experiment. All H${}_{\textrm{para}}\to$H${}_{\textrm{ortho}}$ and
H${}_{\textrm{ortho}}\to$H${}_{\textrm{para}}$ scattering events are inelastic
in nature and necessarily cause a spin flip. A neutron with energy well below
14.7 meV does not have sufficient energy to undergo an inelastic spin-flip
scatter from a parahydrogen molecule such that all interactions with
parahydrogen do not change the neutron spin, unless the neutron has been
upscattered above the 14.7 meV threshold. An
H${}_{\textrm{ortho}}\to$H${}_{\textrm{ortho}}$ scattering event may cause a
spin-flip if the molecular spin is orthogonal to the neutron spin, which
follows the spin flip probability in equation 12. The Young-Koppel gas model
as well as the Keinert-Sax liquid model were used separately in order to
estimate the contributions of these scattering processes to the depolarization
of neutrons for the NPDGamma experiment.

####  2.1.1 Young-Koppel free gas model treatment

Young and Koppel [19] give analytical expressions for the differential
scattering cross sections for orthohydrogen and parahydrogen gases with the
inclusion of spin correlations, vibrations, and rotations. The momentum
transfer is given by $\vec{K}=\vec{k}-\vec{k}^{\prime}$, and the recoil energy
of the molecule is given by,

$E_{\textrm{rec}}=\frac{m_{\textrm{n}}}{m_{\textrm{m}}}\big{(}E_{\textrm{n}}-2\sqrt{E_{\textrm{n}}E_{\textrm{n}}^{\prime}}(\hat{k}\cdot\hat{k}^{\prime})+E_{\textrm{n}}^{\prime}\big{)}=\frac{K^{2}\hbar^{2}}{4m_{\textrm{n}}},$
(13)

where $E_{\textrm{n}}$ and $E_{\textrm{n}}^{\prime}$ are the initial and final
neutron energies, $\hat{k}$ and $\hat{k}^{\prime}$ are the initial and final
neutron velocity vectors, $m_{\textrm{n}}$ is the neutron mass, and
$m_{\textrm{m}}$ is the mass of the molecule. The energy balance for the
scattering event is given by,

$\Delta
E_{\textrm{n}}=E_{\textrm{n}}+E_{\textrm{rot}}(j_{\textrm{i}})-(E_{\textrm{n}}^{\prime}+E_{\textrm{vib}}(n)+E_{\textrm{rot}}(j_{\textrm{f}})+E_{\textrm{rec}}),$
(14)

where $E_{\textrm{vib}}(n)$ is the vibrational energy of the molecule in the
$n$th mode. At low temperatures, the molecular vibrational states can be
ignored because molecules are in the vibrational ground state and cold
neutrons do not have sufficient energy to excite higher order states. The Y-K
model ignores intermolecular interactions and the cross section is given by
the molecular “self”-term. The simplified cross sections (for $J=0,1$ )for the
spin isomer transitions in hydrogen can be expressed analytically as follows
from Zoppi [21][22],

$\begin{array}[]{ll}\frac{d^{2}\sigma}{d\Omega dE^{\prime}}|_{o\to
o}=&2C(a_{\textrm{coh}}^{2}+\frac{2}{3}a_{\textrm{inc}}^{2})(J_{0}(\alpha)^{2}+2J_{2}(\alpha)^{2}),\\\
\frac{d^{2}\sigma}{d\Omega dE^{\prime}}|_{o\to
p}=&-6Ca_{\textrm{inc}}^{2}J_{1}(\alpha)^{2},\\\ \frac{d^{2}\sigma}{d\Omega
dE^{\prime}}|_{p\to p}=&2Ca_{\textrm{coh}}^{2}J_{0}(\alpha)^{2},\\\
\frac{d^{2}\sigma}{d\Omega dE^{\prime}}|_{p\to
o}=&-\frac{2}{3}Ca_{\textrm{inc}}^{2}J_{1}(\alpha)^{2},\\\ \end{array}$ (15)

where $C=\sqrt{\frac{E_{\textrm{n}}^{\prime}}{\pi
E_{\textrm{n}}E_{\textrm{rec}}k_{B}T}}\exp{\big{[}-\frac{(\Delta
E_{\textrm{n}})^{2}}{4E_{\textrm{rec}}k_{B}T}\big{]}}$ is a factor common to
each, $a_{\textrm{coh}}$ and $a_{\textrm{inc}}$ are the bound coherent and
incoherent nuclear scattering lengths, $J_{l}(\alpha)$ is the spherical Bessel
function of order $l$, $d$ is the equilibrium separation between the atoms in
the molecule, and $\alpha=\frac{1}{2}Kd$. The $3J$ symbol and integral terms
have been evaluated in order to supply concise expressions for the $J=0,1$
transitions, which are the only terms that contribute for cold neutrons and
liquid hydrogen in the vicinity of 15 K.

For the Y-K gas model, the probability of a spin flip interaction for
H${}_{\textrm{para}}\to$H${}_{\textrm{ortho}}$ and
H${}_{\textrm{ortho}}\to$H${}_{\textrm{para}}$ transitions is unity. Cold
neutrons do not have sufficient energy for spin flip interactions from
H${}_{\textrm{para}}\to$H${}_{\textrm{para}}$ scattering events. The
depolarization factor for H${}_{\textrm{ortho}}\to$H${}_{\textrm{ortho}}$
transitions and 5 meV neutrons is determined from the cross section [21] in
equation 15 and is given by

$R=\frac{\sigma_{\textrm{coh}}-\frac{2}{9}\sigma_{\textrm{incoh}}}{\sigma_{\textrm{coh}}+\frac{2}{3}\sigma_{\textrm{incoh}}}=-0.29.$
(16)

The Y-K model only includes hydrogen transitions and does not incorporate
spin-flip scattering from apparatus components.

####  2.1.2 Keinert-Sax kernel model treatment

The basis for hydrogen scattering kernels is the Young and Koppel model, which
accurately describes rotational, vibrational, spin correlations, and free
translations for gaseous hydrogen. However, translational modes in liquid
hydrogen are not free and the Young and Koppel model is not sufficient to
describe the scattering of cold neutrons from liquid hydrogen. The Keinert and
Sax model [23] improves the description of the translational modes compared to
the Young and Koppel model in order to more accurately describe scattering
from liquid hydrogen. The scattering kernels compile for MCNPX from MacFarlane
[24] are based on the hindered translation model from Keinert and Sax and also
incorporate interference between scattered waves from different molecules. The
hydrogen molecules are represented as if there are solid-like clusters of
approximately 20 molecules that diffuse through the liquid. The MacFarlane
model better describes the steep drop in the parahydrogen cross section above
20 meV due to intermolecular spin correlation and the drop around 3-4 meV due
to intramolecular interference [24] when compared to the Young and Koppel gas
model.

Despite the fact that MCNPX does not track spin or determine when a spin
exchange scattering event takes place, the kinematics of a neutron interacting
with orthohydrogen or parahydrogen carries the signature of spin exchange and
this is used to determine whether a spin-flip interaction has taken place. The
first excited state of the hydrogen molecule with $J=1$ is 14.7 meV above the
ground state. Therefore, the signature of the
H${}_{\textrm{ortho}}\to$H${}_{\textrm{para}}$ transition is a substantial
increase in the energy of the neutron with an expected value near 14.7 meV,
while the H${}_{\textrm{para}}\to$H${}_{\textrm{ortho}}$ transition signature
is an energy decrease of the same magnitude.

Figure 2: Parahydrogen scattering kinematics from MCNPX in $\Delta
E_{\textrm{n}}$ and $E_{\textrm{rec}}$ (eq. 13) space (color online). The
color scale indicates probability of an interaction occurring in each
($E_{\textrm{rec}}$, $\Delta E_{\textrm{n}}$) bin. The
H${}_{\textrm{para}}\to$H${}_{\textrm{ortho}}$ to transition is the broad
region in $E_{\textrm{rec}}$ at $\Delta E_{\textrm{n}}=-14.7$ meV

.

Figure 3: Orthohydrogen scattering kinematics from MCNPX in $\Delta
E_{\textrm{n}}$ and $E_{\textrm{rec}}$ (eq. 13) space(color online). The color
scale indicates probability of an interaction occurring in each
($E_{\textrm{rec}}$, $\Delta E_{\textrm{n}}$) bin. The
H${}_{\textrm{ortho}}\to$H${}_{\textrm{para}}$ to transition is the broad
region in $E_{\textrm{rec}}$ at $\Delta E_{\textrm{n}}=14.7$ meV

.

A broad and flat input neutron source that does not match the FNPB energy
spectrum was used to generate a map of the kinematics in $\Delta
E_{\textrm{n}}=E_{\textrm{n}}^{\prime}-E_{\textrm{n}}$ and
$E_{\textrm{rec}}=\frac{K^{2}\hbar^{2}}{4m_{\textrm{n}}}$ space for neutron
scattering from parahydrogen and orthohydrogen. This was done for illustration
and to determine whether the kinematics could be used as a marker for spin-
flip scattering interactions by isolating a region around $\Delta
E_{\textrm{n}}=\pm 15$ meV. Figures 2 and 3, parahydrogen and orthohydrogen
respectively, illustrate two dimensional histograms of the kinematics as
calculated in MCNPX 2.7 using ENDF-VII [15] cross sections.

The regions near $\pm$14.7 meV for each were fit to a 2 dimensional Gaussian
in $\Delta E_{\textrm{n}}$ and $E_{\textrm{rec}}$. A 3$\sigma$ cutoff from
these fits was used in the depolarization calculation such that any scattering
event with $\Delta E_{\textrm{n}}$ and $E_{\textrm{rec}}$ inside 3$\sigma$ is
labeled as a H${}_{\textrm{ortho}}\to$H${}_{\textrm{para}}$ or
H${}_{\textrm{para}}\to$H${}_{\textrm{ortho}}$ and the neutron spin is flipped
with unit probability. If the event is outside of this ellipse, the
interacting species and table 1 are used to determine whether the interaction
is H${}_{\textrm{para}}\to$H${}_{\textrm{para}}$ or
H${}_{\textrm{ortho}}\to$H${}_{\textrm{ortho}}$.

##  3 Monte Carlo calculations of spin-flip scattering

Two separate Monte Carlo models have been developed in order to estimate the
depolarization correction for the NPDGamma experiment. The first model uses
the Y-K gas cross sections in order to estimate the contribution from
orthohydrogen and parahydrogen. For the second, MCNPX 2.7 was modified in
order to provide neutron spin state tracking throughout particle histories
such that the neutron spin could be flipped probabilistically and the signal
contributions from each spin state recorded for each detector.

###  3.1 Gaseous cross-section based simulation

While the Y-K model does not adequately describe the translational modes in
the liquid phase [25], the spin-correlation assumptions in the Y-K model are
valid. The Y-K cross sections in equation 15 are used to directly model the
average polarization of a neutron beam as it passes through a
parahydrogen/orthohydrogen gaseous mixture.

The source neutron beam in this model is initialized with 100% polarization
and is monochromatic with an energy of 5 meV. The mean-free-path of the
neutron is calculated and the neutron is propagated to a new location inside
the hydrogen volume and the event type is chosen probabilistically at that
location according to the value of $R$ in equation 16. If the event is a
capture event, the neutron history is terminated, the average polarization is
recorded, and a new neutron track is generated at the source. If the event is
a scattering event, then the neutron spin is flipped as appropriate, the
kinematics of the collision are determined, and the calculation continues with
a new determination of the mean-free-path and calculation continues until the
neutron track is terminated when it has either escaped the liquid hydrogen
volume or has captured on hydrogen.

###  3.2 MCNPX simulation

MCNPX does not have any internal way of defining and tracking neutron spin.
This capability was implemented by adding a spin state tag for neutrons. As
neutrons propagate in MCNPX, all information regarding the scattering
material, interaction type, outgoing momentum, and possible newly created
particles is determined internally. An additional step was incorporated to
probabilistically change the neutron spin tag from spin up to spin down at
each scattering event. Hydrogen in plastics and rubber components in other
sections of the apparatus is denoted as $H_{\textrm{bound}}$.

Isotope | $P_{\textrm{flip}}$
---|---
H${}_{\textrm{bound}}$ | 0.659(1)
6Li | 0.0002(1)
7Li | 0.43(5)
14N | 0.022(3)
27Al | 0.0033(3)
35Cl | 0.0326(7)
37Cl | 0.00034(4)
55Mn | 0.012(5)
63Cu | 0.00041(9)
65Cu | 0.016(1)
Zn | 0.0011(5)
Table 2: Isotopes with spin flip probabilities determined using Sears cross
sections [26] as implemented in the model along with uncertainties due to
averaging over the neutron energy spectrum. Zinc is applied over all natural
isotopes rather than specific isotopes.

The MCNPX calculation incorporates implicit capture, such that capture events
produce secondary particles and reduce the weight of that neutron and continue
propagating it. The materials in the MCNPX model of the NPDGamma experiment
contain complicated alloys rather than mono-isotopic materials, and MCNPX
handles the bookkeeping for the selection of the isotope involved in a
scattering event. The isotopes relevant to the NPDGamma experiment with non-
zero nuclear spin and an incoherent cross section from Sears [26] are shown in
table 2. In the case of H${}_{\textrm{bound}}$, the spin flip probability
corresponds to the probability of a non-absorption event multiplied by a
factor of $\frac{2}{3}$ due to the lack of separate coherent and incoherent
cross sections in MCNPX. The material definitions in the simulation used
ENDF/B-VII cross sections [15] for the MCNPX calculations.

The source neutron beam begins with 100% polarization and has the energy
spectrum shown in figure 4. This spectrum is representative of the spectrum
determined by the beamline chopper settings [3] without frame overlap. As
neutrons propagate in the model and create $\gamma$-rays, the $\gamma$-rays
inherit the tag indicating the spin of the neutron as the $\gamma$-ray was
created. Additionally, Compton electrons inherit this same tag as they are
created through $\gamma$-ray scattering processes.

Figure 4: Input neutron spectrum used in the MCNPX calculations with
approximate cutoffs for a typical chopped and open beam spectra for the
NPDGamma experiment [3].

.

The average polarization is determined by binning the energy deposition in
each detector by the spin state tag, and a user supplied subroutine performs
these tallies. The geometrically averaged polarization as measured for
detector, $l=0...47$, is given by,

$P^{l}_{\textrm{ave}}=\frac{E^{l}_{\uparrow}-E^{l}_{\downarrow}}{E^{l}_{\uparrow}+E^{l}_{\downarrow}},$
(17)

where $E^{l}_{\uparrow}$ and $E^{l}_{\downarrow}$ are the spin up and spin
down energy deposition tallies respectively. The detector rings, $r$, are
numbered from 1 to 4 with ring 1 being the most upstream ring. Each detector
ring should observe approximately the same average polarization and downstream
rings should observe a slightly lower average polarization because further
downstream detector rings measure signals from neutrons that have penetrated
deeper into the liquid and are therefore more likely to have undergone a spin-
flip scattering event. Additionally, the penetration depth for neutrons is
dependent on the neutron energy, and thus this method of determining the
polarization correction captures both the geometric ring-dependence and energy
dependence of $\Delta^{r}_{\textrm{dep}}(\lambda)$.

##  4 Calculated polarization correction

The MCNPX model provides the versatility to calculate the
$\Delta^{r}_{\textrm{dep}}$ correction for multiple target materials in
addition to hydrogen. In particular, this correction factor was calculated for
auxiliary NPDGamma targets such as the chlorine-35 target and the aluminum
6061 target.

###  4.1 Chlorine-35

The chlorine-35 parity-violating asymmetry is large and well-known, with a
world average of $(2.39\pm 0.136)\times 10^{-5}$ [27][28][29]. A Teflon target
filled with carbon tetrachloride was used as a systematic check on the
apparatus components. Natural chlorine is 75.8% chlorine-35 and 24.2%
chlorine-37[30], and chlorine-35 has a significant spin flip probability (see
table 2). The polarization correction for the chlorine target is $0.9840\pm
0.0003$. Since the chlorine target is thin rather than comparable to the
length of the detector arrays, there is no ring dependence to the polarization
correction for chlorine.

###  4.2 Aluminum-27

The apparatus components, including the hydrogen target vessel and vacuum
windows, are constructed of aluminum 6061 and neutron capture on these
components can carry a parity violating signal that is a dilution factor in
the measurement of the hydrogen asymmetry. The polarization correction for the
apparatus aluminum is calculated using the MCNPX model and is dependent on the
orthohydrogen concentration in the hydrogen vessel. The results are shown in
figure 5 and tabulated for each ring in table 3.

There is a detector ring dependence seen in 5 that is the result of the
geometry of the target vessel and the fact that it is distributed over
$\approx$ 30 cm. Downstream detectors are more likely than upstream detectors
to detect $\gamma$-rays from neutrons that travel further through the liquid
hydrogen volume before being captured. Neutrons that travel further through
the liquid hydrogen are more likely to have scattered before being captured
and therefore are more likely to have experienced a spin-flip scattering event
than neutrons that captured a short distance into the target.

Figure 5: $\Delta^{r}_{\textrm{dep}}$ for neutrons that capture on aluminum 6061 containing cells as a function of parahydrogen concentration and shown for each detector. Ring | $\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,lower}})$ | $\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,upper}})$
---|---|---
1 | 0.9961(3) | 0.9799(8)
2 | 0.9955(3) | 0.9722(8)
3 | 0.9939(5) | 0.9581(13)
4 | 0.9941(6) | 0.9556(19)
Total | 0.9953(2) | 0.9724(5)
Table 3: Polarization correction for the apparatus aluminum at the lower and
upper limit orthohydrogen concentrations for each ring and averaged over the
detector array.

###  4.3 Hydrogen

Similar to the aluminum calculation, the $\Delta^{r}_{\textrm{dep}}$
correction for hydrogen has a detector ring dependence due to being a
distributed volume. $\Delta^{r}_{\textrm{dep}}$ is sensitive to the
orthohydrogen concentration ($x_{\textrm{o}}$) and calculations were performed
at the upper and lower bound orthohydrogen concentrations that were determined
along with the measurement of the parahydrogen scattering cross-section [17],
which are 0.15% and 0.015%, respectively. Since there was no in situ method of
determining the $x_{\textrm{o}}$ for the NPDGamma experiment, the
$x_{\textrm{o}}$ distribution is assumed to be flat between the upper and
lower bounds, with an uncertainty given by
$\frac{1}{\sqrt{12}}(x_{\textrm{o,lower}}-x_{\textrm{o,upper}})$. Similarly,
the polarization correction is given by

$\Delta^{r}_{\textrm{dep,ave}}=\frac{1}{2}\left[\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,lower}})+\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,upper}})\right],$
(18)

with uncertainty equal to

$\sigma^{r}_{\textrm{dep,ave}}=\sqrt{\frac{1}{12}}\left[\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,lower}})-\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,upper}})\right].$
(19)

The neutron energy dependence of $\Delta^{r}_{\textrm{dep}}(\lambda)$ was
calculated using a series of calculations with monochromatic neutron source
terms and varying the energy. The results of this calculation are shown in
figure 6, where the results from the four rings have been averaged together
for each neutron energy. The small deviation from 100% polarization at zero
orthohydrogen concentration is due to the upscattering of neutrons above the
spin-flip threshold as well as the combined kinetic energy of the molecules
and neutron above the spin-flip threshold.

Figure 6: Polarization correction averaged over the detector rings as a
function of initial neutron energy.

Simulations with the MCNPX model were performed at multiple orthohydrogen
concentrations within this range and the results are shown in figure 7. The
polarization upon capture is a linear function of orthohydrogen concentration
within this range. Calculations with monochromatic neutron beams were also
performed and show that the polarization correction decreases with increasing
energy within the typical NPDGamma neutron spectrum range, which is below the
14.7 meV rotational transition energy. The average contribution to
$\Delta^{r}_{\textrm{dep}}$ for hydrogen due to materials other than
orthohydrogen and parahydrogen in table 2 is $2\times 10^{-4}$ and is
independent of orthohydrogen concentration, since the majority of this
contribution is due to spin-flip scattering that takes place upstream of the
target vessel in vacuum windows and other apparatus components.

Figure 7: Polarization correction for each ring as a function of orthohydrogen concentration. Ring | $\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,lower}})$ | $\Delta^{r}_{\textrm{dep}}(x_{\textrm{o,upper}})$
---|---|---
1 | $0.99275(11)$ | $0.94712(28)$
2 | $0.99233(9)$ | $0.94304(24)$
3 | $0.99144(10)$ | $0.93735(25)$
4 | $0.99083(13)$ | $0.93417(33)$
Total | $0.99187(21)$ | $0.94053(56)$
Table 4: Ring dependent polarization correction for liquid hydrogen at the
lower and upper limit orthohydrogen concentrations for each ring and averaged
over the detector array.

The uncertainty due to the unknown orthohydrogen concentration is much larger
than both the uncertainty due to the energy dependence and the ring
dependence, thus the polarization correction is averaged over the energy
spectrum and over all detector rings to produce a single value that applies to
all detectors. The resulting polarization correction from the MCNPX
calculation is between 0.99187 and 0.94053 with an average of $0.966\pm
0.015$. By comparison, the Y-K model uses a monochromatic neutron beam at 5
meV and yields an average polarization between $0.97$ and $0.9959$ with an
average of $0.983\pm 0.007$. In order to be conservative, we take the maximum
and minimum of the results from the Y-K and MCNPX results and assign the
bounds on $\Delta_{\textrm{dep}}$ to be between $0.94053$ and $0.9959$ and
take the probability distribution between these to be flat. Therefore, the
multiplicative polarization correction to the NPDGamma experiment,
$\Delta_{\textrm{dep}}$, is determined to be $0.968\pm 0.016$. This
uncertainty is comparable to that of other multiplicative corrections for the
NPDGamma experiment, namely the geometrical factors calculation [31] and the
measured beam polarization [10] and significantly smaller than the statistical
uncertainty in the measurement.

##  5 Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Basic Energy Sciences under contract number DE-
AC05-00OR22725. This work was also supported by DE-FG02-03ER41258.

## References

  * [1] D. Blyth, J. Fry, N. Fomin, R. Alarcon, L. Alonzi, E. Askanazi, S. Baeßler, S. Balascuta, L. Barrón-Palos, A. Barzilov, J. D. Bowman, N. Birge, J. R. Calarco, T. E. Chupp, V. Cianciolo, C. E. Coppola, C. B. Crawford, K. Craycraft, D. Evans, C. Fieseler, E. Frlež, I. Garishvili, M. T. W. Gericke, R. C. Gillis, K. B. Grammer, G. L. Greene, J. Hall, J. Hamblen, C. Hayes, E. B. Iverson, M. L. Kabir, M. Maldonado-Velázquez, Y. Masuda, J. Mei, R. Milburn, P. E. Mueller, M. Musgrave, H. Nann, I. Novikov, D. Parsons, S. I. Penttila, D. Pocanic, A. Ramirez-Morales, M. Root, A. Salas-Bacci, S. Santra, S. Schröder, E. Scott, P.-N. Seo, E. I. Sharapov, F. Simmons, W. M. Snow, A. Sprow, J. Stewart, E. Tang, Z. Tang, X. Tong, D. J. Turkoglu, R. Whitehead, W. S. Wilburn, First Observation of P-odd $\gamma$ Asymmetry in Polarized Neutron Capture on Hydrogen, Physical Review Letters 13 (24) (2018) 14.
  * [2] M. T. Gericke, R. Alarcon, S. Balascuta, L. Barrón-Palos, C. Blessinger, J. D. Bowman, R. D. Carlini, W. Chen, T. E. Chupp, C. Crawford, S. Covrig, M. Dabaghyan, N. Fomin, S. J. Freedman, T. R. Gentile, R. C. Gillis, G. L. Greene, F. W. Hersman, T. Ino, G. L. Jones, B. Lauss, M. Leuschner, W. R. Lozowski, R. Mahurin, Y. Masuda, J. Mei, G. S. Mitchell, S. Muto, H. Nann, S. a. Page, S. I. Penttilä, W. D. Ramsay, A. Salas-Bacci, S. Santra, M. Sharma, P.-N. Seo, E. I. Sharapov, T. B. Smith, W. M. Snow, W. S. Wilburn, V. Yuan, Measurement of parity-violating gamma-ray asymmetry in the capture of polarized cold neutrons on protons, Physical Review C 83 (1) (2011) 15505\.
  * [3] N. Fomin, G. L. Greene, R. R. Allen, V. Cianciolo, C. Crawford, T. M. Tito, P. R. Huffman, E. B. Iverson, R. Mahurin, W. M. Snow, Fundamental neutron physics beamline at the spallation neutron source at ORNL, Nucl. Instrum. Methods Phys. Res., Sect. A 773 (2015) 45–51.
  * [4] S. Balascuta, R. Alarcon, S. Baeßler, G. Greene, A. Mietke, C. Crawford, R. Milburn, S. Penttila, J. Prince, J. Schädler, The implementation of a super mirror polarizer at the SNS fundamental neutron physics beamline, Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 671 (2012) 137–143.
  * [5] M. T. Gericke, J. D. Bowman, R. D. Carlini, T. E. Chupp, K. P. Coulter, M. Dabaghyan, M. Dawkins, D. Desai, S. J. Freedman, T. R. Gentile, R. C. Gillis, G. L. Greene, F. W. Hersman, T. Ino, G. L. Jones, M. Kandes, B. Lauss, M. Leuschner, W. R. Lozowski, R. Mahurin, M. Mason, Y. Masuda, G. S. Mitchell, S. Muto, H. Nann, S. A. Page, S. I. Penttilä, W. D. Ramsay, S. Santra, P.-N. Seo, E. I. Sharapov, T. B. Smith, W. M. Snow, W. S. Wilburn, V. Yuan, H. Zhu, Upper bounds on parity-violating gamma-ray asymmetries in compound nuclei from polarized cold neutron capture, Physical Review C 74 (6) (2006) 065503.
  * [6] J. Fry, R. Alarcon, R. Allen, E. Askanazi, S. Balascuta, L. Barron-Palos, S. Baeßler, A. Barzilov, C. Blessinger, D. Blyth, J. D. Bowman, J. R. Calarco, T. E. Chupp, C. E. Coppola, C. Crawford, K. Craycraft, M. Dabaghyan, D. Evans, J. Favela, C. Fieseler, N. Fomin, W. Fox, S. Freedman, E. Frlež, C. Fu, C. Garcia, I. Garishvili, M. T. Gericke, R. C. Gillis, K. Grammer, G. L. Greene, J. Hamblen, C. Hayes, F. W. Hersman, T. Ino, E. B. Iverson, G. L. Jones, L. Kabir, S. Kucucker, B. Lauss, Y. Li, R. Mahurin, M. Maldonado-Velazquez, M. McCrea, Y. Masuda, J. Mei, R. Milburn, G. S. Mitchell, P. Mueller, S. Muto, M. Musgrave, H. Nann, I. Novikov, S. Page, D. Parsons, D. Počanić, S. I. Penttilä, W. D. Ramsay, A. Salas-Bacci, S. Santra, P.-N. N. Seo, E. Sharapov, M. Sharma, F. Simmons, T. Smith, W. M. Snow, J. Stuart, E. Tang, Z. Tang, J. Thomison, T. Tong, J. Vanderwerp, S. Waldecker, W. S. Wilburn, W. Xu, V. Yuan, Y. Zhang, Status of the NPDGamma experiment, Hyperfine Interactions 238 (1).
  * [7] P.-N. Seo, L. Barrón-Palos, J. D. Bowman, T. E. Chupp, C. Crawford, M. Dabaghyan, M. Dawkins, S. J. Freedman, T. Gentile, M. T. Gericke, R. C. Gillis, G. L. Greene, F. W. Hersman, G. L. Jones, M. Kandes, S. Lamoreaux, B. Lauss, M. B. Leuschner, R. Mahurin, M. Mason, J. Mei, G. S. Mitchell, H. Nann, S. A. Page, S. I. Penttilä, W. D. Ramsay, A. Salas-Bacci, S. Santra, M. Sharma, T. B. Smith, W. M. Snow, W. S. Wilburn, H. Zhu, High-efficiency resonant rf spin rotator with broad phase space acceptance for pulsed polarized cold neutron beams, Physical Review Special Topics - Accelerators and Beams 11 (8) (2008) 84701.
  * [8] S. Santra, L. Barrn Palos, C. Blessinger, J. D. Bowman, T. E. Chupp, S. Covrig, C. Crawford, M. Dabaghyan, J. Dadras, M. Dawkins, M. T. Gericke, W. Fox, R. C. Gillis, M. B. Leuschner, B. Lozowski, R. Mahurin, M. Mason, J. Mei, H. Nann, S. I. Penttila, A. Salas-Bacci, M. Sharma, W. M. Snow, W. S. Wilburn, A liquid parahydrogen target for the measurement of a parity-violating gamma asymmetry in n→+p→d+$\gamma$, Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 620 (2-3) (2010) 421–436.
  * [9] M. T. Gericke, C. Blessinger, J. D. Bowman, R. C. Gillis, J. Hartfield, T. Ino, M. Leuschner, Y. Masuda, G. S. Mitchell, S. Muto, H. Nann, S. A. Page, S. I. Penttilä, W. D. Ramsay, P.-N. Seo, W. M. Snow, J. Tasson, W. S. Wilburn, A current mode detector array for $\gamma$ ray asymmetry measurements, Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 540 (2-3) (2005) 328–347.
  * [10] M. Musgrave, S. Baeßler, S. Balascuta, L. Barrón-Palos, D. Blyth, J. Bowman, T. Chupp, V. Cianciolo, C. Crawford, K. Craycraft, N. Fomin, J. Fry, M. Gericke, R. Gillis, K. Grammer, G. Greene, J. Hamblen, C. Hayes, P. Huffman, C. Jiang, S. Kucuker, M. McCrea, P. Mueller, S. Penttilä, W. Snow, E. Tang, Z. Tang, X. Tong, W. Wilburn, Measurement of the absolute neutron beam polarization from a supermirror polarizer and the absolute efficiency of a neutron spin rotator for the NPDGamma experiment using a polarized 3 He neutron spin-filter, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 895 (2018) 19–28.
  * [11] J. J. Szymanski, J. D. Bowman, P. P. J. Delheij, C. M. Frankle, J. Knudson, S. Penttilä, S. J. Seestrom, S. H. Yoo, V. W. Yuan, X. Zhu, Ion chamber system for neutron flux measurements, Nucl. Instrum. Methods Phys. Res., Sect. A 340 (3) (1994) 564–571.
  * [12] D. M. Dennison, A Note on the Specific Heat of the Hydrogen Molecule, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 115 (771) (1927) 483–486.
  * [13] L. Barrón-Palos, R. Alarcon, S. Balascuta, C. Blessinger, J. D. Bowman, T. E. Chupp, S. Covrig, C. B. Crawford, M. Dabaghyan, J. Dadras, M. Dawkins, W. Fox, M. T. Gericke, R. C. Gillis, B. Lauss, M. B. Leuschner, B. Lozowski, R. Mahurin, M. Mason, J. Mei, H. Nann, S. I. Penttilä, W. D. Ramsay, A. Salas-Bacci, S. Santra, P. N. Seo, M. Sharma, T. Smith, W. M. Snow, W. S. Wilburn, V. Yuan, Determination of the parahydrogen fraction in a liquid hydrogen target using energy-dependent slow neutron transmission, Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 659 (1) (2011) 579–586.
  * [14] Sigma-Aldridge Corp., 3050 Spruce Street, St. Louis, MO 63103, USA, Product Number 371254-250G.
  * [15] M. Chadwick, M. Herman, P. Obložinský, M. Dunn, Y. Danon, A. Kahler, D. Smith, B. Pritychenko, G. Arbanas, R. Arcilla, R. Brewer, D. Brown, R. Capote, A. Carlson, Y. Cho, H. Derrien, K. Guber, G. Hale, S. Hoblit, S. Holloway, T. Johnson, T. Kawano, B. Kiedrowski, H. Kim, S. Kunieda, N. Larson, L. Leal, J. Lestone, R. Little, E. McCutchan, R. MacFarlane, M. MacInnes, C. Mattoon, R. McKnight, S. Mughabghab, G. Nobre, G. Palmiotti, A. Palumbo, M. Pigni, V. Pronyaev, R. Sayer, A. Sonzogni, N. Summers, P. Talou, I. Thompson, A. Trkov, R. Vogt, S. van der Marck, A. Wallner, M. White, D. Wiarda, P. Young, ENDF/B-VII.1 Nuclear Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields and Decay Data, Nuclear Data Sheets 112 (12) (2011) 2887–2996.
  * [16] S. F. Mughabghab, Atlas of Neutron Resonances, 5th Edition, Elsevier, New York, 2006.
  * [17] K. B. Grammer, R. Alarcon, L. Barron-Palos, D. Blyth, J. D. Bowman, J. Calarco, C. Crawford, K. Craycraft, D. Evans, N. Fomin, J. Fry, M. Gericke, R. C. Gillis, G. L. Greene, J. Hamblen, C. Hayes, S. Kucuker, W. S. Wilburn, Measurement of the scattering cross section of slow neutrons on liquid parahydrogen from neutron transmission, Phys. Rev. B 180301 (18) (2015) 1–6.
  * [18] D. B. Pelowitz, MCNPX User’s Manual, Los Alamos National Laboratory, 2011.
  * [19] J. Young, J. Koppel, Slow Neutron Scattering by Molecular Hydrogen and Deuterium, Phys. Rev. 135 (3A) (1964) A603—-A611.
  * [20] R. Moon, T. Riste, W. Koehler, Polarization Analysis of Thermal-Neutron Scattering, Phys. Rev. 181 (2) (1969) 920–931.
  * [21] M. Zoppi, Neutron scattering of homonuclear diatomic liquids, Physica B: Condensed Matter 183 (3) (1993) 235–246.
  * [22] M. Celli, N. Rhodes, A. K. Soper, M. Zoppi, The total neutron cross section of liquid para-hydrogen, J. Phys.: Condens. Matter 11 (50) (1999) 10229–10242.
  * [23] J. Keinert, J. Sax, Investigation of neutron scattering dynamics in liquid hydrogen and deuterium for cold neutron sources, Kerntechnik 51 (1987) 19.
  * [24] R. E. MacFarlane, Cold-moderator scattering kernel methods, Tech. rep., Los Alamos National Laboratory (LANL), Los Alamos, NM (aug 1998).
  * [25] R. J. Elliott, W. M. Hartmann, Theory of neutron scattering from liquid and solid hydrogen, Proceedings of the Physical Society 90 (3) (1967) 671–680.
  * [26] V. F. Sears, Neutron scattering lengths and cross sections, Neutron News 3 (3) (1992) 26–37.
  * [27] M. Avenier, G. Bagieu, H. Benkoula, J. F. Cavaignac, A. Idrissi, D. H. Koang, B. Vignon, R. Wilson, Parity non-conservation in the radiative capture of polarized neutrons by35Cl, Nuclear Physics, Section A 436 (1) (1985) 83–92.
  * [28] V. Vesna, É. Kolomenski, V. Lobashev, V. Nazarenko, A. Pirozhkov, L. Smotritski, Y. Sobolev, N. Titov, Observation of parity nonconservation in the integrated $\gamma$ spectrum from (n, $\gamma$) reactions in Cl, Br, Cd, Sn, and La nuclei, Soviet Journal of Experimental and Theoretical Physics Letters 36 (1982) 209.
  * [29] G. S. Mitchell, C. S. Blessinger, J. D. Bowman, T. E. Chupp, K. P. Coulter, M. Gericke, G. L. Jones, M. B. Leuschner, H. Nann, S. A. Page, S. I. Penttilä, T. B. Smith, W. M. Snow, W. S. Wilburn, A measurement of parity-violating gamma-ray asymmetries in polarized cold neutron capture on 35Cl, 113Cd, and 139La, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 521 (2-3) (2004) 468–479.
  * [30] J. Meija, T. B. Coplen, M. Berglund, W. A. Brand, P. De Bièvre, M. Gröning, N. E. Holden, J. Irrgeher, R. D. Loss, T. Walczyk, T. Prohaska, Isotopic compositions of the elements 2013 (IUPAC Technical Report), Pure and Applied Chemistry 88 (3) (2016) 293–306.
  * [31] K. Grammer, D. Blyth, J. Bowman, N. Fomin, G. Greene, M. Musgrave, E. Tang, Z. Tang, Monte Carlo calculation and verification of the geometrical factors for the NPDGamma experiment, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 903 (2018) 302–308.