# Modeling Electrical Resistance Drift with Ultrafast Saturation of OTS
Selectors

Yiğit Demirağ
Department of Electrical and Electronics Engineering
École Polytechnique Fédérale de Lausanne (EPFL)
Lausanne, Switzerland
<EMAIL_ADDRESS>
&Ekmel Özbay
Department of Electrical and Electronics Engineering
Department of Physics
Bilkent University
Ankara, Turkey
<EMAIL_ADDRESS>
&Yusuf Leblebici
Department of Electrical and Electronics Engineering
École Polytechnique Fédérale de Lausanne (EPFL)
Lausanne, Switzerland
<EMAIL_ADDRESS>

###### Abstract

Crossbar array architecture is an essential design element for densely
connected Non-Volatile Memory (NVM) applications. To overcome intrinsic sneak
current problem of crossbar arrays, each memory unit is serially attached to a
selector unit with highly nonlinear current-voltage (I-V) characteristics.
Recently, Ovonic Threshold Switching (OTS) materials are preferred as
selectors due to their fabrication compatibility with PRAM, MRAM or ReRAM
technologies; however, OTS selectors suffer from the temporal drift of its
threshold voltage. First, based on Poole-Frenkel conduction, we present time
and temperature dependent model that predicts temporally evolving I-V
characteristics, including threshold voltage of OTS selectors. Second, we
report an ultrafast saturation (at $\sim 10^{3}$ seconds) of the drift and
extend the model to predict the time of drift saturation. Our model shows
excellent agreement with OTS devices fabricated with 8 nm technology node at
25 °C and 85 °C ambient temperatures. The proposed model plays a significant
role in understanding OTS device internals and the development of reliable
threshold voltage jump table.

_K_ eywords OTS selectors $\cdot$ Drift problem $\cdot$ Device modeling
$\cdot$ Neuromorphic hardware

##  1 Introduction

An OTS material is a thin-film two terminal amorphous chalcogenide, whom
electrical conductivity can rapidly change from the high resistive state (HRS)
to low resistive state (LRS) by applying large potential exceeding a specific
threshold voltage ($V_{th}$). The conductivity difference between HRS and LRS
can be as high as 4 orders of magnitude; nevertheless, the device immediately
switches back when the applied potential is cut [1]. Its high ON/OFF current
ratio and fast switching make OTS a promising candidate material for selector
applications.

OTS is a chalcogenide material; hence it can crystallize. However, once
crystallization starts, it is not feasible within device operation range to
initiate melting and recover amorphous state. Therefore, crystallized OTS
selector units are always assumed non-operational. As a solution, OTS
materials can be carefully selected to have lower ionicity and higher
hybridization, which lead to more directed covalent bonds to significantly
slow crystallization process [2].

OTS selectors also perfectly match to physical and electrical scaling
properties of NVM technologies. Owing to thin film compatibility with mature
metallization techniques and CMOS stack support, densely connected 3D crossbar
arrays with OTS selectors have been demonstrated [3].

##  2 The Drift Problem of OTS Selectors

The major problem of OTS technology is that the electrical conductivity of the
selector decreases over time (Fig. 1(a)), called the drift problem. We observe
that this conductivity decrease is not consistent, but saturates in time.
Whether the conductance is drifting or already saturated, application of any
potential higher than the threshold voltage ($V>V_{th}$), resets the drift and
revert to selector’s initial HRS level. Saturation of the drift is a rarely
reported physical phenomenon in the literature, but critically important for
understanding and developing OTS technology. In our OTS devices, we observe an
ultrafast drift saturation (Fig. 1(b)) which takes places at least 2-3 orders
of magnitude faster than reported drift measurements [4, 5].

The main problem is due to the increase of $V_{th}$ as conductance drifts.
READ and WRITE operations require a known $V_{th}$ level of OTS. If $V_{th}$
increases and applied READ/WRITE pulses could not pass threshold value, then
the selector device stays in HRS, hence READ/WRITE operations fail. One
practical solution is to determine new $V_{th}$ with applying various prior
READ voltages and detecting the threshold voltage value, which certain current
level is reached [6]. However, this solution requires additional support
circuitry and increases the power consumption of the device; therefore, it is
unfavorable. On the other hand, physical modeling of time-dependent resistance
drift may lead to efficient solutions.

Figure 1: (a) Experimental (dot) and simulated (line) resistivity measurements
of the same OTS device with 800 nm2 contact area at different ambient
temperatures. The READ voltage values are 2.48 V and 2.64 V respectively for
25 °C and 85 °C experiments. The resistivity difference is due to thermally
activated Poole-Frenkel behavior. (b) Resistance drift saturation
measurements. The saturation of drift is faster at 85 °C due to faster
annihilation of the defects.

Modeling the drift behavior of OTS selectors is of vital importance for two
reasons. First, a validated model can be useful for developing reliable time
and temperature dependent jump-table of $V_{th}$. Second, a physically
grounded model can provide a comprehensive understanding of temporally
evolving non-measurable material properties such as activation energy
$(E_{a})$ and inter-trap distance $(\Delta z)$ to investigate the fabricated
material in more detail. Although there exist models capturing the drift
behavior, these methods are either validated only for phase change memory
(PCM) on a short range of time and ambient temperature or unable to predict
the saturation of drift due to employing simple power-law-like models [7, 8] .

Figure 2: Temporally evolving activation energy $(E_{a})$ and inter-trap
distance $(\Delta z)$ calculated by the drift model. Upper and lower plots are
generated for OTS selectors whom time and temperature dependent resistivity
measurements are shown at Fig. 1

##  3 Modeling OTS Drift and Its Saturation

###  3.1 Time and Temperature Dependent Drift Model

The physical phenomena governing the resistance drift on amorphous
chalcogenide materials is yet to be fully understood. Raty et al. [9], Gabardi
et al. [10] and Zipoli et al. [11] have recently provided a significant
insight into the microscopic picture of the the drift mechanism. Using ab-
initio simulations, it was found that there exist energetically unstable
homopolar bonds and defects in melt-quenched amorphous. As these unstable
defects naturally transform into lower-energy structures with time, the
distance between intrinsic traps increase. The structure evolves into a more
crystalline-like state, however without the necessary long-range order (Fig.
3).

Figure 3: From left to right: Homopolar bonds in amorphous OTS are
energetically unstable hence naturally disappear with time. This mechanism is
called structural relaxation (SR) and affects activation energy ($E_{a}$) and
inter-trap distance ($\Delta z$). Change of mainly these two causes the drift
of device resistivity and threshold voltage as the electrical transport model
of OTS materials follow the Poole-Frenkel conduction.

To model the kinetics of the structural relaxation, we started with a drift
model developed by Le Gallo et al., which has been previously validated on
GeTe and GST [5]. In this model, the bond network state of chalcogenide is
denoted with an order parameter $\Sigma$. $\Sigma$ is a normalized parameter
between 0 (low-ordered highly stressed state) and 1 (ideal, energetically
favorable relaxed state). Whenever $V>V_{th}$ is applied to OTS selector,
amorphous network state resets and stressed with initial distance
$\Sigma(0)=\Sigma_{0}$ from the equilibrium. As network collectively relaxed
through more favorable states with time, the energy barrier required to
overcome, $E_{b}$, monotonically increases and it is assumed to be linearly
dependent to $\Sigma$:

$E_{b}(t)=E_{s}(1-\Sigma(t)),$ (1)

where $E_{s}$ is the final energy barrier to reach the most relaxed state at
$\Sigma=0$. With an Arrhenius type temperature-dependence, this relaxation
occurs at the rate of $r(t)=v_{0}\exp\left(-E_{b}(t)/k_{B}T\right)$, where
$v_{0}$ is an attempt-to-relax frequency, $k_{B}$ is the Boltzmann constant.
After that the evolution of $E_{b}(t)$ can be calculated by:

$\frac{d\Sigma(t)}{dt}=-v_{0}\Delta_{\Sigma}\exp\left(-\frac{E_{b}}{k_{B}T(t)}\right).$
(2)

At a constant temperature, Eq. 2 can be solved analytically to track the
progress of structural relaxation, such that

$\Sigma(t)=-\frac{k_{B}T}{E_{s}}\log\left(\frac{t+\tau_{0}}{\tau_{1}}\right),$
(3)

where
$\tau_{1}=\left(k_{B}T/\nu_{0}\Delta_{\Sigma}\mathrm{E}_{\mathrm{s}}\right)\exp\left(E_{s}/k_{B}T\right)$
and $\tau_{0}=\tau_{1}\exp\left(-\Sigma_{0}E_{s}/k_{B}T\right)$. Once $\Sigma$
is calculated by Eq. 3, an empirical linear relationship between structural
relaxation, activation energy and inter-trap distance can be written as:

$\begin{array}[]{c}{E_{a0}(t)=E^{*}-\alpha\Sigma(t),}\\\ {\Delta
z(t)=s_{0}/\Sigma(t),}\end{array}$ (4)

where $E^{*}$ is the activation energy at the equilibrium, $\alpha$ and
$s_{0}$ are the constants linking change in $\Sigma$ to material properties.
Finally, the temperature dependence of activation energy is assumed to follow
the Varshni effect $\left(E_{a}=E_{a0}-\xi T^{2}\right)$, as the optical
bandgap of the material depends on the temperature [12].

###  3.2 Subthreshold Electrical Transport Model

Amorphous chalcogenide materials are known to follow Poole-Frenkel
subthreshold conductivity behavior [13]. The Poole-Frenkel effect suggests
that thermal excitation and strong electric field release trapped carriers
from ionizable defect centers, which are believed to create Coulomb potential.
In this work, we used a previously developed field and temperature dependent
3D Poole-Frenkel emission model with field independent mobility [14]. We first
calculated the potential profile between defect centers located at $r=0$ and
$r=\Delta z$ in all directions using Eq. 5, where $\beta$ is the Poole-Frenkel
constant, $e$ the electronic charge, $\theta$ the direction of escape relative
to applied E-field $F$.

$V(r,\theta,F)=-eFr\cos(\theta)-\frac{\beta^{2}}{4e}\left(\frac{1}{r}+\frac{1}{\Delta
z-r}\right)+\frac{\beta^{2}}{e\Delta z}$ (5)

Figure 4: The potential profile between OTS defect centers are calculated via
Eq. 5. The amount of energy barrier may increase (left) or decrease (right)
significantly due to the angle $(\theta)$ between strong electric field
$3\times 10^{8}$ V/m and the direction of escape.

The potential profile between Coulombic defect centers separated by $\Delta z$
is shown for OTS material in Fig. 4. The energy barrier lowering due to the
Poole-Frenkel effect then can be calculated by

$E_{\mathrm{PF}}(F,\theta)=-\max_{r}V(r,\theta,F).$ (6)

Finally, assuming Boltzmann statistics, we calculated the subthreshold
electrical conductivity of the selector with:

$\sigma=e\mu\frac{K}{4\pi}\int_{0}^{\pi}\exp\left(-\frac{E_{a}-E_{PF}(F,\theta)}{k_{B}T}\right)2\pi\sin(\theta)d\theta.$
(7)

###  3.3 Predicting Saturation Time of Electrical Resistance Drift

In our experiments, we observed an unusually ultrafast saturation of the
electrical resistance drift. Moreover, the drift saturation point changes as a
function of the ambient temperature.

In the strong form of the drift model proposed by Le Gallo et al. [15], the
evolution of subthreshold electrical resistance can be predicted; however, it
falls short predicting a drift saturation point. To extend the previous model
to predict the saturation time, we hypothesize that identical devices at
different ambient temperatures that saturate at different times, eventually
converge to the same $E_{a}$ and $\Delta z$ at the time of the saturation.
This hypothesis requires $\Sigma(t)$ to be the same and constant for identical
devices at different temperatures after the saturation time, $t_{SAT}$:

$\frac{d\Sigma(t)}{dt}=-r(t)\Delta_{\Sigma}=0\hskip
28.45274pt\mathrm{for}~{}~{}t>t_{SAT}.$ (8)

During the training of the model, structural relaxation parameters are tuned
according to this constraint given at Eq. 8. $t_{SAT}$ values used in the
training are experimentally gathered as the saturation times of $V_{th}$ (see
Fig. 1 (b)). Figure 5 shows the evolution of $\Sigma(t)$ for a trained model.

Figure 5: Calculated structural relaxation constant $\Sigma$ given for two
ambient temperatures 25 °C and 85 °C.

###  3.4 Predicting Time Evolution of Threshold Voltage

As the electrical resistance of the selector drifts with time the threshold
voltage also drifts. To predict the time evolution of $V_{th}$, a mere time
and temperature dependent subthreshold electrical model would not suffice. The
model requires an extension to explain the moment of threshold switching for
OTS.

To account the sudden increase of conductivity during threshold switching, we
combined Poole-Frenkel subthreshold transportation with Okuto-Crowell impact
ionization. Okuto-Crowell impact ionization is an empirical model which is
based on electron-avalanche multiplication effect due to the high electric
field ($\sim 4\times 10^{7}$ V/m) inside the OTS material [16]. With this
extension illustrated in Fig. 6 (a), we demonstrated a successful prediction
of time-evolution of $V_{th}$ of OTS selectors at different ambient
temperatures (Fig. 6 (b)).

Figure 6: (a) Illustration of Poole-Frenkel (blue) conduction model and Okuto-
Crowell impact ionization model (red). Together two models are sufficient for
describing threshold switching of OTS. (b) Prediction of $V_{th}$ for T=25 °C,
85 °C and $t=300$ s.

##  4 Physically-Realistic Parameter Optimization

To adjust model parameters according to experimental measurements, the model
I/O is matched with experimental conditions. The implemented model takes the
same control inputs with the fabricated device (voltage, ambient temperature)
and returns the same measurable quantity (resistivity). Figure 7 shows that
fabricated OTS device can be modeled as a black box whom physical
characteristics are represented by a set of parameters $\theta$.

Figure 7: For OTS device, the only control parameters are ambient temperature
and applied voltage.

Figure 8 shows the proposed drift saturation model with 17 parameters. To
successfully optimize these model parameters to match the fabricated device,
we consider two requirements. First, in an ideal situation, the proposed model
and the fabricated device must output the same resistivity level when applied
the same voltage and ambient temperature. Therefore, the aim is to minimize
the difference between measured resistivity of the fabricated device,
$R(\theta)$, and the modeled resistivity, $R(\hat{\theta})$, by tuning model
parameters, $\theta$.

Figure 8: Model parameters used in this work. The temporal nature of
electrical conductance is due to structural relaxation model, which affects
$\Delta z$ and $E_{A}$ of OTS. Poole-Frenkel model and impact ionization
combined, define the electrical behavior of the OTS and enable tracking of
time-dependent $V_{th}$.

Second, the model parameters with physical correspondence must stay within
their physically-realistic ranges. To limit every parameter with different
upper (UB) and lower bounds (LB) as in Eq. 9, several optimization methods
could be used e.g., simulated annealing, evolutionary or gradient-based search
algorithms. We utilized simulated annealing for its easy implementation and
despite its computationally-heavy search, it successfully minimized the loss
function inside of Eq. 9 with physically-realistic parameters [17].

$\displaystyle
x=\arg\min_{\theta}\frac{1}{2}\sum\left(R(t)-\hat{R}_{\theta}(t)\right)^{2},\text{
subject to }LB_{\text{i}}<\theta_{i}<UB_{i}.$ (9)

##  5 Conclusion

We reported an ultrafast saturation phenomenon (at $\sim 10^{3}$ seconds) of
resistance drift on OTS materials, which are promising selector candidates in
the next generation NVM (PRAM, MRAM and ReRAM) crossbar technologies. An
electrical transport model is proposed to describe time and temperature
dependent OTS I-V characteristics. The model based on structural relaxation,
Poole-Frenkel conduction, and impact ionization, is shown to be in close
agreement with our devices fabricated with 8 nm node technology and tested at
25 °C and 85 °C ambient temperatures for $\sim 10^{4}$ seconds. The models and
physical parameters (including $E_{a}$ and $\Delta z$) provide valuable
insight into non-measurable material properties. With the support of drift
saturation and $V_{th}$ prediction, our model may play a significant role in
the development of reliable $V_{th}$ jump tables.

## References

  * [1] Stanford R Ovshinsky.  Reversible Electrical Switching Phenomena in Disordered Structures.  21(20):6.
  * [2] A. Velea, K. Opsomer, W. Devulder, J. Dumortier, J. Fan, C. Detavernier, M. Jurczak, and B. Govoreanu.  Te-based chalcogenide materials for selector applications.  7(1).
  * [3] DerChang Kau, Stephen Tang, Ilya V. Karpov, Rick Dodge, Brett Klehn, Johannes A. Kalb, Jonathan Strand, Aleshandre Diaz, Nelson Leung, Jack Wu, Sean Lee, Tim Langtry, Kuo-wei Chang, Christina Papagianni, Jinwook Lee, Jeremy Hirst, Swetha Erra, Eddie Flores, Nick Righos, Hernan Castro, and Gianpaolo Spadini.  A stackable cross point Phase Change Memory.  pages 1–4. IEEE.
  * [4] A. Pirovano, A.L. Lacaita, A. Benvenuti, F. Pellizzer, and R. Bez.  Electronic Switching in Phase-Change Memories.  51(3):452–459.
  * [5] Manuel Le Gallo, Abu Sebastian, Daniel Krebs, Milos Stanisavljevic, and Evangelos Eleftheriou.  The complete time/temperature dependence of I-V drift in PCM devices.  pages MY–1–1–MY–1–6. IEEE.
  * [6] A Sebastian, M Le Gallo, W W Koelmans, N Papandreou, H Pozidis, and E Eleftheriou.  Multi-level storage in phase-change memory devices.  page 6.
  * [7] SangBum Kim, Byoungil Lee, Mehdi Asheghi, Fred Hurkx, John P. Reifenberg, Kenneth E. Goodson, and H.-S. Philip Wong.  Resistance and Threshold Switching Voltage Drift Behavior in Phase-Change Memory and Their Temperature Dependence at Microsecond Time Scales Studied Using a Micro-Thermal Stage.  58(3):584–592.
  * [8] Nicola Ciocchini, Marco Cassinerio, Davide Fugazza, and Daniele Ielmini.  Modeling of Threshold-Voltage Drift in Phase-Change Memory (PCM) Devices.  59(11):3084–3090.
  * [9] Jean Yves Raty, Wei Zhang, Jennifer Luckas, Chao Chen, Riccardo Mazzarello, Christophe Bichara, and Matthias Wuttig.  Aging mechanisms in amorphous phase-change materials.  6(1).
  * [10] S. Gabardi, S. Caravati, G. C. Sosso, J. Behler, and M. Bernasconi.  Microscopic origin of resistance drift in the amorphous state of the phase-change compound GeTe.  92(5).
  * [11] Federico Zipoli, Daniel Krebs, and Alessandro Curioni.  Structural origin of resistance drift in amorphous GeTe.  93(11).
  * [12] Y P VARSHNI.  TEMPERATURE DEPENDENCE OF THE ENERGY GAP IN SEMICONDUCTORS.  page 6.
  * [13] Robert M. Hill.  Hopping conduction in amorphous solids.  24(192):1307–1325.
  * [14] Manuel Le Gallo, Matthias Kaes, Abu Sebastian, and Daniel Krebs.  Subthreshold electrical transport in amorphous phase-change materials.  17(9):093035.
  * [15] Manuel Le Gallo, Daniel Krebs, Federico Zipoli, Martin Salinga, and Abu Sebastian.  Collective Structural Relaxation in Phase-Change Memory Devices.  page 1700627.
  * [16] Y. Okuto and C. R. Crowell.  Threshold energy effect on avalanche breakdown voltage in semiconductor junctions.  18(2):161–168.
  * [17] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller.  Equation of State Calculations by Fast Computing Machines.  21(6):1087–1092.


# Magnetization Plateau of the Distorted Diamond Spin Chain

Yuki Ueno1, Tomosuke Zenda1, Yuta Tachibana1, Kiyomi Okamoto1 and Tôru
Sakai1,2 1Graduate School of Material Science, University of Hyogo, Hyogo
678-1297, Japan
2National Institutes for Quantum and Radiological Science and Technology
(QST), SPring-8, Hyogo 679-5148, Japan

(Received September 1, 2019)

###### Abstract

The frustrated quantum spin system on the distorted diamond chain lattice is
investigated using the numerical diagonalization of finite-size clusters and
the level spectroscopy analysis. In the previous work this system was revealed
to exhibit the 1/3 magnetization plateau due to two different mechanisms
depending on the coupling parameters, and the phase diagram at the 1/3
magnetization was obtained. In the present work it is found that the 1/3
magnetization plateau vanishes for sufficiently large $XY$-like coupling
anisotropy. The phase diagram based on the level spectroscopy analysis is also
presented.

##  I Introduction

The magnetization plateau is one of interesting phenomena in the field of the
strongly correlated electron systems. It possibly appears as the quantization
of magnetization, when the one-dimensional quantum spin system satisfies the
necessary condition

$\displaystyle S-m={\rm integer},$ (1)

where $S$ is the total spin and $m$ is the magnetization per unit celloshikawa
. The S=1/2 distorted diamond spin chainokamoto0  is a strongly frustrated
quantum spin system which exhibits the 1/3 magnetization plateau according to
the condition. This system was proposed as a good theoretical model of the
compound Cu3(CO3)2(OH)2, called azuritekikuchi . Actually the magnetization
measurement of the azurite detected a clear magnetization plateau at 1/3 of
the saturation magnetization. The theoretical study by the numerical exact
diagonalizationokamoto1  suggested that the 1/3 magnetization plateau is
induced by two different mechanisms. One is based on the ferromagnetic
mechanism and the other is due to the formation of singlet dimers and free
spins. In the case of the azurite the 1/3 plateau is supposed to be due to the
latter mechanism. In some previous theoretical worksokamoto2 ; sakai1 ; sakai2
; ito  on quantum spin systems it was reported that the magnetization plateau
disappears in the presence of sufficiently large $XY$-like (easy-plane)
coupling anisotropy. Thus in this paper we introduce the $XY$-like anisotropy
to the $S=1/2$ distorted diamond spin chain and consider the stability of the
1/3 magnetization plateau against the anisotropy in the case of both
mechanisms. We note that the anisotropy inversion phenomena was found in the
distorted diamond chain with the $XXZ$ anisotropyoka-ichi ; tokuno . For this
purpose we use the numerical exact diagonalization of finite-size clusters and
the level spectroscopy analysis to find the quantum phase transition between
the plateau and the no-plateau phases.

##  II Model

We investigate the model described by the Hamiltonian

$\displaystyle{\cal H}={\cal H}_{0}+{\cal H}_{\rm Z}$ (2) $\displaystyle{\cal
H}_{0}=J_{1}\sum_{j=1}^{N/3}\left[(\mbox{\boldmath$S$}_{3j-1}\cdot\mbox{\boldmath$S$}_{3j})_{\lambda}+(\mbox{\boldmath$S$}_{3j}\cdot\mbox{\boldmath$S$}_{3j+1})_{\lambda}\right]+J_{2}\sum_{j=1}^{N/3}(\mbox{\boldmath$S$}_{3j+1}\cdot\mbox{\boldmath$S$}_{3j+2})_{\lambda}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+J_{3}\sum_{j=1}^{N/3}\left[(\mbox{\boldmath$S$}_{3j-2}\cdot\mbox{\boldmath$S$}_{3j})_{\lambda}+(\mbox{\boldmath$S$}_{3j}\cdot\mbox{\boldmath$S$}_{3j+2})_{\lambda}\right]$
(3) $\displaystyle{\cal H}_{\rm Z}=-H\sum_{l=1}^{N}S_{l}^{z}$ (4)
$\displaystyle(\mbox{\boldmath$S$}_{l}\cdot\mbox{\boldmath$S$}_{m})_{\lambda}\equiv
S_{l}^{x}S_{m}^{x}+S_{l}^{y}S_{m}^{y}+\lambda S_{l}^{z}S_{m}^{z}$ (5)

where $\mbox{\boldmath$S$}_{j}$ is the spin-1/2 operator, $J_{1}$, $J_{2}$,
$J_{3}$ are the coupling constants of the exchange interactions, and $\lambda$
is the coupling anisotropy. The schematic picture of the model is shown in
Fig. 1. In this paper we consider the case of the $XY$-like (easy-plane)
anisotropy, namely $\lambda<1$.

Figure 1: The model of the $S=1/2$ distorted diamond spin chain.

$N$ is the number of spins and $L$ is defined as the number of the unit cells,
namely $N=3L$. For $L$-unit systems, the lowest energy of ${\cal H}_{0}$ in
the subspace where $\sum_{j}S_{j}^{z}=M$, is denoted as $E(L,M)$. The reduced
magnetization $m$ is defined as $m=M/M_{\rm s}$, where $M_{\rm s}$ denotes the
saturation of the magnetization, namely $M_{\rm s}=3L/2$ for this system.
$E(L,M)$ is calculated by the Lanczos algorithm under the periodic boundary
condition (${\bf S}_{N+1}={\bf S}_{1}$) and the twisted boundary condition
($S^{x,y}_{N+1}=-S^{x.y}_{1},S^{z}_{N+1}=S^{z}_{1}$), for $L=$4, 6 and 8.
Under the twisted boundary condition we calculate the lowest energy $E_{{\rm
TBC}P=+}(L,M)$ ($E_{{\rm TBC}P=-}(L,M)$) in the subspace where the parity is
even (odd) with respect to the lattice inversion at the twisted bond.

##  III Magnetization plateau

In the isotropic coupling case $\lambda=1$ the previous theoretical study
indicates that the model (2) exhibits the 1/3 magnetization plateau due to two
different mechanisms depending on the parameters $J_{1}$, $J_{2}$, $J_{3}$
okamoto2 . Particularly $J_{2}$ plays a key role. The plateau is based on the
ferrimagnet-like mechanism for smaller $J_{2}$, while it is based on the
structure where the singlet dimer lies on the $J_{2}$ bond and the other spins
are free for larger $J_{2}$. The small $J_{2}$ plateau and the large $J_{2}$
one are denoted as the plateaux A and B, respectively. The second-order
quantum phase transition occurs on the boundary between the two plateau
phases. In this paper the quantum phase transition from the plateau phase to
the no-plateau one with respect to the $XY$-like coupling anisotropy is
investigated in the both cases of the plateaux A and B.

##  IV Level spectroscopy analysis and phase diagram

In order to detect the quantum phase transitions among the plateaux-A, -B and
no-plateau phases, the level spectroscopy analysiskitazawa1 ; kitazawa2  is
one of the best methods. According to this analysis, we should compare the
following three energy gaps;

$\displaystyle\Delta_{2}={E(L,M-2)+E(L,M+2)-2E(L,M)\over 2},$ (6)
$\displaystyle\Delta_{{\rm TBC}P=+}=E_{{\rm TBC}P=+}(L,M)-E(L,M),$ (7)
$\displaystyle\Delta_{{\rm TBC}P=-}=E_{{\rm TBC}P=-}(L,M)-E(L,M).$ (8)

The level spectroscopy method indicates that the smallest gap among these
three gaps for $M=L=M_{\rm s}/3$ determines the phase at $m=1/3$.
$\Delta_{2}$, $\Delta_{{\rm TBC}P=+}$ and $\Delta_{{\rm TBC}P=-}$ correspond
to the no-plateau, plateau-A and plateau-B phases, respectively. Especially,
$\Delta_{{\rm TBC}P=+}$ and $\Delta_{{\rm TBC}P=-}$ directly reflect the
symmetry of two plateau states. The $\lambda$ dependence of these three gaps
for $J_{1}=1.0$, $J_{2}=1.1$ and $J_{3}=0.1$ is shown in Fig. 2 for $L=$4, 6
and 8.

Figure 2: $\lambda$ dependence of the three gaps $\Delta_{2}$(solid line),
$\Delta_{{\rm TBC}P=+}$(dotted line) and $\Delta_{{\rm TBC}P=1}$(dashed line)
for $J_{1}=1.0$, $J_{2}=1.1$ and $J_{3}=0.1$. Black, blue and red lines
correspond to the system sizes $L=$4, 6 and 8, respectively.

Assuming that the finite-size correction is proportional to $1/L^{2}$, we
estimate the phase boundaries in the thermodynamic limit from every level-
cross point. The phase diagram in the $J_{2}$-$\lambda$ plane at $m=1/3$ for
$J_{1}=1.0$ and $J_{3}=0.1$ is shown in Fig. 3.

Figure 3:  Phase diagram at $m=1/3$ for $J_{1}=1.0$ and $J_{3}=0.1$. Crosses,
circles and triangles correspond to the boundaries between the plateau-A and
no-plateau, between the plateau-B and no-plateau, and between the plateaux-A
and -B phases, respectively.

It suggests that both plateaux are so stable against the $XY$-like anisotropy
that the plateaux vanish at quite large negative $\lambda$.

##  V Concluding Remarks

Using the numerical exact diagonalization and the level spectroscopy analysis,
the $S=1/2$ distorted diamond spin chain is investigated. It is found that the
1/3 magnetization plateau vanishes for sufficiently large negative $XY$-like
anisotropy, in both cases of the ferrimagnet-like and the dimer-monomer
plateaux. A typical phase diagram at $m=1/3$ is presented in Fig.3.

One of the remarkable natures of the plateau phase diagram Fig.3 is the
survival of the 1/3 plateau to the $\lambda<0$ region, where the direct
$S^{z}-S^{z}$ interaction is ferromagnetic. The essential mechanism of the
plateau-A is the Lieb-Mattis ferrimagnetismLM . In the simple $S=1/2$ $XXZ$
chain with $\lambda<0$, although the direct $S^{z}-S^{z}$ interaction is
ferromagnetic, the antiferromagnetic $S^{z}-S^{z}$ correlation survives to
$\lambda=-1$ induced by the antiferromagnetic $S^{x}-S^{x}$ and $S^{y}-S^{y}$
correlations. Thus, the plateau A survives to near $\lambda=-1$, where the
ferromagnetic long-range order begins. On the other hand, the plateau-B state
is attributed to the singlet-dimer formation at the $J_{2}$ bonds, which makes
the nodal spins nearly free. This type of singlet dimer survives to
$\lambda=-1/\sqrt{2}=-0.707$ in the $S=1/2$ bond-alternating $XXZ$
chainkohmoto ; oka-sugi . Further, if we neglect $J_{3}$, the present system
becomes a trimerized $S=1/2$ $XXZ$ chain, in which Okamoto and Kitazawaoka-
kita  showed that the no-plateau state appears for $\lambda<-0.729$. These
facts well explains why the boundary of the platea-B state and no-plateau
state lies around $\lambda=0.70\sim 0.75$. Thus both plateaus survive into the
$\lambda<0$ region, and the plateau-A is more stable than the plateau-B as the
$\lambda=-1$ line is approached.

We think that the no-plateau state is essentially the Tomonaga-Luttinger
liquid state in the magnetic field. We can see a sudden decrease of the phase
boundary between the plateau-A and the no-plateau $J_{2}=0.2$. The physical
explanation for this phenomena is a future problem. We believe that our
plateau phase diagram is important for the full understanding of the distorted
diamond chain and will yields important informations if the distorted diamond
chain with the $XXZ$ anisotropy is found or synthesized in future.

## Acknowledgment

This work was partly supported by JSPS KAKENHI, Grant Numbers 16K05419,
16H01080 (J-Physics) and 18H04330 (J-Physics). A part of the computations was
performed using facilities of the Supercomputer Center, Institute for Solid
State Physics, University of Tokyo, and the Computer Room, Yukawa Institute
for Theoretical Physics, Kyoto University.

## References

  * (1) M. Oshikawa, M. Yamanaka and I. Affleck, Phys. Rev. Lett. 78, 1984 (1997).
  * (2) K. Okamoto T. Tonegawa, Y. Takahashi, and M. Kaburagi, J. Phys.: Condens. Matter 11, 10485 (1999).
  * (3) H. Kikuchi, Y. Fujii, M. Chiba, S. Mitsudo, T. Idehara, T. Tonegawa, K. Okamoto, T. Sakai, T. Juwai, and H. Ohta, Phys. Rev. Lett. 94, 227201 (2005).
  * (4) K. Okamoto, T. Tonegawa, and M. Kaburagi, J. Phys.: Condens. Matter, 15, 5979 (2003).
  * (5) K. Okamoto and A. Kitazawa, J. Phys. A: Math. Gen. 32, 4601 (1999).
  * (6) T. Sakai and S. Yamamoto, Phys. Rev. B 60, 4053 (1999).
  * (7) S. Yamamoto and T. Sakai, Phys. Rev. B 62, 3795 (2000).
  * (8) K. Ito, S. Yokoo, K. Okamoto and T. Sakai, Polyhedron 136, 42 (2017).
  * (9) K. Okamoto and Y. Ichikawa, J. Phys. Chem. Solids 63, 1575 (2002).
  * (10) A. Tokuno and K. Okamoto J. Phys. Soc. Jpn. 74, Suppl. 157 (2005)
  * (11) A. Kitazawa, J. Phys. A: Math. Gen. 30, L285 (1997).
  * (12) K. Nomura and A. Kitazawa, J. Phys. A: Math. Gen. 31, 7341 (1998).
  * (13) E. Lieb and D. Mattis, J. Math. Phys. 3, 749 (1962).
  * (14) M. Kohmoto, M. den Nijs and L. P. Kadanoff, Phys. Rev B 24, 5229 (1981).
  * (15) K. Okamoto and T. Sugiyama, J. Phys. Soc. Jpn. 57, 1610 (1988).
  * (16) K. Okamoto and A. Kitazawa, J. Phys.: Math. Gen. 32, 4601 (1999).


# Comparing Bayes factors and hierarchical inference for testing
general relativity with gravitational waves

Maximiliano Isi<EMAIL_ADDRESS>Center for Computational
Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010  Will M.
Farr<EMAIL_ADDRESS>Center for Computational Astrophysics, Flatiron
Institute, 162 5th Ave, New York, NY 10010  Department of Physics and
Astronomy, Stony Brook University, Stony Brook NY 11794, USA  Katerina
Chatziioannou<EMAIL_ADDRESS>Department of Physics, California
Institute of Technology, Pasadena, California 91125, USA  LIGO Laboratory,
California Institute of Technology, Pasadena, California 91125, USA

###### Abstract

In the context of testing general relativity with gravitational waves,
constraints obtained with multiple events are typically combined either
through a hierarchical formalism or though a combined multiplicative Bayes
factor. We show that the well-known dependence of Bayes factors on the
analysis priors in regions of the parameter space without likelihood support
can lead to strong confidence in favor of incorrect conclusions when one
employs the multiplicative Bayes factor. Bayes factors ${\cal{O}}(1)$ are
ambivalent as they depend sensitively on the analysis priors, which are rarely
set in a principled way; additionally, combined Bayes factors
$>{\cal{O}}(10^{3})$ can be obtained in favor of the incorrect conclusion
depending on the analysis priors when many $\mathcal{O}(1)$ Bayes factors are
multiplied, and specifically when the priors are much wider than the
underlying population. The hierarchical analysis that instead infers the
ensemble distribution of the individual beyond-general-relativity constraints
does not suffer from this problem, and generically converges to favor the
correct conclusion. Rather than a naive multiplication, a more reliable Bayes
factor can be computed from the hierarchical analysis. We present a number of
toy models showing that the practice of multiplying Bayes Factors can lead to
incorrect conclusions.

##  I Introduction

With an increasing number of LIGO-Virgo Aasi  _et al._ (2015); Acernese  _et
al._ (2015) gravitational wave (GW) observations, we can leverage the
collective set of measurements to study the properties of the astrophysical
objects that generate GWs Abbott  _et al._ (2021a, b) and the validity of
general relativity (GR) Isi  _et al._ (2019); Abbott  _et al._ (2021c, d). Two
broad and complementary approaches exist for drawing inferences from sets of
detections. The first relies on the posterior distribution for some model and
its continuous parameters whose range of possible values encapsulates the
different physics we would like to study. The second phrases a question of
interest in the language of model selection between two discrete hypotheses to
compute Bayes factors (BFs). The latter approach is common in the context of
testing GR where one introduces a parametrized deviation of the signal as
predicted by GR Li  _et al._ (2012a, b); Agathos  _et al._ (2014); Sampson
_et al._ (2013, 2014); Meidam  _et al._ (2014), but further examples relate to
higher-order modes of the radiation Chatziioannou  _et al._ (2019), GW memory
Hübner  _et al._ (2020), the neutron star equation of state Del Pozzo  _et
al._ (2013); Chatziioannou  _et al._ (2015); Abbott  _et al._ (2020); Ghosh
_et al._ (2021); Isi  _et al._ (2017); Callister  _et al._ (2017),
gravitational lensing Liu  _et al._ (2021); Lo and Magaña Hernandez (2021),
the association between GWs and potential electromagnetic counterparts Ashton
_et al._ (2020), GW ringdowns Ota and Chirenti (2021), and the signal
detection problem in general Veitch and Vecchio (2008, 2010); Cornish and
Littenberg (2015); Isi  _et al._ (2018); Ashton  _et al._ (2019); Pratten and
Vecchio (2021); Cornish  _et al._ (2021).

Although posteriors and BFs are mathematically related, in practice approaches
focusing on one or the other can come to seemingly different, and even
contradictory, conclusions. For example, Ref. Bustillo  _et al._ (2021)
considers the case of testing the no-hair theorem with GW ringdowns and shows
that BFs can favor the incorrect conclusion even in cases where the posterior
has minimal support for it. The alternative of working directly with the
posterior and hierachical inference has been introduced in the context of
tests of GR Isi  _et al._ (2019); Abbott  _et al._ (2021c, d) after the
limitations of using BFs to combine information from multiple events were
pointed out in Zimmerman  _et al._ (2019). Specifically, multipyling BFs
corresponds to assuming that a GR deviation will manifest independently and
distributed according to the underlying prior in each observation Zimmerman
_et al._ (2019). In this paper, we further this line of argument to highlight
issues with the use of BFs in the context of nested models and show that the
hierarchical modeling of population distributions offers a more flexible and
reliable alternative.

BFs, or marginalized-likelihood ratios, provide a succinct way to compare the
likelihood of two models in light of some data. The BF comparing a hypothesis
$\mathcal{H}_{0}$ to another $\mathcal{H}_{1}$ for some data $d$ is given by

$\mathcal{B}^{0}_{1}\equiv\frac{P(d\mid\mathcal{H}_{0})}{P(d\mid\mathcal{H}_{1})}=\frac{\int
p(d\mid\theta_{0},\mathcal{H}_{0})\,p(\theta_{0}|\mathcal{H}_{0})\,\mathrm{d}\theta_{0}}{\int
p(d\mid\theta_{1},\mathcal{H}_{1})\,p(\theta_{1}|\mathcal{H}_{1})\,\mathrm{d}\theta_{1}}\,,$
(1)

where the likelihoods are marginalized over some (potentially different) sets
of parameters $\theta_{0/1}$ for the $\mathcal{H}_{0/1}$ hypotheses
respectively. The definition of the hypotheses encompasses the choice of
parameter priors $p(\theta\mid\mathcal{H})$, as well as any other assumptions
built into the functional form of the likelihoods
$p(d\mid\theta,\mathcal{H})$. A larger value of $\mathcal{B}^{0}_{1}$ (or,
equivalently, its natural logarithm, $\ln\mathcal{B}^{0}_{1}$) indicates a
preference for $\mathcal{H}_{0}$ over $\mathcal{H}_{1}$. When enhanced with
prior weights for each hypothesis, $P(\mathcal{H}_{0/1})$, this returns the
betting odds in favor of one model over the other, conditional on the observed
data—namely,

$\mathcal{O}^{0}_{1}\equiv\frac{P(\mathcal{H}_{0})P(d\mid\mathcal{H}_{0})}{P(\mathcal{H}_{1})P(d\mid\mathcal{H}_{1})}=\frac{P(\mathcal{H}_{0})}{P(\mathcal{H}_{1})}\mathcal{B}^{0}_{1}\,.$
(2)

To avoid expressing an _a priori_ preference for either model, it is common to
set $P(\mathcal{H}_{0})=P(\mathcal{H}_{1})$, and so
$\mathcal{O}^{0}_{1}=\mathcal{B}^{0}_{1}$.

BFs reduce the complicated problem of selecting between two models of reality
to a single number—a feature which lies at the core of its appeal, but also of
its shortcomings. Like other scalar statistics, interpreting this one number
is usually far from straightforward, leading to somewhat _ad hoc_ scales such
as Kass and Raftery (1995). This difficulty is worsened by the fact that BFs
necessarily are affected by _all_ aspects of a model, including those
corresponding to less interesting regions of the parameter space like regions
of the prior space for which the likelihood offers no support. Consequently,
BFs can vary wildly with different choices of prior bounds, which are often
set arbitrarily. Since priors can rarely be set from first principles,
calibrating BFs tends to require large scale injection campaigns Isi  _et al._
(2018)—although this is only possible when the injections themselves can be
designed in a principled way (i.e., when we actually know how to simulate
expected astrophysical distributions). Even when the model is specified
correctly and we can take BFs at face value, the result offers no insight as
to why exactly one model is to be preferred.

All these drawbacks compound when one attempts to combine multiple
observations by multiplying BFs computed with a fixed prior, which enhances
the sensitivity to prior choices. Moreover, naive BF computations from
collections of events impose strong, generally unrealistic assumptions
Zimmerman  _et al._ (2019). Multiplying BFs obtained from individual events
results in a collective BF that assumes the targeted effect (say, a deviation
from GR) manifests independently for each observation. On the other hand,
generating a collective BF from the product of likelihoods from multiple
measurements presumes that the effect manifests identically for all
observations. Neither of these are valid assumptions in general Zimmerman  _et
al._ (2019). In general the degree to which a targeted effect appears
independently versus identically in each observation is something that should
be leared from the data; this insight leads directly to hierarchical modeling
James and Stein (1961); Lindley and Smith (1972); Efron and Morris (1977);
Rubin (1981).

Rather than assuming a fixed and known distribution (e.g., all events are the
same, or all the events are different), a hierarchical analysis works by
inferring the underlying distribution of the parameter whose values encode the
targeted hypotheses. For example, in the context of tests of GR, a parameter
$x$ may represent the magnitude of a GR violation, so that the GR (non-GR)
hypothesis implies $x=0$ ($x\neq 0$); likewise, when searching for GW memory,
this could be the amplitude of the effect in question. Once this parameter is
identified, we can use hierarchical inference to characterize its values
across the observed events—using the collection of measurements holistically
to infer the distribution of $x$. The challenge here lies in choosing a
suitable parametrization for this underlying distribution.

We can circumvent this through a moment expansion: as a first approximation,
we will only be interested in recovering the mean and standard deviation of
the unknown distribution, so we can parametrize it as a Gaussian, whose mean
and standard deviation ($\mu$, $\sigma_{\mathrm{pop}}$) we are to measure from
the data Isi  _et al._ (2019); the null hypothesis will often be constructed
so that it is recovered for $\mu=\sigma_{\mathrm{pop}}=0$. This approach
allows us to study the population of measurements without assuming the
targeted effect manifests either identically or distinctly for all events,
learning more about the population distribution, which now plays the role of a
prior, with each new observation.

In this paper, we compare the BF and hierarchical approaches directly in the
context of multiple GW observations, and argue that BFs are unreliable in any
context in which the prior does not adapt to the observations at hand. We show
that in such context, BFs do not have the right scaling with the number of
events, even in simple situations, due to their inherently strong dependence
on the (fixed) priors in regions of no likelihood support. We show that
hierarchical posteriors do not suffer from such limitations as they rely on
“priors” that are inferred from the data, and have a weaker dependence on
hyperparameter assumptions.

We begin in Sec. II by reviewing some basic properties of BFs obtained from
single observations and their interplay with Occam penalties. Then, in Sec.
III, we study the scaling of BFs when combining multiple observations under
two example distributions for the true parameters under consideration; we show
that this approach can often lead to the incorrect conclusion. We also
consider the same scenarios under the hierarchical approach, showing that it
does not suffer from the same drawbacks. In Sec. IV, we consider a final model
that cannot be obtained as a special case of the distribution assumed by the
hierarchical analysis, showing that the hierarchical approach yields the
correct result even then. We conclude in Sec. V.

##  II Single event: Occam penalty vs goodness-of-fit

Figure 1: Bayes factor computation as a Savage-Dickey density ratio. The BF
($\mathcal{B}$) comparing a model in which some parameter $x$ takes the value
$x=0$ to one in which $x\neq 0$ is given by the ratio of the posterior to the
prior evaluated at $x=0$. The main panel shows the case of a normal posterior
truncated by a broad prior (half-width $\Delta=10$, in blue) and a narrow
prior ($\Delta=1$, in orange). The inset shows the scaling of $\mathcal{B}$ as
a function of $\Delta$ (solid curve); for broad priors,
$\mathcal{B}\propto\Delta$ (dashed line). In this toy example, we assume the
posterior peaks always at $x=0$.

Before tackling the problem of multiple observations, we first review the
behavior of BFs computed from a single event. A typical situation that is
simple to understand analytically is that of nested models, i.e., two models
constructed so that one can be recovered as a special case of the other.
Parametrized tests of GR, for example, typically involve nested hypotheses
wherein the non-GR model is characterized by all the usual GR signal
parameters plus one or more additional variables $x_{\rm nGR}$ that quantify
the deviation from GR Blanchet and Sathyaprakash (1995, 1994); Arun  _et al._
(2006a, b); Yunes and Pretorius (2009); Mishra  _et al._ (2010); Li  _et al._
(2012a, b); Agathos  _et al._ (2014). These parametrizations are usually
constructed so that GR is recovered when the deviation parameters vanish.
Then, to determine whether GR is favored, the data are analyzed with some
broad (typically flat) prior on the deviations in order to compute BFs
comparing $x_{\rm nGR}=0$ to $x_{\rm nGR}\neq 0$.

In that spirit, consider some real-valued parameter $x$ and two related
hypotheses $\mathcal{H}_{x=0}$ and $\mathcal{H}_{x\neq 0}$, respectively
defined to imply $x=0$ and $x\neq 0$, with some prior over $x$ and any other
relevant parameters. Since $x=0$ is a special case of the model in which $x$
is allowed to vary over a broad range including the origin, we say the
hypotheses are nested.111Here we are identifying $\mathcal{H}_{x\neq 0}$ with
the model in which $x$ is allowed to vary freely; this is legitimate because
the point $x=0$ is a set of measure zero, so it does not need to be explicitly
excised from the arbitrary-$x$ model. In that case, the BF comparing the two
is given exactly by the Savage-Dickey ratio Dickey (1971); Verdinelli and
Wasserman (1995),

$\mathcal{B}^{x=0}_{x\neq 0}=\frac{p(x=0\mid d)}{p(x=0)}\,,$ (3)

where $p(x=0\mid d)$ is the marginal posterior and $p(x=0)$ is the prior, both
evaluated at $x=0$. In other words, the BF in favor of $x=0$ is simply the
ratio of the posterior to the prior evaluated at the origin.

We can elucidate the role of the prior in Eq. (3) by considering a specific
functional form. For simplicity, assume the marginal posterior is given by a
standard normal distribution, truncated symmetrically around the origin by
some uniform prior of half-width $\Delta$ (i.e., flat in $-\Delta<x<\Delta$).
Then, the BF can be computed analytically to yield

$\mathcal{B}^{x=0}_{x\neq
0}=\frac{1}{\sqrt{2\pi}}\frac{2\Delta}{\Phi(\Delta)-\Phi(-\Delta)}=\sqrt{\frac{2}{\pi}}\frac{\Delta}{\mathrm{erf}(\Delta/\sqrt{2})}\,,$
(4)

in terms of the standard cumulative distribution function
$\Phi(x)\equiv\left(1+\mathrm{erf}(x/\sqrt{2})\right)/2$ and the error
function
$\mathrm{erf}(x)\equiv(2/\sqrt{\pi})\int_{0}^{x}\exp(-y^{2})\,\mathrm{d}y$.
Since $\mathrm{erf}(\Delta/\sqrt{2})\to 1$ for large $\Delta$, the BF can be
made to favor $x=0$ with arbitrarily-high confidence by sufficiently
broadening the prior—in fact, $\mathcal{B}^{x=0}_{x\neq 0}\propto\Delta$ in
the large $\Delta$ limit. We illustrate this in Fig. 1.

The dependence on the prior range is a general feature not specific to our
example: the same data can produce arbitrary odds in favor of a specific value
of a parameter ($x=0$ here) relative to a model with increasing prior volume
(proxied by $\Delta$). This is related to the concept of the Occam penalty in
Bayesian inference: BFs do not only favor the model that fits the data best,
but also the one that is simplest—where simplicity is defined as a model’s
ability to fit the data without having to significantly constrain its
parameters relative to their a priori allowed values. The interplay between
goodness of fit and Occam penalty creates the possibility of a BF that
strongly favors the incorrect conclusion. In the context of testing GR, if the
theory is indeed violated, then some observation can give $p(x=0\mid d)\ll 1$
in Eq. (3); yet, this can always be countered with a wide enough prior that
makes the Occam penalty $1/p(x=0)$ so large that goodness of fit cannot
overcome it. This is the expected manifestation of the Occam penalty in BFs;
the failure is symptomatic of a mismatch between the implemented prior and the
observer’s expectation.

Since we generally have little a priori information about the nature of
possible beyond-GR effects, and these effects are not strongly constrained by
the likelihood of an individual measurement, a natural inclination is to make
the prior much wider than the likelihood. However, the argument above shows
that a broad prior prevents us from detecting small deviations from GR, which
is an important regime for tests of GR Perkins and Yunes (2022)—either because
GR is close to correct or because of selection effects that disfavor the
detection of signals with morphology far from GR. The same tension arises in
other contexts where BFs are used without a principled prior. In the next
section, we show that this behavior is not unique to single-event analyses but
carries over to combined constraints.

##  III Combining events under a broad prior

We now turn to collections of measurements and show that the “combined”
multiplicative BF does not have the correct scaling in the regime of interest,
with support accumulating in favor of the null hypothesis even when this
conclusion is incorrect. Combining BFs from multiple events will only lead to
the right conclusion when the deviation (e.g., the deviation from GR) is large
enough to be apparent in individual posteriors, negating the need for
combining observations in the first place. We then show that the hierarchical
approach is not susceptible to this issue.

Again, consider a single parameter, ${x}$, that stands in for the magnitude of
the GR deviation or any other effect of interest, and let ${x}=0$ be our null
hypothesis (e.g., GR is correct). We conduct experiments to measure ${x}$, and
assume additive Gaussian noise so that the observed value $x_{\mathrm{obs}}$
is normally distributed about the true value ${x}$ with standard deviation
(measurement noise) $\sigma_{\mathrm{obs}}$:

$p\left(x_{\mathrm{obs}}\mid{x}\right)=\frac{1}{\sqrt{2\pi\sigma_{\mathrm{obs}}^{2}}}\exp\left(-\frac{\left(x_{\mathrm{obs}}-{x}\right)^{2}}{2\sigma_{\mathrm{obs}}^{2}}\right).$
(5)

In the previous section, we considered a simplified case of this model, in
which we had a single measurement with $x_{\mathrm{obs}}=0$ and
$\sigma_{\mathrm{obs}}=1$.

We consider two “populations” for the true values of ${x}$ under repeated
measurements, i.e., the true magnitude of the GR deviation for each observed
GW event. In the first the true value of ${x}$ is fixed to some value
${x}_{0}$ for each measurement so that the population distribution is a Dirac
delta,

$p_{1}\left({x}\right)=\delta\left({x}-{x}_{0}\right).$ (6)

In the second, the value of ${x}$ is randomly distributed with mean zero and
standard deviation $\sigma_{0}$ for each measurement, so that the population
distribution is

$p_{2}\left({x}\right)=\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}\exp\left(-\frac{{x}^{2}}{2\sigma_{0}^{2}}\right).$
(7)

In the language of tests of gravity, the two models recover GR (our null
hypothesis) when ${x}_{0}=0$ or $\sigma_{0}=0$, respectively. In the second
model, the _mean_ deviation from GR vanishes, but the actual variation in a
given measurement fluctuates. In both cases the deviation parameter ${x}$ is
assumed independent and identically distributed (iid) for each measurement, so
there is some prior choice such that multiplying BFs for repeated measurements
is optimal Zimmerman  _et al._ (2019).222The optimal prior is the actual
population from which the true parameters controlling the measurements are iid
draws. However, that prior is unknown in realistic situations because the true
distribution of GR deviations is not known a priori. Rather, in all cases we
assume a measurement is analyzed with a flat prior on $-\Delta<{x}<\Delta$,
following common practice.

###  III.1 Bayes factors

With that prior, the BF between the null and generalized hypotheses (for
concreteness, “GR” and “non-GR”) for a given observation $x_{\mathrm{obs}}$ is
a generalization of Eq. (4)

$\mathcal{B}^{\mathrm{GR}}_{\mathrm{nGR}}=\frac{2\Delta}{\sqrt{2\pi\sigma_{\mathrm{obs}}^{2}}\left(\Phi\left(\frac{\Delta-
x_{\mathrm{obs}}}{\sigma_{\mathrm{obs}}}\right)-\Phi\left(\frac{-\Delta-
x_{\mathrm{obs}}}{\sigma_{\mathrm{obs}}}\right)\right)}\\\
\times\exp\left(-\frac{x_{\mathrm{obs}}^{2}}{2\sigma_{\mathrm{obs}}^{2}}\right).$
(8)

When the prior is very broad, $\Delta\gg x_{\mathrm{obs}}$, and the Gaussian
posterior on ${x}$ is not meaningfully truncated by the prior, then the BF
simplifies to

$\mathcal{B}^{\mathrm{GR}}_{\mathrm{nGR}}\simeq\frac{2\Delta}{\sqrt{2\pi\sigma_{\mathrm{obs}}^{2}}}\exp\left(-\frac{x_{\mathrm{obs}}^{2}}{2\sigma_{\mathrm{obs}}^{2}}\right),$
(9)

corresponding to the $\mathcal{B}^{x=0}_{x\neq 0}\propto\Delta$ limit
discussed in the previous section. Ensuring $\Delta\gg x_{\mathrm{obs}}$ is a
common analysis choice because this prior permits the true value of ${x}$ to
correspond to the observed value $x_{\mathrm{obs}}$ which is the value of
${x}$ that maximizes the likelihood in each observation. A broad prior is also
desired when combining multiple observations in order to accommodate the
expected scatter in the individual likelihoods.

If we choose to combine observations by adding log BFs333In this paper, all
logarithms are natural logarithms. (equivalently, multiplying BFs), it is
sufficient to compute the expected value of an individual log BF under the
true deviations; the expected total log BF will then be the expected
individual log BF times the number of events. The expected value of the log of
Eq. (8) is not expressible in closed form, but can be straightforwardly
computed numerically. In the limit that $\Delta\gg x_{\mathrm{obs}}$, the
expected log BFs under our two populations of GR deviations become

$\displaystyle\left.\left\langle\ln\mathcal{B}^{\mathrm{GR}}_{\mathrm{nGR}}\right\rangle\right|_{{x}_{0}}$
$\displaystyle\simeq\ln\frac{2\Delta}{\sqrt{2\pi\sigma_{\mathrm{obs}}^{2}}}-\frac{\sigma_{\mathrm{obs}}^{2}+{x}_{0}^{2}}{2\sigma_{\mathrm{obs}}^{2}}$
$\displaystyle\simeq\ln\frac{\Delta}{\sigma_{\mathrm{obs}}}-0.23-\frac{\sigma_{\mathrm{obs}}^{2}+{x}_{0}^{2}}{2\sigma_{\mathrm{obs}}^{2}},$
(10)

and

$\displaystyle\left.\left\langle\ln\mathcal{B}^{\mathrm{GR}}_{\mathrm{nGR}}\right\rangle\right|_{\sigma_{0}}$
$\displaystyle\simeq\ln\frac{2\Delta}{\sqrt{2\pi\sigma_{\mathrm{obs}}^{2}}}-\frac{\sigma_{\mathrm{obs}}^{2}+\sigma_{0}^{2}}{2\sigma_{\mathrm{obs}}^{2}}$
$\displaystyle\simeq\ln\frac{\Delta}{\sigma_{\mathrm{obs}}}-0.23-\frac{\sigma_{\mathrm{obs}}^{2}+\sigma_{0}^{2}}{2\sigma_{\mathrm{obs}}^{2}}.$
(11)

From the above, we can see that, whenever the GR deviation is nonzero but
small enough to be undetectable in a single observation
($0<{x}_{0},\sigma_{0}\ll\sigma_{\mathrm{obs}}$), then for $\Delta\gtrsim
2\sigma_{\mathrm{obs}}$ the expected log BF is positive, and evidence
accumulates _in favor_ of GR even though there is a deviation. For deviations
that are marginally detectable in a single observation
${x}_{0},\sigma_{0}\simeq\sigma_{\mathrm{obs}}$, choosing a wide,
uninformative prior $\Delta\gtrsim 3.5\sigma_{\mathrm{obs}}$ (which, recall,
is necessary to ensure that all observations $x_{\mathrm{obs}}$ in a modest-
sized catalog of observations are within the prior range) will result in
evidence that accumulates _against_ modifications to GR.

An exact calculation of the expected log BF without the assumption that
$\Delta\gg x_{\mathrm{obs}}$ appears in Fig. 2. Interestingly, any value of
$\Delta$ will accumulate evidence for GR when ${x}_{0}=0$ or $\sigma_{0}=0$;
but if ${x}_{0}>0$ or $\sigma_{0}>0$, so that the GR model is incorrect, prior
choices that encompass deviation parameter values comparable to those
observed, i.e., $\Delta\simeq\mathrm{few}\times\sigma_{\mathrm{obs}}$,
accumulate evidence for the incorrect GR model unless the deviation parameter
is comparable to or larger than the observational uncertainty. In this regime,
either the BF fails to select the correct theory (non-GR) on average (with
more and more certainty as the number of observations grows); or the deviation
is so large that it is marginally detectable (i.e., “${\sim}1\sigma$”) with a
single observation. Multipliying BFs in this situation is counter-productive,
and the best constraint is achieved with a single measurement.

Figure 2:  The expected log BF in favor of GR from the toy model discussed in
Sec. III (color) as a function of the prior width $\Delta$ (ordinate) and the
true deviation parameter ${x}_{0}$ (abscissa left) or the scatter in the zero-
mean deviation parameter $\sigma_{0}$ (abscissa right). Both quantities are
shown relative to the measurement uncertainty on the deviation parameter
$\sigma_{\mathrm{obs}}$. For any non-zero value of the deviation parameter or
its scatter, there is a prior width that results in incorrectly accumulating
evidence _for_ GR (above dashed line); moreover, for reasonable prior widths
of a few times the measurement uncertainty (so that observed values of the
deviation parameter lie within the prior range and the prior is therefore
broad and uninformative) the deviation parameter must be large enough to be
marginally detectable in a single measurement
(${x}_{0},\sigma_{0}\sim\sigma_{\mathrm{obs}}$) before evidence against the
incorrect GR model accumulates on average over many stacked measurements. Fig.
3 demonstrates this effect for $\Delta=5\sigma_{\mathrm{obs}}$.

Figure 3:  The expected log BF for the GR model when
$\Delta=5\sigma_{\mathrm{obs}}$ (a broad, uninformative prior) versus the
fixed value of the deviation parameter (${x}_{0}$, left) or scatter in the
zero-mean deviation parameters ($\sigma_{0}$, right). For this choice of prior
width, the average log BF favors the incorrect GR model until the deviation
parameter is ${x}_{0},\sigma_{0}\gtrsim 1.3\sigma_{\mathrm{obs}}$, marginally
detectable in a single measurement. Shading shows the expected variation in
the log BF ($\pm 1\sigma$).

Even when the log BF is expected to take the correct sign on average, the
result will vary for any given set of detections. In the regime where
$\Delta\gg x_{\mathrm{obs}}$, we can quantify this through the variance
associated with the means in Eqs. (III.1) and (III.1). For the former, this is
just

$\left.\mathrm{var}\left(\ln\mathcal{B}^{\rm GR}_{\rm
nGR}\right)\right|_{{x}_{0}}\simeq\left(\frac{{x}_{0}}{\sigma_{\mathrm{obs}}}\right)^{2}+\frac{1}{2}\,,$
(12)

while for the latter we have

$\left.\mathrm{var}\left(\ln\mathcal{B}^{\rm GR}_{\rm
nGR}\right)\right|_{\sigma_{0}}\simeq\frac{\left(\sigma_{0}^{2}+\sigma_{\mathrm{obs}}^{2}\right)^{2}}{2\sigma_{\mathrm{obs}}^{4}}\,.$
(13)

These scalings are illustrated in Fig. 3 for $\Delta=5\sigma_{\mathrm{obs}}$.
On the left hand side, the variance asymptotes to $1/2$ for decreasing
$\sigma_{\mathrm{obs}}$ and fixed intrinsic scatter $\sigma_{0}$; on the right
hand side, the BF variance keeps increasing as we make the individual
measurements more precise (decreasing $\sigma_{\mathrm{obs}}$) because for
narrower posteriors the value at $x=0$ varies more drastically from event to
event. In the case where GR is correct (${x}_{0}=\sigma_{0}=0$), the scatter
in the single-event log BF is large enough that it is not extremely uncommon
for moderately large catalogs to yield evidence for the wrong hypothesis, even
when the null hypothesis is correct (Fig. 4).

Figure 4: Log BF for a growing catalog of $N$ observations in the case where
the null hypothesis (GR) is correct, and where
$\Delta=2\sigma_{\mathrm{obs}}$, for 200 simulations (represented by each thin
gray trace). On average, the log BF favors the right hypothesis (blue line),
but there is enough variance in this quantity that some of the individual
catalogs favor the wrong hypothesis, even for a moderately high number of
detections (traces below the dashed line).

The qualitative behavior here can be understood in the context of Zimmerman
_et al._ (2019). The flat prior on $x$ implicitly assumes that each
observation has a true $x$ parameter that is independent of other observations
and uniformly distributed on $\left(-\Delta,\Delta\right)$. If one chooses
$\Delta$ large enough to include all event likelihoods and not truncate them,
then the flat prior for the true parameter implicitly demands that most of the
true deviation parameters are comparable to $\Delta$. If, instead, the
deviation parameters are considerably smaller than the observational
uncertainty $\sigma_{\mathrm{obs}}$—which is the regime where stacking
multiple events _should_ be the most beneficial—then the assumption is so
badly violated that pooling observations prefers the incorrect model where the
deviation parameters are zero to the even-less-correct model where the
deviation parameters are iid uniform on $\left(-\Delta,\Delta\right)$.

The solution to this problem is to allow the assumed population to adapt its
properties to the stacked set of observations, for example by allowing its
location and scale to fit the set, as suggested in Isi  _et al._ (2019). This
effectively constructs a model that better represents our beliefs about the
combined data set.

###  III.2 Hierarchical treatment

We now revisit the two experiments above using a hierarchical model for the
distribution of GR deviations, instead of computing BFs with a uniform prior.
Assuming we are only interested in the first two moments of the distribution,
we parametrize the true deviations as drawn from a Gaussian such that
${x}\sim\mathcal{N}(\mu,\sigma_{\mathrm{pop}})$ Isi  _et al._ (2019). The goal
will be to infer the $\mu$ and $\sigma_{\mathrm{pop}}$ hyperparameters from
the collection of measurements, and to quantify agreement with the null
hypothesis $\mu=\sigma_{\mathrm{pop}}=0$ based on the corresponding 2D
posterior. The two non-GR models we considered above are encompassed within
this parametrization when $(\mu={x}_{0},\,\sigma_{\mathrm{pop}}=0)$ and
$(\mu=0,\,\sigma_{\mathrm{pop}}=\sigma_{0})$.444The next section examines a
non-GR model that is not fully encompassed in the Gaussian hierarchical model.

As before, we assume that the observed value $x_{\mathrm{obs}}$ in an
individual event is normally distributed around the true value per Eq. (5),
and we take the true values themselves to be distributed normally given $\mu$
and $\sigma_{\mathrm{pop}}$. In that case, the likelihood for a given
observation becomes (see Appendix A)

$p(x_{\mathrm{obs}}\mid\mu,\sigma_{\mathrm{pop}},\sigma_{\mathrm{obs}})=\frac{1}{\sqrt{2\pi\sigma_{\mathrm{tot}}^{2}}}\exp\left(-\frac{\left(x_{\mathrm{obs}}-\mu\right)^{2}}{2\sigma_{\mathrm{tot}}^{2}}\right),$
(14)

where
$\sigma_{\mathrm{tot}}^{2}\equiv\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{pop}}^{2}$
is the total variance arising from the combination of statistical uncertainty
and the intrinsic population scatter. Unlike in Eq. (5), the true value ${x}$
for the individual measurement does not appear in this likelihood because we
have marginalized over it and replaced it with the hyperparameters $\mu$ and
$\sigma_{\mathrm{pop}}$.

Figure 5:  The GR quantile $\mathcal{Q}_{\mathrm{GR}}$ (color) obtained in a
hierarchical analysis of data in which there is a fixed deviation for all
events (${x}_{0}$, ordinate left) or a scatter across events ($\sigma_{0}$,
ordinate right), as a function of the number of events ($N$, abscissa). The
reported value of $\mathcal{Q}_{\mathrm{GR}}$ is the median over 50 simulated
catalogs with $N$-events. Low values of $\mathcal{Q}_{\mathrm{GR}}$ indicate
that the null hypothesis ($\mu=\sigma_{\mathrm{pop}}=0$) is disfavored; for
reference, the dashed line marks the point at which GR is disfavored at
$1\sigma$, i.e., $\mathcal{Q}_{\mathrm{GR}}=\exp(-1/2)\approx 0.61$. Unlike in
Fig. 2, given enough $N$ we always detect the deviation, even when
${x}_{0},\sigma_{0}<\sigma_{\mathrm{obs}}$.

Consider the same two true distributions for $x$ above, which depart from GR
as given by Eqs. (6) and (7). We simulate catalogs of detections following
these distributions and analyze them hierarchically with the likelihood of Eq.
(14), and Gaussian priors on $\mu$ and $\sigma_{\mathrm{pop}}$ (with zero mean
and standard deviation equal to $\sigma_{\mathrm{obs}}$, restricting to
positive values for $\sigma_{\mathrm{pop}}$). Although other choices are
possible, it is natural to tie the scale of the hyperprior to the measurement
uncertainty because this is the only scale built into the problem a priori.
Furthermore, this choice is guaranteed to be sufficiently broad in the small-
deviation regime (${x}_{0},\sigma_{0}<\sigma_{\mathrm{obs}}$) in which we are
interested, and smooth enough to accommodate larger deviations if needed.

We quantify agreement with the null hypothesis through the marginal posteriors
for $\mu$ and $\sigma_{\mathrm{pop}}$, as well as the credible level at which
GR is recovered in the 2D posterior for those two quantities (the 2D quantile,
$\mathcal{Q}_{\mathrm{GR}}$ in Abbott  _et al._ (2021c)), defined as

$\mathcal{Q}_{\mathrm{GR}}\equiv\int_{p<p(0,0)}p(\mu,\sigma_{\mathrm{pop}}\mid
x_{\mathrm{obs}},\sigma_{\mathrm{obs}})\,\mathrm{d}\mu\,\mathrm{d}\sigma_{\mathrm{pop}}\,,$
(15)

where the shorthand “$p<p(0,0)$” stands in for values of $\mu$ and
$\sigma_{\mathrm{pop}}$ such that $p(\mu,\sigma_{\mathrm{pop}}\mid
x_{\mathrm{obs}},\sigma_{\mathrm{obs}})<p(\mu=0,\sigma_{\mathrm{pop}}=0\mid
x_{\mathrm{obs}},\sigma_{\mathrm{obs}})$. Given this definition, a value of
$\mathcal{Q}_{\mathrm{GR}}=1$ means that the posterior peaks at the origin
$\mu=\sigma_{\mathrm{pop}}=0$, while $\mathcal{Q}_{\mathrm{GR}}=0$ means that
the posterior offers no support for that point (i.e., a higher value of
$\mathcal{Q}_{\mathrm{GR}}$ implies better agreement with GR). In all cases,
we simulate each catalog of $N$ observations 50 times and report medians over
the ensemble.

Unlike with BFs, there are no values of ${x}_{0}$ or $\sigma_{0}$ for which
the hierarchical analysis converges to the wrong answer given enough
observations. This is apparent from Fig. 5, which shows the value of
$\mathcal{Q}_{\mathrm{GR}}$ as a function of the deviation magnitude
(${x}_{0}/\sigma_{\mathrm{obs}}$ or $\sigma_{0}/\sigma_{\mathrm{obs}}$) and
the number of observations: even for deviations small relative to the
individual-measurement uncertainty, $\mathcal{Q}_{\mathrm{GR}}$ approaches
zero for large $N$—indicating that the posterior offers little support for
$\mu=\sigma_{\mathrm{pop}}=0$, in tension with GR. For small catalogs ($N\leq
10$), the value of $\mathcal{Q}_{\mathrm{GR}}$ is more strongly influenced by
the prior, to a greater or lesser extent depending on the magnitude of the
deviation.

Besides indicating that the data are inconsistent with the null hypothesis,
the hierarchical analysis provides descriptive information about the nature of
the deviation. With enough observations, the measurements of $\mu$ and
$\sigma_{\mathrm{pop}}$ converge to the ${x}_{0}$ and $\sigma_{0}$ values
respectively for the two models with increasing precision for larger $N$. In
Fig. 6, we show this behavior explicitly for two example magnitudes of the GR
deviation. As expected, we can detect larger deviations with fewer detections,
and need more observations to notice a nonzero scatter than a nonzero mean.

Figure 6:  Recovered population mean ($\mu$, left) and standard deviation
($\sigma_{\mathrm{pop}}$, right) as a function of catalog size ($N$,
abscissa), for our two toy models: a fixed deviation ${x}_{0}$ for all events
(left) and a deviation scatter $\sigma_{0}$ across events (right). In each
case, we show two values for the true deviation:
${x}_{0},\sigma_{0}=0.1\,\sigma_{\mathrm{obs}}$ (blue) and
${x}_{0},\sigma_{0}=0.8\,\sigma_{\mathrm{obs}}$ (red). The measurement is
represented by posterior median (solid lines) surrounded by 68% and 90%
highest-density credible bands (shading); we also show the true value (dotted
lines) and the null expectation (dashed gray line). Smaller deviations require
more observations to be detected—for example, we only need $N\gtrsim 3$ events
to notice $\mu>0$ at $1\sigma$ (68% credibility) if
${x}_{0}=0.8\,\sigma_{\mathrm{obs}}$, but $N\gtrsim 100$ if
${x}_{0}=0.1\,\sigma_{\mathrm{obs}}$.

Prior choices play a lesser role in the hierarchical approach than in BF
computations. The hierarchical analysis takes as input the likelihoods for
individual events, so the prior used to initially analyze the data is largely
irrelevant as long as it offers the likelihood ample support (e.g.,
$\Delta\gg\sigma_{\mathrm{obs}}$ in the notation of the previous section). The
choice of prior is, instead, transferred to the $\mu$ and
$\sigma_{\mathrm{pop}}$ hyperparameters; however, any reasonably smooth choice
will work assuming we have a enough events for the hierarchical measurement to
be informative. If observations are not sufficiently numerous, the result will
be influenced by the $\mu$ and $\sigma_{\mathrm{pop}}$ prior.

Figure 7:  Scaling of the uncertainty in the recovered hyperparameters $\mu$
(left) and $\sigma_{\mathrm{pop}}^{2}$ (right) as a function of catalog size
$N$, computed for a case in which the true values are $\mu=0$ and
$\sigma_{\mathrm{pop}}=\sigma_{\mathrm{obs}}/2$. A solid line marks the median
of the posterior variance computed over 50 simulated catalogs for any given
$N$, while surrounding bands enclose the corresponding 68% highest-density
interval. For low $N$, the variances
$\operatorname{\mathrm{var}}\left(\mu/\sigma_{\mathrm{obs}}\right)$ and
$\operatorname{\mathrm{var}}\left(\sigma_{\mathrm{pop}}^{2}/\sigma_{\mathrm{obs}}^{2}\right)$
are dominated by the hyperprior; as $N$ increases, the data become more
informative and the variances approach the analytic prediction of Eqs. (16)
and (17) (red dashed lines). Dotted vertical lines mark the expected number of
observations needed for the likelihood of Eq. (14) to become narrower than the
prior, as dictated by Eqs. (18) and (19) with
$\sigma_{\mathrm{prior}}=\sigma_{\mathrm{obs}}$, but rounded up to the closest
integer.

In the case of our idealized examples, we can analytically predict the number
of detections needed for the hierarchical measurement to be informative. As we
show in Appendix A, the variance of the marginalized hierarchical likelihood
is expected to scale as

$\operatorname{\mathrm{var}}\left({\mu}\right)=\frac{\sigma_{\mathrm{obs}}^{2}+\sigma_{0}^{2}}{N}$
(16)

for the inferred population mean and

$\operatorname{\mathrm{var}}\left({\sigma_{\mathrm{pop}}^{2}}\right)=\frac{2}{N}\left(\sigma_{\mathrm{obs}}^{2}+\sigma_{0}^{2}\right)^{2}$
(17)

for the inferred population variance, assuming the true population variance is
$\sigma_{0}^{2}$ and irrespective of the true mean. As a rule of thumb, the
hierarchical measurement will become informative once the characteristic width
of the likelihood of Eq. (14) becomes smaller than the scale imposed by the
prior. With hyperpriors of scale $\sigma_{\mathrm{prior}}$, this implies that
the $\mu$ measurement should start becoming informative (in the sense that we
obtain a hyperposterior narrower than the prior) once we accumulate

$N\gtrsim\frac{\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{pop}}^{2}}{\sigma_{\mathrm{prior}}^{2}}=\left(\frac{\sigma_{\mathrm{tot}}}{\sigma_{\mathrm{prior}}}\right)^{2}$
(18)

measurements; meanwhile, for $\sigma_{\mathrm{pop}}$, the equivalent
requirement is

$N\gtrsim\frac{\left(\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{pop}}^{2}\right)^{2}}{\sigma_{\mathrm{prior}}^{4}}=\left(\frac{\sigma_{\mathrm{tot}}}{\sigma_{\mathrm{prior}}}\right)^{4},$
(19)

as we show in Appendix A.

Figure 8:  As in Fig. 6, except the true population follows the fencepost
model of Sec. IV, by which the true deviation is
${x}=\pm\sigma_{\mathrm{obs}}/2$ with equal probability for either sign. Even
though this distribution cannot be expressed as the limit of a Gaussian, the
hierarchical analysis infers the correct values of $\mu=0$ and
$\sigma_{\mathrm{pop}}=\sigma_{\mathrm{obs}}/2$.

We demonstrate these scalings in Fig. 7, where we show the variance in the
inferred population mean and variance from simulated populations of
measurements with $\sigma_{0}=\sigma_{\mathrm{obs}}/2$ and no mean, and
setting $\sigma_{\mathrm{prior}}=\sigma_{\mathrm{obs}}$ for concreteness. For
increasing number of measurements $N$, the posterior converges to the true
values ($\mu=0,\sigma_{\mathrm{pop}}=\sigma_{\mathrm{obs}}/2$). For small
catalogs, i.e., for $N$ comparable or smaller than the thresholds above, the
average uncertainty in these measurements is broad but smaller than expected
simply from Eqs. (16) and (17) because it is dominated by the prior. As the
number of detections increases, the posterior variance becomes well described
by Eqs. (16) and (17), meaning that the data become informative and the
likelihood dominates over the prior.

Alternately, we may ask how many measurements would be required to establish a
non-vanishing population mean or variance. This requires
$\operatorname{\mathrm{var}}\left(\mu\right)\lesssim x_{0}^{2}$ or
$\operatorname{\mathrm{var}}\left(\sigma_{\mathrm{pop}}^{2}\right)\lesssim\sigma_{0}^{4}$.
The former would imply

$N\gtrsim\left(\frac{\sigma_{\mathrm{tot}}}{x_{0}}\right)^{2}$ (20)

and the latter

$N\gtrsim 2\left(\frac{\sigma_{\mathrm{tot}}}{\sigma_{0}}\right)^{4}.$ (21)

##  IV Fencepost model

In the examples above, the simulated populations could be perfectly reproduced
as special cases of the hierarchical population model—that is, there existed a
choice of $\mu$ and $\sigma_{\mathrm{pop}}$ for which the hierarchical model
reduced exactly to the true distribution we simulated. Of course, we do not
necessarily expect this to be the case in reality: a deviation from GR could
manifest as a nontrivial function of the source parameters and the coupling
constants in the theory; similarly, for other effects such as memory, it is
not realistic to expect the true population of parameters to be fully
described by a simple Gaussian if the null hypothesis is incorrect. In light
of that, one might worry that the hierarchical method only outperformed BFs in
the above examples because the hierarchical model was able to match the true
population exactly, and that this gain would fail to materialize in realistic
situations. Yet, as we show in this section, this is not the case: the
hierarchical method is more robust than products of BFs even when the
underlying population cannot be fit exactly by a Gaussian.

Consider a situation in which the true deviation parameter is either
${x}=\pm{x}_{0}$ for some ${x}_{0}$ and with equal probability for both signs.
The true distribution in this “fencepost” model is simply the sum of two delta
functions,

$p_{3}\left({x}\right)=\frac{1}{2}\left[\delta\left({x}-{x}_{0}\right)+\delta\left({x}+{x}_{0}\right)\right],$
(22)

and, therefore, has zero mean and standard deviation $\sigma_{0}=|{x}_{0}|$.
This population cannot be described as a limiting case of a Gaussian
distribution; nevertheless, we can always analyze it hierarchically with the
Gaussian likelihood of Eq. (14), and expect to recover the correct values for
the population mean and spread (namely, $\mu=0$ and
$\sigma_{\mathrm{pop}}=|{x}_{0}|$) given enough detections Isi  _et al._
(2019). In Fig. 8, we show this explicitly for simulated measurements in which
${x}_{0}=\sigma_{\mathrm{obs}}/2$.

On the other hand, when presented with the fencepost model, BF computations
suffer from the same problems already identified above: with a broad prior
relative to the measurement uncertainty ($\Delta>\sigma_{\mathrm{obs}}$), the
combined BF for multiple observations will necessarily converge to the wrong
answer (i.e., favor of the null hypothesis) unless
${x}_{0}\gtrsim\sigma_{\mathrm{obs}}$, in which case the deviation is
detectable in individual observations. Similarly to Fig. 2, in Fig. 9 we show
the scaling of the expected BF accumulated from fencepost-model observations,
as a function of $x_{0}$ and $\Delta$.

In Fig. 10, we compare the hierarchical and Bayes-factor results for a
progressively-higher number of simulated measurements from a fencepost
population with ${x}_{0}=\sigma_{\mathrm{obs}}/2$. As the number of
measurements increases, the hierarchical model disfavors the null hypothesis
more strongly, with $\mathcal{Q}_{\mathrm{GR}}$ vanishing rapidly for
increasing $N$. On the other hand, a BF computed with
$\Delta=5\sigma_{\mathrm{obs}}$ grows exponentially in favor of the null
hypothesis, yielding a spectacularly incorrect result.

Figure 9:  Average log BF in favor of GR accumulated for each observation in
the “fencepost” model described in Sec. IV as a function of the prior width
$\Delta$ and the location of the true deviations at $\pm{x}_{0}$. The black
lines indicate the parameter values chosen for the numerical experiment
illustrated in Fig. 10.

Figure 10:  Accumulation of credibility (left) or log BF (right) for GR versus
catalog size $N$ in numerical experiments corresponding to the parameter
choices indicated by the black lines in Fig. 9. The solid line gives the
median over 100 catalog realizations at each catalog size $N$ while the band
shows the range of the 16th to 84th percentile values. The credibility of GR
(left) is the fraction of posterior mass for $\mu$ and $\sigma$ in our
hierarchical model that lies at a lower posterior density than the GR values
$\mu=\sigma=0$ when the model is fit to a catalog of observations whose true
and observed deviations are drawn from the “fencepost” model described in Sec.
IV with ${x}=\pm\sigma_{\mathrm{obs}}/2$. The log Bayes factor (right) is the
sum of the log BFs for each observation in the catalog using a flat prior on
the true deviation $-\Delta<{x}<\Delta$ with $\Delta=5\sigma_{\mathrm{obs}}$.
The hierarchical model correctly finds that there is little credibility for GR
once the catalog size is a few hundred; the accumulated BF, on the other hand,
becomes very confident in the incorrect GR model even at small catalog sizes.

##  V Conclusions

Although conceptually appealing in idealized situations, the use of BFs to
aggregate information from multiple observations presents difficulties in
practice. Their apparent simplicity in reducing a complex model selection
problem to a single number hides an opaque dependence strict and unrealistic
population assumptions. Unless priors (aka, the “population”) adapt to the
observations at hand, BFs are difficult to interpret—a problem that is
compounded when multiplying such BFs from a catalog of observations. Even when
priors are adequate, the result on its own provides no insight as to _why_ a
model is to be preferred over another. This and related problems have been
widely discussed in the statistics literature (e.g., Lotfi  _et al._ , 2022,
and references therein), but not extensively in the context of GWs and testing
GR.

In this paper, we have examined BFs and hierarchical posteriors as two
commonly-used alternatives to derive information from collections of GW
detections in order to decide between two models, e.g., the presence or
absence of beyond-GR effects in the detected waveforms. We furthered arguments
in Zimmerman  _et al._ (2019); Isi  _et al._ (2019) to show that, without a
principled way to set priors, BFs are an unreliable tool for this task. We
demonstrated this with three examples in which the value of some parameter $x$
encoding the effect in question (e.g., a deviation from GR) follows different
distributions deviating from the null hypothesis.

We have found that, when the truth does not conform to the null model, the
usual approach of multiplying single-event BF converges to the incorrect
answer for an increasing number of observations, except in a regime where the
targeted effect is discernible in individual observations (thus negating the
need for combining events in the first place). On the other hand, hierarchical
modeling of the underlying population leads to identification of the
appropriate priors (aka, the “population model”) and converges to the correct
answer. We established this in the context of nested models for which GR can
be recovered as a special case of the beyond-GR model (i.e., $x_{0}=0$);
however, the issue of sensitivity to the prior width will still be present in
non-nested models (e.g., Laghi  _et al._ , 2021) where it will require a
different solution.

In principle, BFs could be computed after the hierarchical population
inference Lotfi  _et al._ (2022) or between different population models Abbott
_et al._ (2021a), but we here show that they are unreliable without this step.
Even then, it is not possible to evade the core problem of prior dependence
when computing BFs, no matter how many levels of inference are applied: the BF
computation based on the highest level of inference in a hierarchical model
will still depend on the choice of priors on that level, reducing the problem
once again to the choice of a prior distribution that is difficult to
establish in a principled way. This issue is devistatingly acute for the
approach that multiplies Bayes factors with a simple, fixed prior because each
observation contributes an additional prior factor.

###### Acknowledgements.

We are grateful to Tyson Littenberg for insightful feedback on this
manuscript. During part of this work, M.I. was supported by NASA through the
NASA Hubble Fellowship grant No. HST-HF2-51410.001-A awarded by the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. The Flatiron Institute is a division of the Simons Foundation,
supported through the generosity of Marilyn and Jim Simons. We thank Gregorio
Carullo for comments on non-nested models. This paper carries LIGO document
number LIGO-P2200099. Software: matplotlib Hunter (2007), julia Bezanson  _et
al._ (2012), Turing.jl Ge  _et al._ (2018), Plots.jl Breloff (2022).

##  Appendix A Expectation value and variance of hierarchical parameters

Consider $i=0\dots N-1$ measurements of parameters $x_{i}$ whose true values,
$\mu_{i}$, are drawn from a normal distribution
$\mu_{i}\sim\mathcal{N}(\mu,\,\sigma_{\mathrm{pop}})$; further assume each
measurement is unbiased, i.e., $\langle x_{i}\rangle=\mu_{i}$, where the angle
brackets denote a noise average, and that it is well represented by a Gaussian
distribution such that $x_{i}\sim\mathcal{N}(\mu_{i},\,\sigma_{i})$, where
$\sigma_{i}$ is the measurement uncertainty.555In the main text, we used
slightly different notation: instead of
$\left(x_{i},\mu_{i},\sigma_{i}\right)$ we had
$\left(x_{\mathrm{obs}},{x},\sigma_{\mathrm{obs}}\right)$; the former is
slightly more succinct, which will be helpful here given the increased number
of mathematical expressions.

The joint likelihood for $\mu$ and $\sigma_{\mathrm{pop}}$ conditional on the
$N$ uncertainties $\sigma_{i}$ can be obtained by marginalizing over the true
values $\mu_{i}$:

$\displaystyle p(x_{i}\mid\mu,\sigma_{\mathrm{pop}},\sigma_{i})$
$\displaystyle=\int
p(x_{i}\mid\mu_{i},\sigma_{i})\,p(\mu_{i}\mid\mu,\sigma_{\mathrm{pop}})\,\mathrm{d}\mu_{i}$
$\displaystyle=\frac{1}{2\pi\sigma_{i}\sigma_{\mathrm{pop}}}\int
e^{-\frac{\left(x_{i}-\mu_{i}\right)^{2}}{2\sigma_{i}^{2}}}e^{-\frac{\left(\mu_{i}-\mu\right)^{2}}{2\sigma_{\mathrm{pop}}^{2}}}\,\mathrm{d}\mu_{i}$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma^{2}_{\mathrm{tot},i}}}e^{-\frac{\left(x_{i}-\mu\right)^{2}}{2\sigma_{\mathrm{tot},i}^{2}}}\,,$
(23)

where $\sigma_{\mathrm{tot},i}^{2}=\sigma_{i}^{2}+\sigma_{\mathrm{pop}}^{2}$
is the total variance for the $i$th measurement. For the full set of $N$
measurements $\\{x_{i}\\}$, then, the hierarchical likelihood is

$p(\\{x_{i}\\}\mid\mu,\sigma_{\mathrm{pop}},\\{\sigma_{i}\\})=\prod_{i=0}^{N-1}\frac{1}{\sqrt{2\pi\sigma_{\mathrm{tot},i}^{2}}}e^{-\frac{\left(x_{i}-\mu\right)^{2}}{2\sigma_{\mathrm{tot},i}^{2}}}.$
(24)

The maximum-likelihood estimators for $\mu$ and $\sigma_{\mathrm{pop}}$,
denoted $\hat{\mu}$ and $\hat{\sigma}_{\mathrm{pop}}$ respectively, can be
found by enforcing

$\left.\partial_{\mu}\ln\left[p(\\{x_{i}\\}\mid\mu,\sigma_{\mathrm{pop}},\\{\sigma_{i}\\})\right]\right|_{\mu=\hat{\mu}}=0\,,$
(25)

where $\partial_{\mu}$ denotes the partial derivative
$\frac{\partial}{\partial\mu}$, and similar for $\sigma_{\mathrm{pop}}$. We
find that $\hat{\mu}$ and $\hat{\sigma}_{\mathrm{pop}}$ satisfy

$\hat{\mu}=\frac{\sum x_{i}w_{i}}{\sum w_{i}}\,,$ (26)

and

$\sum w_{i}-\sum w_{i}^{2}\left(x_{i}-\hat{\mu}\right)^{2}=0$ (27)

for
$w_{i}\equiv\sigma_{\mathrm{tot},i}^{-2}=\left(\sigma_{i}^{2}+\hat{\sigma}_{\mathrm{pop}}^{2}\right)^{-1}$.

We cannot solve this most general (heteroskedastic) case for $\hat{\mu}$ and
$\hat{\sigma}_{\mathrm{pop}}$ in closed form, so we specialize to the
(homoskedastic) case where all measurements have similar error,
$\sigma_{i}=\sigma_{\mathrm{obs}}$ for all $i$. With this simplification, the
above relations reduce to

$\hat{\mu}=\frac{1}{N}\sum x_{i}\,,$ (28)

and

$\hat{\sigma}_{\mathrm{pop}}^{2}=\frac{1}{N}\sum\left(x_{i}-\hat{\mu}\right)^{2}-\sigma_{\mathrm{obs}}^{2}\,.$
(29)

As might be expected, the maximum-likelihood estimate of the population mean
is simply the sample mean, and the inferred population variance corresponds to
the variance in the data that cannot be accounted for by the statistical
uncertainty in each individual measurement.

We can go one step further and compute the uncertainty in these estimators,
which we can take as a proxy for the width of the marginal likelihoods. The
uncertainty in $\hat{\mu}$ is straightforward to compute, since this is just a
linear combination of independent random variables $x_{i}$ with known variance
$\sigma_{\mathrm{tot}}^{2}\equiv\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{pop}}^{2}$,
hence

$\operatorname{\mathrm{var}}\left(\hat{\mu}\right)=\frac{\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{pop}}^{2}}{N}\,.$
(30)

Obtaining
$\operatorname{\mathrm{var}}\left(\hat{\sigma}_{\mathrm{pop}}^{2}\right)$ is
less straightforward, but we can do so by writing
$\operatorname{\mathrm{var}}\left(\hat{\sigma}_{\mathrm{pop}}^{2}\right)=F^{-1}$,
in terms of the corresponding Fisher element666We carry out this calculation
for $\sigma_{\mathrm{pop}}^{2}$ instead of $\sigma_{\mathrm{pop}}$ because the
latter is irregular at the origin.

$F\equiv-\left\langle\left.\partial^{2}_{\sigma_{\mathrm{pop}}^{2}}\ln\left[p(\\{x_{i}\\}\mid\mu,\sigma_{\mathrm{pop}},\sigma_{\mathrm{obs}})\right]\right|_{\sigma_{\mathrm{pop}}^{2}=\hat{\sigma}_{\mathrm{pop}}^{2}}\right\rangle.$
(31)

Doing the math, we find

$\operatorname{\mathrm{var}}\left(\hat{\sigma}_{\mathrm{pop}}^{2}\right)=2\frac{\left(\sigma_{\mathrm{obs}}^{2}+\hat{\sigma}_{\mathrm{pop}}^{2}\right)^{2}}{N}\
.$ (32)

In the main text, we quoted these results for
$\operatorname{\mathrm{var}}\left(\hat{\mu}\right)$ and
$\operatorname{\mathrm{var}}\left(\hat{\sigma}_{\mathrm{pop}}^{2}\right)$ in
Eqs. (16) and (17) respectively.

We can compare the width of the hierarchical likelihood, as proxied by the
estimator variances above, to some typical scale of interest in the problem.
Below we consider the scale imposed by the $\mu$ and $\sigma_{\mathrm{pop}}$
hyperpriors to estimate the number of events before the likelihood become
informative with respect to the prior. We chose the hyperpriors to be
Gaussians with scale $\sigma_{\mathrm{prior}}$, restricting to positive
$\sigma_{\mathrm{pop}}$ values, i.e.,

$p(\mu\mid\sigma_{\mathrm{prior}})=\frac{1}{\sqrt{2\pi\sigma_{\mathrm{prior}}^{2}}}\exp\left(-\frac{\mu^{2}}{2\sigma_{\mathrm{prior}}^{2}}\right)$
(33)

for the mean, and

$p(\sigma_{\mathrm{pop}}\mid\sigma_{\mathrm{prior}})=\begin{cases}\sqrt{\frac{2}{\pi\sigma_{\mathrm{prior}}^{2}}}\exp\left(-\frac{\sigma_{\mathrm{pop}}^{2}}{2\sigma_{\mathrm{prior}}^{2}}\right)\hfill(\sigma_{\mathrm{pop}}\geq
0)\\\ 0\hfill(\sigma_{\mathrm{pop}}<0)\end{cases}$ (34)

for the standard deviation, where the difference in normalization arises from
the $\sigma_{\mathrm{pop}}\geq 0$ truncation. The prior variances in our
example are, thus,
$\operatorname{\mathrm{var}}\left(\mu\right)=\sigma_{\mathrm{prior}}^{2}$ for
the mean, and
$\operatorname{\mathrm{var}}\left(\sigma_{\mathrm{pop}}^{2}\right)=2\sigma_{\mathrm{prior}}^{4}$
for the variance (obtained through direct computation). For concreteness, in
the main text we set $\sigma_{\mathrm{prior}}=\sigma_{\mathrm{obs}}$, since
that is the only scale intrinsic to the measurement.

We can now directly compute the number of observations required for the
likelihood to achieve comparable widths by equating these variances to the
likelihood variances from Eqs. (30) and (32) above. The result is

$N\gtrsim\frac{\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{pop}}^{2}}{\sigma_{\mathrm{prior}}^{2}}$
(35)

for $\mu$, and

$N\gtrsim\left(\frac{\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathrm{pop}}^{2}}{\sigma_{\mathrm{prior}}^{2}}\right)^{2}=\left(\frac{\sigma_{\mathrm{tot}}}{\sigma_{\mathrm{prior}}}\right)^{4}$
(36)

for $\sigma_{\mathrm{pop}}$. We quoted these results in Eqs. (18) and (19) in
the main text. We require at least this many measurements before the
uncertainty in the population variance can be smaller than the measurement
uncertainty. Until that point, the $\sigma_{\mathrm{pop}}$ posterior will be
dominated by the prior.

The results for the $N$ thresholds quoted above hinge on the specific choice
of prior for $\mu$ and $\sigma_{\mathrm{pop}}$. In the main text, we justified
our decision to set the scale of those priors based on $\sigma_{\mathrm{obs}}$
by noting that this is the only intrinsic scale to the problem, and should
always be sufficiently broad as long as the deviation from GR is not visible
in a single detection—the regime in which we are interested in the first
place. Had we chosen to increase the prior variance by some factor, then the
$N$ thresholds would decrease by the same factor, i.e., we need fewer
detections to gain information relative to a broad (less informative) prior
than a narrower prior. Either way, the result converges to the right answer as
we accumulate more observations.

## References

  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific Collaboration), Classical Quantum Gravity 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese _et al._ (Virgo Collaboration), Classical Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 913, L7 (2021a), arXiv:2010.14533 [astro-ph.HE] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2111.03634 [astro-ph.HE] .
  * Isi  _et al._ (2019) M. Isi, K. Chatziioannou, and W. M. Farr, Phys. Rev. Lett. 123, 121101 (2019), arXiv:1904.08011 [gr-qc] .
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 103, 122002 (2021c), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ (2021d) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021d), arXiv:2112.06861 [gr-qc] .
  * Li  _et al._ (2012a) T. G. F. Li, W. Del Pozzo, S. Vitale, C. Van Den Broeck, M. Agathos, J. Veitch, K. Grover, T. Sidery, R. Sturani, and A. Vecchio, Phys. Rev. D 85, 082003 (2012a), arXiv:1110.0530 [gr-qc] .
  * Li  _et al._ (2012b) T. G. F. Li, W. Del Pozzo, S. Vitale, C. Van Den Broeck, M. Agathos, J. Veitch, K. Grover, T. Sidery, R. Sturani, and A. Vecchio, J. Phys. Conf. Ser. 363, 012028 (2012b), arXiv:1111.5274 [gr-qc] .
  * Agathos  _et al._ (2014) M. Agathos, W. Del Pozzo, T. G. F. Li, C. Van Den Broeck, J. Veitch, and S. Vitale, Phys. Rev. D 89, 082001 (2014), arXiv:1311.0420 [gr-qc] .
  * Sampson  _et al._ (2013) L. Sampson, N. Cornish, and N. Yunes, Phys. Rev. D 87, 102001 (2013), arXiv:1303.1185 [gr-qc] .
  * Sampson  _et al._ (2014) L. Sampson, N. Cornish, and N. Yunes, Phys. Rev. D 89, 064037 (2014), arXiv:1311.4898 [gr-qc] .
  * Meidam  _et al._ (2014) J. Meidam, M. Agathos, C. Van Den Broeck, J. Veitch, and B. S. Sathyaprakash, Phys. Rev. D 90, 064009 (2014).
  * Chatziioannou  _et al._ (2019) K. Chatziioannou _et al._ , Phys. Rev. D 100, 104015 (2019), arXiv:1903.06742 [gr-qc] .
  * Hübner  _et al._ (2020) M. Hübner, C. Talbot, P. D. Lasky, and E. Thrane, Phys. Rev. D 101, 023011 (2020), arXiv:1911.12496 [astro-ph.HE] .
  * Del Pozzo  _et al._ (2013) W. Del Pozzo, T. G. F. Li, M. Agathos, C. Van Den Broeck, and S. Vitale, Phys. Rev. Lett. 111, 071101 (2013), arXiv:1307.8338 [gr-qc] .
  * Chatziioannou  _et al._ (2015) K. Chatziioannou, K. Yagi, A. Klein, N. Cornish, and N. Yunes, Phys. Rev. D 92, 104008 (2015), arXiv:1508.02062 [gr-qc] .
  * Abbott  _et al._ (2020) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Class. Quant. Grav. 37, 045006 (2020), arXiv:1908.01012 [gr-qc] .
  * Ghosh  _et al._ (2021) S. Ghosh, X. Liu, J. Creighton, I. Magaña Hernandez, W. Kastaun, and G. Pratten, (2021), arXiv:2104.08681 [gr-qc] .
  * Isi  _et al._ (2017) M. Isi, M. Pitkin, and A. J. Weinstein, Phys. Rev. D 96, 042001 (2017), arXiv:1703.07530 [gr-qc] .
  * Callister  _et al._ (2017) T. Callister, A. S. Biscoveanu, N. Christensen, M. Isi, A. Matas, O. Minazzoli, T. Regimbau, M. Sakellariadou, J. Tasson, and E. Thrane, Phys. Rev. X 7, 041058 (2017), arXiv:1704.08373 [gr-qc] .
  * Liu  _et al._ (2021) X. Liu, I. M. Hernandez, and J. Creighton, Astrophys. J. 908, 97 (2021), arXiv:2009.06539 [astro-ph.HE] .
  * Lo and Magaña Hernandez (2021) R. K. L. Lo and I. Magaña Hernandez, (2021), arXiv:2104.09339 [gr-qc] .
  * Ashton  _et al._ (2020) G. Ashton, K. Ackley, I. M. n. Hernandez, and B. Piotrzkowski, (2020), arXiv:2009.12346 [astro-ph.HE] .
  * Ota and Chirenti (2021) I. Ota and C. Chirenti, (2021), arXiv:2108.01774 [gr-qc] .
  * Veitch and Vecchio (2008) J. Veitch and A. Vecchio, Class. Quantum Grav. 25, 184010 (2008), arXiv:0807.4483 [gr-qc] .
  * Veitch and Vecchio (2010) J. Veitch and A. Vecchio, Phys.Rev. D81, 062003 (2010), arXiv:0911.3820 [astro-ph.CO] .
  * Cornish and Littenberg (2015) N. J. Cornish and T. B. Littenberg, Classical Quantum Gravity 32, 135012 (2015), arXiv:1410.3835 [gr-qc] .
  * Isi  _et al._ (2018) M. Isi, R. Smith, S. Vitale, T. J. Massinger, J. Kanner, and A. Vajpeyi, Phys. Rev. D 98, 042007 (2018), arXiv:1803.09783 [gr-qc] .
  * Ashton  _et al._ (2019) G. Ashton, E. Thrane, and R. J. E. Smith, Phys. Rev. D 100, 123018 (2019), arXiv:1909.11872 [gr-qc] .
  * Pratten and Vecchio (2021) G. Pratten and A. Vecchio, Phys. Rev. D 104, 124039 (2021), arXiv:2008.00509 [gr-qc] .
  * Cornish  _et al._ (2021) N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse, Phys. Rev. D 103, 044006 (2021), arXiv:2011.09494 [gr-qc] .
  * Bustillo  _et al._ (2021) J. C. Bustillo, P. D. Lasky, and E. Thrane, Phys. Rev. D 103, 024041 (2021), arXiv:2010.01857 [gr-qc] .
  * Zimmerman  _et al._ (2019) A. Zimmerman, C.-J. Haster, and K. Chatziioannou, Phys. Rev. D 99, 124044 (2019), arXiv:1903.11008 [astro-ph.IM] .
  * Kass and Raftery (1995) R. E. Kass and A. E. Raftery, J. Am. Statist. Assoc. 90, 773 (1995).
  * James and Stein (1961) W. James and C. Stein, in  _Proc. 4th Berkeley Sympos. Math. Statist. and Prob., Vol. I_ (Univ. California Press, Berkeley, Calif., 1961) pp. 361–379.
  * Lindley and Smith (1972) D. V. Lindley and A. F. M. Smith, Journal of the Royal Statistical Society. Series B (Methodological) 34, 1 (1972).
  * Efron and Morris (1977) B. Efron and C. Morris, Scientific American 236, 119 (1977).
  * Rubin (1981) D. B. Rubin, Journal of Educational Statistics 6, 377 (1981).
  * Blanchet and Sathyaprakash (1995) L. Blanchet and B. Sathyaprakash, Phys. Rev. Lett. 74, 1067 (1995).
  * Blanchet and Sathyaprakash (1994) L. Blanchet and B. Sathyaprakash, Class. Quant. Grav. 11, 2807 (1994).
  * Arun  _et al._ (2006a) K. G. Arun, B. R. Iyer, M. S. S. Qusailah, and B. S. Sathyaprakash, Phys. Rev. D 74, 024006 (2006a), arXiv:gr-qc/0604067 .
  * Arun  _et al._ (2006b) K. G. Arun, B. R. Iyer, M. S. S. Qusailah, and B. S. Sathyaprakash, Classical Quantum Gravity 23, L37 (2006b), arXiv:gr-qc/0604018 .
  * Yunes and Pretorius (2009) N. Yunes and F. Pretorius, Phys. Rev. D 80, 122003 (2009), arXiv:0909.3328 [gr-qc] .
  * Mishra  _et al._ (2010) C. K. Mishra, K. G. Arun, B. R. Iyer, and B. S. Sathyaprakash, Phys. Rev. D 82, 064010 (2010), arXiv:1005.0304 [gr-qc] .
  * Li  _et al._ (2012a) T. Li, W. Del Pozzo, S. Vitale, C. Van Den Broeck, M. Agathos, _et al._ , Phys.Rev. D85, 082003 (2012a), arXiv:1110.0530 [gr-qc] .
  * Li  _et al._ (2012b) T. G. F. Li, W. Del Pozzo, S. Vitale, C. Van Den Broeck, M. Agathos, J. Veitch, K. Grover, T. Sidery, R. Sturani, and A. Vecchio,  _Gravitational waves. Numerical relativity - data analysis. Proceedings, 9th Edoardo Amaldi Conference, Amaldi 9, and meeting, NRDA 2011, Cardiff, UK, July 10-15, 2011_ , J. Phys. Conf. Ser. 363, 012028 (2012b), arXiv:1111.5274 [gr-qc] .
  * Agathos  _et al._ (2014) M. Agathos, W. Del Pozzo, T. G. F. Li, C. Van Den Broeck, J. Veitch, and S. Vitale, Phys. Rev. D 89, 082001 (2014), arXiv:1311.0420 [gr-qc] .
  * Dickey (1971) J. M. Dickey, The Annals of Mathematical Statistics 42, 204 (1971).
  * Verdinelli and Wasserman (1995) I. Verdinelli and L. Wasserman, Journal of the American Statistical Association 90, 614 (1995).
  * Perkins and Yunes (2022) S. Perkins and N. Yunes, (2022), arXiv:2201.02542 [gr-qc] .
  * Lotfi  _et al._ (2022) S. Lotfi, P. Izmailov, G. Benton, M. Goldblum, and A. G. Wilson, arXiv e-prints , arXiv:2202.11678 (2022), arXiv:2202.11678 [cs.LG] .
  * Laghi  _et al._ (2021) D. Laghi, G. Carullo, J. Veitch, and W. Del Pozzo, Classical and Quantum Gravity 38, 095005 (2021), arXiv:2011.03816 [gr-qc] .
  * Hunter (2007) J. D. Hunter, Computing In Science & Engineering 9, 90 (2007).
  * Bezanson  _et al._ (2012) J. Bezanson, S. Karpinski, V. B. Shah, and A. Edelman, arXiv e-prints (2012), 10.48550/arXiv.1209.5145, arXiv:1209.5145 [cs.PL] .
  * Ge  _et al._ (2018) H. Ge, K. Xu, and Z. Ghahramani, in  _International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain_ (2018) pp. 1682–1690.
  * Breloff (2022) T. Breloff, “Plots.jl,”  (2022).