We know that zones IV-1 and II-1 have the same convergence performance because
of the odd function property in the sliding mode. Therefore, combing (115),
(117) and (123), for the zones IV-1 (where, $-e_{1}\left(t_{c}\right)\geq
e_{2}\left(t_{c}\right)>0$) and II-1 (where, $e_{1}\left(t_{c}\right)\geq-
e_{2}\left(t_{c}\right)>0$), we get the non-overshooting convergence
conditions, when $e_{1}\left(t_{c}\right)e_{2}\left(t_{c}\right)<0$ and
$\left|e_{1}\left(t_{c}\right)\right|\geq\left|e_{2}\left(t_{c}\right)\right|$,
as follows:

$\displaystyle k_{1}$ $\displaystyle\in$
$\displaystyle\left(\frac{\left|e_{2}\left(t_{c}\right)\right|}{\left|e_{1}\left(t_{c}\right)\right|},\infty\right)$
(160) $\displaystyle\bar{k}_{2}$ $\displaystyle>$
$\displaystyle\max\left\\{k_{1}\left|e_{2}\left(t_{c}\right)\right|,\frac{k_{1}^{2}}{3}\left[\left|e_{1}\left(t_{c}\right)\right|+\sqrt{e_{1}^{2}\left(t_{c}\right)+3\left(\frac{e_{2}\left(t_{c}\right)}{k_{1}}\right)^{2}}\right]\right\\}$
(161)

Furthermore, considering the disturbance $d\left(t\right)$, the conditions of
parameters selection, when $e_{1}\left(t_{c}\right)e_{2}\left(t_{c}\right)<0$
and
$\left|e_{1}\left(t_{c}\right)\right|\geq\left|e_{2}\left(t_{c}\right)\right|$,
are expressed as follows:

$\displaystyle k_{1}$ $\displaystyle\in$
$\displaystyle\left(\frac{\left|e_{2}\left(t_{c}\right)\right|}{\left|e_{1}\left(t_{c}\right)\right|},\infty\right)$
(162) $\displaystyle k_{2}$ $\displaystyle>$
$\displaystyle\max\left\\{k_{1}\left|e_{2}\left(t_{c}\right)\right|+L_{d},\frac{k_{1}^{2}}{3}\left[\left|e_{1}\left(t_{c}\right)\right|+\sqrt{e_{1}^{2}\left(t_{c}\right)+3\left(\frac{e_{2}\left(t_{c}\right)}{k_{1}}\right)^{2}}\right]+L_{d}\right\\},$
(163)

Therefore, the sliding mode is non-overshooting stable for
$t\in\left[t_{c},t_{c}+t_{s}\right)$, and $\sigma\left(t\right)=0$ holds for
$t\geq t_{c}+t_{s}$. Then, the linear convergence law
$\dot{e}_{1}\left(t\right)=-k_{1}e_{1}\left(t\right)$ makes
$\underset{t\rightarrow\infty}{\lim}e_{1}\left(t\right)=0$ without overshoot.
In addition, from $e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)=0$ for $t\geq
t_{c}+t_{s}$, we get
$\underset{t\rightarrow\infty}{\lim}e_{2}\left(t\right)=0$.

In general, for $e_{1}\left(t_{c}\right)$ and $e_{2}\left(t_{c}\right)$ in
zones IV-1 and II-1, when $k_{1}$ and $k_{2}$ are selected from (126) and
(127), the system is exponentially stable, and no overshoot exists for the
variable $e_{1}\left(t\right)$. This confirms the convergence curves of
sliding variables $e_{1}\left(t\right)$, $e_{2}\left(t\right)$ and function
$\sigma\left(t\right)$ described in Fig. 14 (b)-1 and (b)-2 for the zones IV-1
and II-1, respectively.

_3) For zone III, i.e., III-1:_ $-e_{2}\left(t_{c}\right)\geq-
e_{1}\left(t_{c}\right)\geq 0$_and III-2:_
$-e_{1}\left(t_{c}\right)>-e_{2}\left(t_{c}\right)\geq 0$

From Figure 13, we know that the corresponding symmetrical zone of III is zone
I; and the corresponding symmetrical zone of III-2 is zone I-2.

Figure 17: Arranged trajectories of $e_{1}\left(t\right)$ and
$\sigma\left(t\right)$ for $e_{1}\left(t\right)$ non-overshooting convergence
for $t\in\left[t_{c},t_{c}+t_{s}\right)$ in range III-1:
$-e_{2}\left(t_{c}\right)\geq-e_{1}\left(t_{c}\right)\geq 0$ and range III-2:
$-e_{1}\left(t_{c}\right)>-e_{2}\left(t_{c}\right)\geq 0$.

For zone III, we can get the conditions of non-overshooting convergence for
$t\in\left[t_{c},t_{c}+t_{s}\right)$:

1) $\bar{k}_{2}>k_{1}\left|e_{2}\left(t\right)\right|$ [finite-time stability
to get convergence law $\dot{e}_{1}\left(t\right)=-k_{1}e_{1}\left(t\right)$];

2) the sliding function
$\sigma\left(t\right)=e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)<0$ for
$t\in\left[t_{c},t_{c}+t_{s}\right)$, and
$\sigma\left(t_{c}+t_{s}\right)=e_{2}\left(t_{c}+t_{s}\right)+k_{1}e_{1}\left(t_{c}+t_{s}\right)=0$;

3) $e_{1}\left(t\right)<0$ [ $e_{1}\left(t\right)$ does not go beyond zero].

In the following, we will determine $k_{1}$ and $k_{2}$ to satisfy these
conditions.

From (90), for $e_{1}\left(t\right)$ (a segment of a parabola) for
$t\in\left[t_{c},t_{c}+t_{s}\right)$, we get its axis of symmetry

$-\frac{b_{1}}{2a_{1}}=-\frac{e_{2}\left(t_{c}\right)}{\bar{k}_{2}}\geq 0$
(164)

because $e_{2}\left(t_{c}\right)\leq 0$. Also, its vertex satisfies

$\frac{4a_{1}c_{1}-b_{1}^{2}}{4a_{1}}=\frac{2\bar{k}_{2}e_{1}\left(t_{c}\right)-e_{2}^{2}\left(t_{c}\right)}{2\bar{k}_{2}}<0$
(165)

because $e_{1}\left(t_{c}\right)\leq 0$ (See $e_{1}\left(t\right)$ in Figure
17(a)).

From (91), for
$\sigma\left(t\right)=e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)$ (the
segment of a parabola) for $t\in\left[t_{c},t_{c}+t_{s}\right)$, we get its
axis of symmetry:

$-\frac{b}{2a}=-\frac{\bar{k}_{2}+k_{1}e_{2}\left(t_{c}\right)}{k_{1}\bar{k}_{2}}<0$
(166)

because of the finite-time convergence condition

$\bar{k}_{2}>k_{1}\left|e_{2}\left(t_{c}\right)\right|$ (167)

In order to make
$\sigma\left(t\right)=e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)<0$ for
$t\in\left[t_{c},t_{c}+t_{s}\right)$ (See Figure 17(b)), its vertex needs to
satisfy

$\frac{4ac-b^{2}}{4a}=\frac{2k_{1}\bar{k}_{2}\left(e_{2}\left(t_{c}\right)+k_{1}e_{1}\left(t_{c}\right)\right)-\left(\bar{k}_{2}+k_{1}e_{2}\left(t_{c}\right)\right)^{2}}{2k_{1}\bar{k}_{2}}<0$
(168)

To satisfy (132), it should be that
$\sigma\left(t_{c}\right)=e_{2}\left(t_{c}\right)+k_{1}e_{1}\left(t_{c}\right)<0$.
Therefore, for zones III, we get the $k_{1}$ condition as

$k_{1}>0$ (169)

Comparing the settling time $t_{c}+t_{s}$ in (97) and
$\left.t\right|_{e_{1}\left(t\right)=0}$ in (94), we get

$t_{c}+t_{s}<\left.t\right|_{e_{1}\left(t\right)=0}$ (170)

Therefore, we get

$e_{1}\left(t\right)<0\text{ for }t\in\left[t_{c},t_{c}+t_{s}\right)$ (171)

For $e_{2}\left(t\right)$, when $t=t_{c}+t_{s}$, from (89) and (97), we get

$\left.e_{2}\left(t\right)\right|_{t=t_{c}+t_{s}}=e_{2}\left(t_{c}\right)-\bar{k}_{2}t_{s}=-\frac{\bar{k}_{2}}{k_{1}}+\sqrt{e_{2}^{2}\left(t_{c}\right)+\left(\frac{\bar{k}_{2}}{k_{1}}\right)^{2}-2\bar{k}_{2}e_{1}\left(t_{c}\right)}$
(172)

For $t\in\left[t_{c},t_{c}+t_{s}\right)$, considering the finite-time
convergence condition, we need

$\bar{k}_{2}>k_{1}\max\left\\{\left|e_{2}\left(t\right)\right|\right\\}=k_{1}\left|\left.e_{2}\left(t\right)\right|_{t=t_{c}+t_{s}}\right|$
(173)

Combining (136) and (137), we get

$\bar{k}_{2}>k_{1}\left|\left.e_{2}\left(t\right)\right|_{t=t_{c}+t_{s}}\right|=k_{1}\left|-\frac{\bar{k}_{2}}{k_{1}}+\sqrt{e_{2}^{2}\left(t_{c}\right)+\left(\frac{\bar{k}_{2}}{k_{1}}\right)^{2}-2\bar{k}_{2}e_{1}\left(t_{c}\right)}\right|$
(174)

Therefore, $\bar{k}_{2}$ should satisfy

$\bar{k}_{2}>\frac{k_{1}^{2}}{3}\left[\left|e_{1}\left(t_{c}\right)\right|+\sqrt{e_{1}^{2}\left(t_{c}\right)+3\left(\frac{e_{2}\left(t_{c}\right)}{k_{1}}\right)^{2}}\right]$
(175)

We know that zones III and I have the same convergence performance because of
the odd function property in the sliding mode. Therefore, combing (131), (133)
and (139), for the zones III and I, we get the non-overshooting convergence
conditions:

$\displaystyle k_{1}$ $\displaystyle\in$ $\displaystyle\left(0,\infty\right)$
(176) $\displaystyle\bar{k}_{2}$ $\displaystyle>$
$\displaystyle\max\left\\{k_{1}\left|e_{2}\left(t_{c}\right)\right|,\frac{k_{1}^{2}}{3}\left[\left|e_{1}\left(t_{c}\right)\right|+\sqrt{e_{1}^{2}\left(t_{c}\right)+3\left(\frac{e_{2}\left(t_{c}\right)}{k_{1}}\right)^{2}}\right]\right\\}$
(177)

Furthermore, considering the disturbance $d\left(t\right)$, the conditions of
parameters selection are expressed as follows:

$\displaystyle k_{1}$ $\displaystyle\in$ $\displaystyle\left(0,\infty\right)$
(178) $\displaystyle k_{2}$ $\displaystyle>$
$\displaystyle\max\left\\{k_{1}\left|e_{2}\left(t_{c}\right)\right|+L_{d},\frac{k_{1}^{2}}{3}\left[\left|e_{1}\left(t_{c}\right)\right|+\sqrt{e_{1}^{2}\left(t_{c}\right)+3\left(\frac{e_{2}\left(t_{c}\right)}{k_{1}}\right)^{2}}\right]+L_{d}\right\\}$
(179)

Therefore, the sliding mode is non-overshooting stable for
$t\in\left[t_{c},t_{c}+t_{s}\right)$, and $\sigma\left(t\right)=0$ holds for
$t\geq t_{c}+t_{s}$. Then, the linear convergence law
$\dot{e}_{1}\left(t\right)=-k_{1}e_{1}\left(t\right)$ makes
$\underset{t\rightarrow\infty}{\lim}e_{1}\left(t\right)=0$ without overshoot.
In addition, from $e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)=0$ for $t\geq
t_{c}+t_{s}$, we get
$\underset{t\rightarrow\infty}{\lim}e_{2}\left(t\right)=0$.

In general, for $e_{1}\left(t_{c}\right)$ and $e_{2}\left(t_{c}\right)$ in
zones III and I, when $k_{1}$ and $k_{2}$ are selected from (142) and (143),
the system is exponentially stable, and no overshoot exists for the variable
$e_{1}\left(t\right)$. This confirms the convergence curves of sliding
variables $e_{1}\left(t\right)$, $e_{2}\left(t\right)$ and function
$\sigma\left(t\right)$ described in Fig. 13 (c)-1 and (c)-2 for the zones
III-1 and I-1, and Fig. 13 (d)-1 and (d)-2 for zones III-2 and I-2,
respectively.

Finally, combing parameter selection conditions (110)-(111), (126)-(127), and
(142)-(143) in the different zones, we can get the parameter conditions (19)
and (20) for non-overshooting stable system in Theorem 3.1.

C. Determination of $e_{1c}$ and $e_{2c}$ for bounded system gain and without
overshoot

From the expression of $k_{2}$ in (20), in order to make the system gain
bounded, we need to reduce $\left|e_{1}\left(t_{c}\right)\right|$,
$\left|e_{2}\left(t_{c}\right)\right|$ and
$\frac{e_{2}^{2}\left(t_{c}\right)}{2\left|e_{1}\left(t_{c}\right)\right|}$.
For the first subsystem, according to the robust and non-overshooting
reachability, we get $e_{1c}=e_{1}(t_{c})$ and $e_{2c}=e_{2}(t_{c})$ when
$t=t_{c}$. Therefore, we can select the bounded $e_{1c}$ and $e_{2c}$ for
bounded $\left|e_{1}\left(t_{c}\right)\right|$ and
$\left|e_{2}\left(t_{c}\right)\right|$. Furthermore, we hope to get the fast
convergence as the mode in the zones II-1 and IV-1, and the up-bound of
$k_{2}$ from the system gain limitation is considered. In order to get the
non-overshooting stability with the bounded system gain, we select

$\displaystyle e_{1c}$ $\displaystyle\in$
$\displaystyle\left(0,k_{2M}-L_{d}\right)$ $\displaystyle e_{2c}$
$\displaystyle\in$
$\displaystyle\left(e_{1c},\sqrt{\left(k_{2M}-L_{d}\right)e_{1c}}\right]$
(180)

where, $k_{2M}$ is the up-bound of $k_{2}$. In fact, due to
$e_{1c}=e_{1}(t_{c})$ and $e_{2c}=e_{2}(t_{c})$ when $t=t_{c}$, we get the
parameter $k_{1}\in\left(0,\frac{e_{2c}}{e_{1c}}\right)$ for the second
subsystem. Therefore, there exists $\beta_{11}\in(0,1)$ such that
$k_{1}=\beta_{11}\frac{e_{2c}}{e_{1c}}$. In (20), we know that

$k_{1}e_{2c}\leq k_{2M}-L_{d}$ (181)

From $k_{1}=\beta_{11}\frac{e_{2c}}{e_{1c}}$ and (145), we can get

$\beta_{11}\frac{e_{2c}^{2}}{e_{1c}}\leq k_{2M}-L_{d}$ (182)

In (20), we also know that

$\frac{e_{2c}^{2}}{2e_{1c}}\leq k_{2M}-L_{d}$ (183)

Combing (146) and (147), we can select

$\max\left\\{\beta_{11},\frac{1}{2}\right\\}\frac{e_{2c}^{2}}{e_{1c}}\leq
k_{2M}-L_{d}$ (184)

Because the system gain limitation $0<k_{2}\leq k_{2M}$, i.e., $k_{2M}$ is the
maximum implementation of $k_{2}$. Because $\beta_{11}\in(0,1)$, we get
$\max\left\\{\beta_{11},\frac{1}{2}\right\\}<1$. Therefore, from
$e_{1c}<e_{2c}$, for (148), we select

$e_{1c}<e_{2c}\leq\sqrt{\left(k_{2M}-L_{d}\right)e_{1c}}$ (185)

i.e., $e_{1c}^{2}<\left(k_{2M}-L_{d}\right)e_{1c}$ and
$e_{1c}<e_{2c}\leq\sqrt{\left(k_{2M}-L_{d}\right)e_{1c}}$. Then, we get (144).

This concludes the proof. $\blacksquare$

_Proof of Theorem 3.2:_

(i) Firstly, we consider $\left|e_{1}\left(t\right)\right|>e_{1c}$.

Case one: If $e_{1}\left(t\right)>e_{1c}$, for (24), we get

$\frac{d\left(e_{2}\left(t\right)+e_{2c}\right)}{dt}=-k_{c}\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\right)\right]+d(t)$
(186)

A Lyapunov function candidate is selected as

$V_{c}=\frac{1}{2}\left(e_{2}\left(t\right)+e_{2c}\right)^{2}$ (187)

Then, taking the derivative for $V_{c}$, we get

$\displaystyle\dot{V}_{c}$ $\displaystyle=$
$\displaystyle\left(e_{2}\left(t\right)+e_{2c}\right)\left\\{-k_{c}\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\right)\right]+d(t)\right\\}$
(188) $\displaystyle\leq$
$\displaystyle-\left(k_{c}\left|\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\right)\right]\right|-L_{d}\right)\left|e_{2}\left(t\right)+e_{2c}\right|$
$\displaystyle=$
$\displaystyle-\sqrt{2}\left(k_{c}\left|\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\right)\right]\right|-L_{d}\right)V_{c}^{\frac{1}{2}}$

From (152), if
$\left|\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\right)\right]\right|>\frac{L_{d}}{k_{c}}$,
then $\dot{V}_{c}\left(t\right)<0$, and
$\left|e_{2}\left(t\right)+e_{2c}\right|$ decreases, and we get

$\left|\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\right)\right]\right|\leq\frac{L_{d}}{k_{c}}$
(189)

i.e.,

$\left|1-\frac{2}{\text{e}^{2\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\right)}+1}\right|\leq\frac{L_{d}}{k_{c}}$
(190)

Therefore, there exists a finite time $t_{c1}>0$, for $t\geq t_{c1}$, such
that

$\left|e_{2}\left(t\right)+e_{2c}\right|\leq\frac{1}{2\rho_{c}}\ln\frac{k_{c}+L_{d}}{k_{c}-L_{d}}$
(191)

Because $\rho_{c}\gg\frac{1}{2}\ln\frac{k_{c}+L_{d}}{k_{c}-L_{d}}$, we get
$\left|e_{2}\left(t\right)+e_{2c}\right|\ll 1$. There exists
$\left|O\left(1/\rho\right)\right|\leq\frac{1}{2\rho_{c}}\ln\frac{k_{c}+L_{d}}{k_{c}-L_{d}}$,
such that $e_{2}\left(t\right)=-e_{2c}+O\left(1/\rho\right)$. From
$\dot{e}_{1}=e_{2}$, we get $\dot{e}_{1}=-e_{2c}+O\left(1/\rho_{c}\right)$ for
$t\geq t_{c1}$. Therefore, there exists a finite time $t_{c2}>t_{c1}>0$, for
$t\geq t_{c2}$, such that $e_{1}\left(t\right)\leq e_{1c}$.

We know that $k_{c}>L_{d}$. Therefore, there exists a finite time $t_{c1}>0$,
for $t\geq t_{c1}$, such that $e_{2}\left(t\right)=-e_{2c}$. According to
$\dot{e}_{1}=e_{2}$, we get that $\dot{e}_{1}\left(t\right)=-e_{2c}$ for
$t\geq t_{c1}$. Therefore, there exists a finite time $t_{c}>0$, for $t\geq
t_{c}$, such that $e_{1}\left(t\right)\leq e_{1c}$.

Case two: If $e_{1}\left(t\right)<-e_{1c}$, for (24), we get

$\frac{d\left(e_{2}\left(t\right)-e_{2c}\right)}{dt}=-k_{c}\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)-e_{2c}\right)\right]+d(t)$
(192)

Similar method to case $e_{1}\left(t\right)>e_{1c}$, when we select the
Lyapunov function candidate as
$V_{c}=\frac{1}{2}\left(e_{2}\left(t\right)-e_{2c}\right)^{2}$, there exists a
finite time $t_{c}>0$, for $t\geq t_{c}$, such that $e_{1}\left(t\right)\geq-
e_{1c}$.

Combining cases one and two, we can get that
$\left|e_{1}\left(t\right)\right|\leq e_{1c}$ for $t\geq t_{c}$.

In the following, we discuss the sliding mode system when
$\left|e_{1}\left(t\right)\right|\leq e_{1c}$.

For the smoothed sliding mode (24), a Lyapunov function candidate is selected
as

$V_{2}\left(t\right)=\frac{1}{2}\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]^{2}$
(193)

Then, taking the derivative for $V_{2}\left(t\right)$, we get

$\displaystyle\dot{V}_{2}$ $\displaystyle=$
$\displaystyle\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]\left\\{-k_{2}\text{tanh}\left[\rho\left(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right)\right]-d(t)+k_{1}e_{2}\left(t\right)\right\\}$
(194) $\displaystyle\leq$ $\displaystyle-
k_{2}\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]\left|\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]\right|\text{sign}\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]$
$\displaystyle+(k_{1}\left|e_{2}\left(t\right)\right|+L_{d})\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right|$
$\displaystyle=$
$\displaystyle-\left(k_{2}\left|\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]\right|-k_{1}\left|e_{2}\left(t\right)\right|-L_{d}\right)\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right|$
$\displaystyle=$
$\displaystyle-\left(k_{2}\left|\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]\right|-k_{1}\left|e_{2}\left(t\right)\right|-L_{d}\right)V^{\frac{1}{2}}\left(t\right)$

From (158), if
$\left|\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]\right|>\frac{k_{1}\left|e_{2}\left(t\right)\right|+L_{d}}{k_{2}}$,
then $\dot{V}_{2}\left(t\right)<0$, and
$\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right|$ decreases. We know
that, the function
$\left|\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]\right|$
is the the monotonically increasing function about
$\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right|$. Therefore,
$\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right|$ decreases until

$\left|\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]\right|\leq\frac{k_{1}\left|e_{2}\left(t\right)\right|+L_{d}}{k_{2}}$
(195)

We define

$\max\left\\{\left|e_{2}\left(t\right)\right|\right\\}=\max\left\\{\left|e_{2}\left(t_{c}\right)\right|,\frac{k_{1}}{3}\left[\left|e_{1}\left(t_{c}\right)\right|+\sqrt{e_{1}^{2}\left(t_{c}\right)+3\left(\frac{e_{2}\left(t_{c}\right)}{k_{1}}\right)^{2}}\right]\right\\}\overset{\text{define}}{=}e_{2\max}$
(196)

For (159), from (160), we get

$\left|1-\frac{2}{\text{e}^{2\rho\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]}+1}\right|\leq\frac{k_{1}e_{2\max}+L_{d}}{k_{2}}$
(197)

Function $e_{2}+k_{1}e_{1}<0$ holds for $t\in\left[t_{c},t_{c}+t_{s}\right)$
in the zones (II-2, IV-1, III-1 and III-2). Therefore, (161) can be expressed
by

$\frac{2}{\text{e}^{2\rho\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]}+1}-1\leq\frac{k_{1}e_{2\max}+L_{d}}{k_{2}}$
(198)

Then, it follows that

$\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right|\leq\frac{1}{2\rho}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}$
(199)

From the relation $\dot{e}_{1}\left(t\right)=e_{2}\left(t\right)$ in the
sliding mode (24), for (163), there exists a function $\beta\left(t\right)$,
where
$\left|\beta\left(t\right)\right|\leq\frac{1}{2\rho}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}$,
such that the following convergence law holds:

$\dot{e}_{1}\left(t\right)=-k_{1}e_{1}\left(t\right)+\beta\left(t\right)$
(200)

The solution to the convergence law (164) is

$e_{1}\left(t\right)=\left(\int_{t_{c}}^{t}\beta\left(\tau\right)\text{e}^{k_{1}\tau}d\tau\right)\text{e}^{-k_{1}t}$
(201)

Therefore, we get

$\left|e_{1}\left(t\right)\right|\leq\left|\beta\left(\tau\right)\right|\left(\int_{t_{c}}^{t}\text{e}^{k_{1}\tau}d\tau\right)\text{e}^{-k_{1}t}\leq\frac{1}{2\rho
k_{1}}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}\left(1-\text{e}^{-k_{1}\left(t-t_{c}\right)}\right)$
(202)

Then, it follows that

$\underset{t\rightarrow\infty}{\lim}\left|e_{1}\left(t\right)\right|\leq\frac{1}{2\rho
k_{1}}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}$ (203)

Because
$\rho\gg\frac{1}{2k_{1}}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}$
is selected, the up-bound of
$\underset{t\rightarrow\infty}{\lim}\left|e_{1}\left(t\right)\right|$ is
sufficiently small, and $\frac{1}{2\rho
k_{1}}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}\ll 1$
holds. For $e_{2}\left(t\right)$, from (163) and (166), we have

$\displaystyle\left|e_{2}\left(t\right)\right|$ $\displaystyle=$
$\displaystyle\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)-k_{1}e_{1}\left(t\right)\right|\leq\left|e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right|+k_{1}\left|e_{1}\left(t\right)\right|$
(204) $\displaystyle\leq$
$\displaystyle\frac{1}{2\rho}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}+\frac{1}{2\rho}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}\left(1-\text{e}^{-k_{1}\left(t-t_{c}\right)}\right)$

Therefore, we get

$\underset{t\rightarrow\infty}{\lim}\left|e_{2}\left(t\right)\right|\leq\frac{1}{\rho}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}$
(205)

Because
$\rho\gg\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}$ is
selected, the up-bound of
$\underset{t\rightarrow\infty}{\lim}\left|e_{2}\left(t\right)\right|$ is
sufficiently small, and
$\frac{1}{\rho}\ln\frac{k_{2}+k_{1}e_{2\max}+L_{d}}{k_{2}-k_{1}e_{2\max}-L_{d}}\ll
1$ holds.

(ii) Specially, when $\rho$ is selected large enough, i.e.,
$\rho\rightarrow+\infty$, the sliding variable up-bounds in (167) and (169)
approach to zero, and

$\underset{t\rightarrow\infty}{\lim}\underset{\rho\rightarrow+\infty}{\lim}e_{1}\left(t\right)=0\text{
and
}\underset{t\rightarrow\infty}{\lim}\underset{\rho\rightarrow+\infty}{\lim}e_{2}\left(t\right)=0$

From $\underset{\rho\rightarrow+\infty}{\lim}$tanh$\left(\rho\cdot
x\right)=$sign$(x)$, the tanh-function-based sliding mode (24) becomes the
ideal 2-sliding mode (18). This concludes the proof. $\blacksquare$

_Proof of Theorem 5.1:_

Define $e_{1}\left(t\right)=x_{d}\left(t\right)-x_{1}\left(t\right)$ and
$e_{2}\left(t\right)=\dot{x}_{d}\left(t\right)-x_{2}\left(t\right)$. Then, the
error system is

$\displaystyle\dot{e}_{1}\left(t\right)$ $\displaystyle=$ $\displaystyle
e_{2}\left(t\right)$ $\displaystyle\dot{e}_{2}\left(t\right)$ $\displaystyle=$
$\displaystyle-h(t)-u(t)+\ddot{x}_{d}\left(t\right)+\delta\left(t\right)$
(206)

The desired stable sliding mode (18) in Theorem 3.1 is selected, where,
$d(t)=\ddot{x}_{d}\left(t\right)+\delta\left(t\right)$. In order to turn the
error system (170) into the sliding mode (18), we select

$\displaystyle\dot{e}_{2}\left(t\right)$ $\displaystyle=$
$\displaystyle-h(t)-u(t)+\ddot{x}_{d}\left(t\right)+d\left(t\right)$ (209)
$\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l}-k_{c}\text{sign}\left[e_{2}\left(t\right)+e_{2c}\text{sign}\left(e_{1}\left(t\right)\right)\right]+\ddot{x}_{d}\left(t\right)+\delta\left(t\right),\text{
if }\left|e_{1}\left(t\right)\right|>e_{1c};\\\
-k_{2}\text{sign}\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]+\ddot{x}_{d}\left(t\right)+\delta\left(t\right),\text{
if }\left|e_{1}\left(t\right)\right|\leq e_{1c}\end{array}\right.$

Therefore, we get the controller as follows:

$u(t)=\left\\{\begin{array}[]{l}k_{c}\text{sign}\left[e_{2}\left(t\right)+e_{2c}\text{sign}\left(e_{1}\left(t\right)\right)\right]-h(t),\text{
if }\left|e_{1}\left(t\right)\right|>e_{1c}\\\
k_{2}\text{sign}\left[e_{2}\left(t\right)+k_{1}e_{1}\left(t\right)\right]-h(t),\text{
if }\left|e_{1}\left(t\right)\right|\leq e_{1c}\end{array}\right.$ (210)

Thus, for the uncertain system (39), when the controller (40) is selected, the
system is stable, and $x_{1}$ tracking $x_{d}\left(t\right)$ is non-
overshooting. This concludes the proof. $\blacksquare$

_Proof of Theorem 5.2:_

Define $e_{1}\left(t\right)=x_{d}\left(t\right)-x_{1}\left(t\right)$ and
$e_{2}\left(t\right)=\dot{x}_{d}\left(t\right)-x_{2}\left(t\right)$. Then, the
error system is

$\displaystyle\dot{e}_{1}\left(t\right)$ $\displaystyle=$ $\displaystyle
e_{2}\left(t\right)$ $\displaystyle\dot{e}_{2}\left(t\right)$ $\displaystyle=$
$\displaystyle-h(t)-u(t)+\ddot{x}_{d}\left(t\right)+\delta\left(t\right)$
(211)

The desired stable sliding mode (24) in Theorem 3.2 is selected, where,
$d(t)=\ddot{x}_{d}\left(t\right)+\delta\left(t\right)$. In order to turn the
error system (173) into the sliding mode (24), we select

$\displaystyle\dot{e}_{2}$ $\displaystyle=$
$\displaystyle-h(t)-u(t)+\ddot{x}_{d}\left(t\right)+d\left(t\right)$ (214)
$\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l}-k_{c}\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\text{sign}\left(e_{1}\left(t\right)\right)\right)\right]+\ddot{x}_{d}\left(t\right)+\delta\left(t\right),\text{
if }\left|e_{1}\left(t\right)\right|>e_{1c};\\\
-k_{2}\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]+\ddot{x}_{d}\left(t\right)+\delta\left(t\right),\text{
if }\left|e_{1}\left(t\right)\right|\leq e_{1c}\end{array}\right.$

Therefore, we get the controller as follows:

$u(t)=\left\\{\begin{array}[]{l}k_{c}\text{tanh}\left[\rho_{c}\left(e_{2}\left(t\right)+e_{2c}\text{sign}\left(e_{1}\left(t\right)\right)\right)\right]-h(t),\text{
if }\left|e_{1}\left(t\right)\right|>e_{1c};\\\
k_{2}\text{tanh}\left[\rho(e_{2}\left(t\right)+k_{1}e_{1}\left(t\right))\right]-h(t),\text{
if }\left|e_{1}\left(t\right)\right|\leq e_{1c}\end{array}\right.$ (215)

Thus, for the uncertain system (39), when the controller (44) is selected, the
system is stable, and $x_{1}$ tracking $x_{d}\left(t\right)$ is non-
overshooting. This concludes the proof. $\blacksquare$


# Model exploration in gravitational-wave astronomy with the
maximum population likelihood

Ethan Payne<EMAIL_ADDRESS>Department of Physics, California Institute of
Technology, Pasadena, California 91125, USA  LIGO Laboratory, California
Institute of Technology, Pasadena, California 91125, USA  School of Physics
and Astronomy, Monash University, VIC 3800, Australia  OzGrav: The ARC Centre
of Excellence for Gravitational-Wave Discovery, Clayton, VIC 3800, Australia
Eric Thrane<EMAIL_ADDRESS>School of Physics and Astronomy, Monash
University, VIC 3800, Australia  OzGrav: The ARC Centre of Excellence for
Gravitational-Wave Discovery, Clayton, VIC 3800, Australia

###### Abstract

Hierarchical Bayesian inference is an essential tool for studying the
population properties of compact binaries with gravitational waves. The basic
premise is to infer the unknown prior distribution of binary black hole and/or
neutron star parameters such component masses, spin vectors, and redshift.
These distributions shed light on the fate of massive stars, how and where
binaries are assembled, and the evolution of the Universe over cosmic time.
Hierarchical analyses model the binary black hole population using a prior
distribution conditioned on hyper-parameters, which are inferred from the
data. However, a misspecified model can lead to faulty astrophysical
inferences. In this paper we answer the question: given some data, which prior
distribution––from the set of all possible prior distributions––produces the
largest possible population likelihood? This distribution (which is not a true
prior) is  – $\pi$ (pronounced “pi stroke”), and the associated maximum
population likelihood is  – $\mathcal{L}$ (pronounced “L stroke”). We
postulate that the structure of  – $\pi$ is a linear superposition of delta
functions. We show how  – $\pi$ and  – $\mathcal{L}$ can be used for model
exploration/criticism. We apply this  – $\mathcal{L}$ formalism to study the
population of binary black hole mergers observed in LIGO–Virgo–KAGRA’s third
Gravitational-Wave Transient Catalog. Based on our results, we discuss
possible improvements for gravitational-wave population models.

##  I Motivation

Bayesian inference has become a mainstay of modern scientific data analysis as
a means of analysing signals in noisy observations. This procedure determines
the posterior distributions for parameters given one or more model. In order
to study the population properties of a set of uncertain observations, a
hierarchical Bayesian framework can be employed. The basic idea is to model
the population using a conditional prior $\pi(\theta|\Lambda,M)$, which
describes, for example, the distribution of black hole masses
$\\{m_{1},m_{2}\\}\in\theta$ given some hyper-parameters $\Lambda$, which
determine the shape of the prior distribution. Here, $M$ denotes the choice of
model. One then carries out Bayesian inference using a “population likelihood”

$\displaystyle{\cal
L}(d|\Lambda,M)=\prod_{i}^{N}\frac{1}{\xi(\Lambda)}\int\differential\theta_{i}\,{\cal
L}(d_{i}|\theta_{i})\pi(\theta_{i}|\Lambda,M),$ (1)

where ${\cal L}(d_{i}|\theta_{i})$ is the likelihood for data associated with
event $i$ given parameters $\theta_{i}$, and $\xi(\Lambda)$ is the detected
fraction for a choice of hyper-parameters. Meanwhile, $N$ is the total number
of observations. For an overview of hierarchical modeling in gravitational-
wave astronomy including selection effects, see Refs. Thrane and Talbot
(2019); Vitale  _et al._ (2022); Mandel  _et al._ (2019).

The LIGO-Virgo-KAGRA (LVK) Collaboration’s third gravitational-wave transient
catalog (GWTC-3) Abbott  _et al._ (2021a) contains the cumulative set of
observations of $N=69$ confident binary black-hole mergers 111We adopt the
threshold utilized in (Abbott  _et al._ , 2021c) of a false-alarm-rate
$<1\,\textrm{yr}^{-1}$. detected by the LVK (Aasi  _et al._ , 2015; Acernese
_et al._ , 2015; Akutsu  _et al._ , 2021). Additional detection candidates
have been put forward by independent groups Olsen  _et al._ (2022); Nitz  _et
al._ (2021); Zackay  _et al._ (2021); Venumadhav  _et al._ (2020); Zackay  _et
al._ (2019). Hierarchical inference is employed to study the population
properties these merging binary black holes; see, e.g., Refs. Abbott  _et al._
(2019, 2021b, 2021c); Roulet  _et al._ (2021); Farr  _et al._ (2017); Talbot
and Thrane (2018, 2017); Callister  _et al._ (2021); Fishbach  _et al._
(2022); Biscoveanu  _et al._ (2021, 2022); Vitale  _et al._ (2017); Stevenson
_et al._ (2017); Miller  _et al._ (2020); Galaudage  _et al._ (2021); Fishbach
_et al._ (2018); Edelman  _et al._ (2022a, b); Golomb and Talbot (2022a).
These analyses have revealed a number of exciting results, such as the
surprising excess rate of mergers with a primary black hole mass of $\sim
35\,M_{\odot}$ (Abbott  _et al._ , 2021b), and the evolution of the binary
merger rate with redshift (Abbott  _et al._ , 2021c), to name just two.

However, Bayesian inference has its limitations. One can use Eq. (1) in order
to infer the distribution of binary black hole parameters—given some model;
and one can compare the marginal likelihoods of two models to see which one
better describes the data. However, Bayesian inference does not tell us if any
of the models we are using are suitable descriptions of the data. While all
models for the distribution of binary black hole parameters are likely to be
imperfect, some may be adequate for describing our current dataset 222Here, we
paraphrase the aphorism attributed to statistician, George Box: “all models
are wrong, but some are useful.”. When a model fails to capture some salient
feature of the data, it is said to be “misspecified” (Romero-Shaw  _et al._ ,
2022; Gelman  _et al._ , 2013). Some effort has been made to assess the
suitability of gravitational-wave models, both qualitatively and
quantitatively; see, e.g., Abbott  _et al._ (2021b, c); Romero-Shaw  _et al._
(2022); Essick  _et al._ (2022). However, the idea of “model
criticism”—testing the suitability of Bayesian models—is still being developed
within the context of gravitational-wave astronomy and beyond.

Hierarchical Bayesian inference studies often depend upon parametric models.
Modelers design parameterizations in order to capture the key features of the
astrophysical distributions. However, one must still worry about “unknown
unknowns”—features which do not occur to the modeler to add. For example,
recent studies Abbott  _et al._ (2021b, c); Callister  _et al._ (2022); Tong
_et al._ (2022) find a sub-population of binary black holes merge with spin
vectors that are misaligned with respect to the orbital angular momentum axis.
However, the degree to which the spins are misaligned might be model
dependent. In Refs. Abbott  _et al._ (2021b, c); Callister  _et al._ (2022),
the inferred minimum spin tilt is confidently $\gtrsim 90^{\circ}$. In
contrast, Refs. Roulet  _et al._ (2021); Galaudage  _et al._ (2021); Tong  _et
al._ (2022) argue this signature could be due to a lack of flexibility in LVK
models to account for a sub-population of black holes with negligible spin
magnitude, finding support for misalignment at smaller minimum tilt angles.
The inferred population distribution of spin misalignment has important
consequences for understanding the formation channels of binary black-hole
channels. This debate highlights how astrophysical inferences can be affected
by model design.

In order to help alleviate some of the issues arising from model
misspecification in Bayesian inference, we present a framework for assessing
the suitability of a model. This framework is built around the concept of the
maximum population likelihood – $\mathcal{L}$ (pronounced “L stroke”)—the
largest possible value of ${\cal L}(d|\Lambda)$ in Eq. (1), maximized over all
possible choices of population model $\pi(\theta|\Lambda)$ independent of the
choice of parameterization. The “prior” distribution, which yields this
maximum is $\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$
(pronounced “pi stroke”). It is not a true prior because it is determined by
the data.

The  – $\mathcal{L}$ framework is useful for several reasons. First, the
numerical value of  – $\mathcal{L}$ is an upper bound on the population
likelihood. We can compare the maximum likelihood for a specific model

$\mathcal{L}_{\textrm{max}}(M)=\max_{\Lambda\sim p(\Lambda|d)}{\cal
L}(d|\Lambda,M)$ (2)

to  – $\mathcal{L}$ . Often in Bayesian model selection, the Bayesian evidence
values (${\cal Z}_{i}$) of two hypotheses can be used to determine the extent
to which one model is preferred over the other. A typical threshold chosen to
rule out one model in favor of another is that $\ln({\cal Z}_{1}/{\cal
Z}_{2})>8$ (Jeffreys, 1961). In a similar vein, if
$\ln(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}/\mathcal{L}_{\textrm{max}}(M))\lesssim
8$, we can be sure the model M is not badly misspecified since there is no
second model $M^{\prime}$ that can be written down with that will yield a
statistically significant improvement. We emphasize that a model which does
not satisfy this condition is not necessarily misspecified.

Second, the  – $\mathcal{L}$ framework can be used to quantitatively assess if
a model $M$ is misspecified. By generating synthetic data from $M$, one can
generate the expected distribution of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$.
In this paper, we show how one can compare the observed values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
to the expected distribution in order to determine the extent to which $M$ is
misspecified—and the way in which it is misspecified.

Third, the  – $\mathcal{L}$ framework can be used for “model
exploration”—providing clues of where in parameter space unmodeled features
might be lurking. By comparing
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ with the prior
from our phenomenological model $\pi(\theta|M)$, one can see if the
phenomenological model is capturing key structure present in  – $\pi$ and use
the comparison to design new models to test on forthcoming datasets.

The remainder of this paper is organized as follows. In Sec. II, we introduce
the  – $\mathcal{L}$ formalism, illustrating key features with a simple toy
model. In Sec. III, we show how the formalism can be used for model criticism.
In Sec. IV, we apply the formalism to study the population properties of
merging binary black holes observed by the LVK. Our concluding remarks are
presented in Sec. V.

##  II The maximum population likelihood  – $\mathcal{L}$

###  II.1 Preliminaries

We begin with a brief review of Bayesian hierarchical inference with a
parametric model. Our starting point is the population likelihood (copied here
from Eq. (1)):

${\cal
L}(d|\Lambda,M)=\prod_{i}^{N}\frac{1}{\xi(\Lambda)}\int\differential\theta_{i}\,{\cal
L}(d_{i}|\theta_{i})\pi(\theta_{i}|\Lambda,M).$ (3)

Here, $\mathcal{L}(d_{i}|\theta_{i})$ is the likelihood of event-$i$ data
$d_{i}$ given parameters $\theta_{i}$. The quantity
$\pi(\theta_{i}|\Lambda,M)$ is a conditional prior for $\theta_{i}$ given
hyper-parameters for some population model $M$, which describes the shape of
the prior distribution. The term $\xi(\Lambda)$ accounts for selection
effects; for example, high-mass systems are typically easier to detect than
low-mass systems. It is the detectable fraction of the population given the
model given hyper-parameters $\Lambda$

$\xi(\Lambda)=\int\differential\theta\,p_{\textrm{det}}(\theta)\pi(\theta|\Lambda,M).$
(4)

Here, $p_{\textrm{det}}(\theta)$ is the detection probability of an
observation with parameters $\theta$.

###  II.2 The maximum population likelihood  – $\mathcal{L}$

The maximum population likelihood  – $\mathcal{L}$ is obtained by taking Eq.
(3) and maximizing over all possible prior distributions $\pi(\theta)$. Thus,
– $\mathcal{L}$ is an upper bound (or supremum) on the set of likelihoods from
all possible choices of models for $\pi(\theta)$ such that

$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}\equiv\mathcal{L}({d}|\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})\geq\mathcal{L}({d}|\Lambda,M),$
(5)

for all models $M$. The “prior” distribution that yields  – $\mathcal{L}$ is
denoted

$\displaystyle\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$
(8)

(pronounced “pi stroke”). It is not a true prior because the distribution
which maximizes the population likelihood in Eq. (3) depends on the data. One
should therefore refer to  – $\pi$ as a pseudo-prior. The associated model is
denoted  – $M$ (pronounced “M stroke”). Combining this notation into a single
equation, we have

$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}\equiv\prod_{i=1}^{N}\frac{1}{\xi(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})}\int\differential\theta_{i}\,\mathcal{L}(d_{i}|\theta_{i})\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta_{i}).$
(9)

###  II.3 Calculating  – $\pi$ : special cases

Having introduced the concept of  – $\mathcal{L}$ and  – $\pi$ , the natural
next question is: given data $d$, how does one calculate these quantities?
Before answering this question, we study three special cases where we can work
out  – $\pi$ from intuition. This discussion will help sharpen our instincts
for the more general solution that follows. Readers looking for to skip to the
punchline may wish to skip this subsection.

####  II.3.1 A single measurement

For the first case, we consider a single measurement ($N=1$) with a unimodal
likelihood function ${\cal L}(d|\theta)$, which is maximal when the parameter
$\theta$ is equal to the maximum likelihood value $\widehat{\theta}$. For the
sake of simplicity, we ignore selection effects so that
$\xi(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})=1$. In this case,  –
$\mathcal{L}$ in Eq. (9) is clearly maximized if the prior support is entirely
concentrated at $\widehat{\theta}$. Thus,  – $\pi$ is a delta function

$\displaystyle\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)=\delta(\theta-\widehat{\theta}),$
(12)

which yields

$\displaystyle\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}=$
$\displaystyle\int\differential\theta\,{\cal
L}(d|\theta)\,\delta(\theta-\widehat{\theta})$ (15) $\displaystyle=$
$\displaystyle{\cal L}(d|\widehat{\theta}).$ (16)

This result is intuitive: the prior that maximizes the population likelihood
is the one that concentrates all its support at the maximum-likelihood value
of $\theta$.

####  II.3.2 $N$ signals in the high-SNR Limit

For the second case, we consider a scenario in which the data consists of $N$
observations carried out in the high-SNR limit. In this limit, the likelihood
of the data for each measurement $d_{i}$ given some parameter $\theta$
approaches a delta function

$\displaystyle{\cal
L}(d_{i}|\theta_{i})=\delta(\theta_{i}-\widehat{\theta}_{i}),$ (17)

located at the maximum-likelihood value $\widehat{\theta}_{i}$. We assume that
each measurement is distinct so that no two maximum-likelihood values
$\widehat{\theta}_{i}$ are exactly the same. Again, for the sake of
simplicity, we ignore selection effects so that
$\xi(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})=1$, though, the
argument here holds even if we relax this assumption. Equation (9) becomes

$\displaystyle\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}=\prod_{i=1}^{N}\int\differential\theta_{i}\,\delta(\theta_{i}-\widehat{\theta}_{i})\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta).$
(22)

The population likelihood is maximized when  – $\pi$ is a sum of delta
functions peaking at the set of $\\{\widehat{\theta}_{i}\\}$:

$\displaystyle\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)=$
$\displaystyle\sum_{k=1}^{N}w_{k}\,\delta(\theta-\widehat{\theta}_{k})$ (25)
$\displaystyle w_{k}=$ $\displaystyle 1/N.$ (26)

This solution for  – $\pi$ ensures that there is maximal prior support at
every likelihood peak. Obviously, the population likelihood is not maximized
if any prior probability density is wasted to values of $\theta$ where all the
likelihood functions are zero. Choosing an equal weight for each delta
function $w_{i}=1/N$ produces the largest possible population likelihood.

We illustrate this case in Fig. 1(a) using high-SNR, toy-model data drawn from
a mean-zero, unit-variance Gaussian distribution. In the top panel, we plot
the set of $N=10$ maximum likelihood points $\\{\widehat{\theta}_{i}\\}$. In
the lower panel, we “plot” the
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ for these ten
data points. We put the word “plot” in quotation marks because, technically,
we are not plotting
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$, which goes to
infinity, but rather we are plotting the weights $w_{k}$ (Eq. (26)), which
allows us to see the relative weight given to each delta function—something
that will prove useful below. Throughout the paper, when we refer to plots of
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$, it should be
understood that we are actually plotting representations of
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ using the
weights $w_{k}$. Finally, note that each peak in the distribution of
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ matches up
with one of the maximum likelihood points in the upper panel.

Figure 1:  Examples of the distribution
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ described in
Subsections II.2-II.4. Each column represents a different dataset. The top-
panel dots show the set of $N=10$ maximum-likelihood estimates
$\\{\widehat{\theta}_{i}\\}$. The top-panel horizontal lines represent error
bars; (in the first column they are too small to see). The bottom panels show
the distribution of
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ associated
with each data set. The left-hand column (a) represents data in the high-SNR
limit so that the likelihood functions for each measurement approach delta
functions (this is why the error bars are not visible). In this case,
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ consists of
$N$ delta functions, each associated with one of the maximum likelihood points
$\widehat{\theta}_{i}$. In the middle column (b), we are no longer in the
high-SNR limit, but the maximum likelihood points are all assumed to be
identical with $\widehat{\theta}_{i}=0$. In this case,
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ consists of
one delta function peaking at $\theta=0$. In the right-hand column (c), the
data are not in the high-SNR limit, and each $\widehat{\theta}_{i}$ is random.
In this case, $\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$
consists of $n=3$ delta functions, each with different heights.

####  II.3.3 $N$ identical measurements

For the third case, we consider a set of $N$ observations. This time, we do
not assume the high-SNR limit, but we assume that every measurement has the
same maximum-likelihood value of $\widehat{\theta}$. This case is highly
contrived—one does not typically work with multiple identical measurements—but
the example is nonetheless helpful for illustrative purposes. In this case,
the integral in Eq. (9) is maximized when the prior support is entirely
concentrated at $\widehat{\theta}$ (where all of the likelihood functions
peak), so that  – $\pi$ is a single delta function:

$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)=\delta(\theta-\widehat{\theta}),$
(27)

while

$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}=\prod_{i=1}^{N}\mathcal{L}(d_{i}|\widehat{\theta}).$
(28)

This scenario is demonstrated in Fig. 1(b). The top panel shows the set of
$N=10$ maximum-likelihood points $\\{\widehat{\theta}_{i}\\}$, all with the
same value. The horizontal lines represent the error bars for each
measurement, which we draw from a uniform distribution on the interval
$(0.01,1)$. In the lower panel, we plot
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ these ten data
points. This time, since every measurement is identical,
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ is a single
delta function peaking at $\theta=0$.

From these three examples, we observe a pattern: in each case,
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ can be written
as a weighted sum of delta functions. Numerical studies further lead us to
believe that $\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$
is always given by a sum of delta functions. We explore this hypothesis in the
next subsection.

###  II.4 The general form of  – $\pi$

In the previous subsection, we studied three special cases where the solution
to Eq. (9) was a weighted sum of delta functions. Based on numerical studies,
we believe this result is true in general, regardless of the form of the
likelihood ${\cal L}(d|\theta)$ and the selection effect term
$p_{\text{det}}(\theta)$. We therefore conjecture that, for any set of
observations, $\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$
is always of the form,

$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)=\sum_{k=1}^{n}w_{k}\,\delta(\theta-\theta_{k}),$
(29)

where $w_{k}$ are weights which sum to unity

$\sum_{k=1}^{n}w_{k}=1.$ (30)

We conjecture that the number of delta function is always less than or equal
to the number of measurements in all but the most pathological of cases (e.g.,
multimodal distributions with equal likelihoods) so that

$\displaystyle n\leq N.$ (31)

The ratio

$\displaystyle\mathcal{I}\equiv n/N,$ (32)

is a measure of the “informativeness” of the data. It compares the typical
likelihood width to the scatter in the astrophysical distribution. In the
high-SNR limit, $\mathcal{I}=1$ since a delta function is required for every
data point (see Fig. 1(a)). The other limiting case is, $\mathcal{I}=1/N$,
which happens when the likelihood for each measurement overlaps (see Fig.
1(b)).

Using this conjecture, we now consider a variation on the toy-model problems
discussed in the earlier subsections. In particular, we consider finite-SNR
data drawn from our Gaussian, toy-model distribution. Using Eqs. (29-30) as an
ansatz, we calculate
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ for $N=10$
random data points. The maximum likelihood values $\widehat{\theta}_{i}$ are
drawn from a mean-zero, unit-variance Gaussian and the error bars are drawn
from a uniform distribution on the interval $(0.01,1)$. The results of this
calculation are shown in Fig. 1(c). The top panel shows the data, represented
by the maximum-likelihood values $\\{\widehat{\theta}_{i}\\}$, which are
arranged from bottom to top in increasing order. The horizontal lines show the
uncertainty for each measurement. In the bottom panel, we show
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ for this
dataset. It consists of just $n=3$ delta functions of varying heights
($\mathcal{I}=0.3$). The exact weights, locations, and number of delta
functions are not obvious; we obtain them numerically by maximising Eq. (29)
subject to Eq. (30) using the “combined” method described below in Subsection
II.5. Comparing the red data points with error bars to the turquoise
representation of
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$, one can see
that every data point can be plausibly associated with at least one of the
delta functions.

Assuming the validity of our conjecture, we can write down a general
expression for  – $\mathcal{L}$ :

$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}=\prod_{i=1}^{N}\frac{1}{\xi(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})}\sum_{k=1}^{n}w_{k}\,\mathcal{L}(d_{i}|\theta_{k}),$
(33)

where

$\xi(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})=\sum_{k=1}^{n}w_{k}\,p_{\textrm{det}}(\theta_{k}).$
(34)

Given Eqs. (33) and (34), the problem of calculating
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}$
reduces to the problem of simply finding the locations and weights of $n$
delta functions. In Section II.5, we explore three different approaches to
this problem. One of these algorithms will allow us to test our conjecture and
show that it is born out numerically.

###  II.5 Computing  – $\pi$

In this subsection, we consider three techniques that can be applied to
compute
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}$:
optimization, iterative grid, and stochastic methods. We show that a combined
approach, which uses a grid-based approach to guess a solution, which is
subsequently refined through optimization performs the best out of the
algorithms we tried. Meanwhile, the stochastic approach allows us to test our
conjecture with minimal assumptions.

Figure 2:  Demonstration of different methods for calculating
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}},\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}$.
Each panel shows the results for a different number of measurements with
$N=10$ on top, $N=100$ in the middle, and $N=1000$ on the bottom. The black
distribution is the true distribution $\pi(\theta)$ used to generate the data.
The colored spikes show the reconstructed distribution
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ as determined
by different methods. Cyan is for the “combined” technique, which uses the
iterative grid to obtain a first guess that is refined with the optimization
method. Meanwhile, orange is for the grid-based technique by itself and gray
is for the stochastic method.

####  II.5.1 Optimization

The first approach we consider is to use an optimization algorithm subject to
the constraint in Eq. (30). We use Scipy’s trust-constr optimization
implementation (Virtanen  _et al._ , 2020; Conn  _et al._ , 2000). We find
this approach fails to find the correct global maximum of Eq. (33) once the
number of peaks $n$ becomes large. However, this issue can be resolved if a
sufficiently close guess to the true shape of
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ can be made.
Fortunately, the iterative-grid approach can be used to supply this initial
guess.

####  II.5.2 Iterative grid

The second approach we consider is to iteratively place delta functions on a
fixed grid. There are two steps: the greedy addition of many delta functions,
and the removal of no-longer-useful delta functions. In the first step, we
first attempt to place a delta function with a fixed height at each grid point
and evaluate Eq. (33) (with appropriate normalization of the distribution). We
determine which of all possible delta function additions produces the highest
population likelihood. We then vary the height of this delta function between
zero and twice initial height in order to obtain an updated guess for
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$. The addition
of delta functions is repeated, reducing the initial height by a factor at
each iteration. After many iterations, we then attempt to remove no-longer-
useful delta functions to further increase the population likelihood. We
repeat this procedure five times, iteratively adding 30 delta functions with
varying heights at each iteration. After these iterations,  – $\mathcal{L}$ is
usually well-converged for the problems we are studying. In some iterations,
this procedure adds support to preexisting delta functions. This is how the
approach “corrects” under-supported delta functions.

This method has a significant advantage over generic constrained optimization
techniques as the procedure does not require the optimization of individual
parameters governing the delta functions through the $\\{\theta_{k},w_{k}\\}$
space. However, we find that this method is improved by pairing it with
optimization. The most accurate optimization of the maximum population
likelihood and structure of the distribution occurs when we utilize grid-based
approximation to inform the starting location and weights for the constrained
optimization. This allows for the grid-based approximation to find the region
of parameter space where  – $\mathcal{L}$ is nearly maximal. The constrained
optimization then purifies the delta function structure and slightly increases
the maximum population likelihood. The combined method is used for all the
maximum population likelihood computations in Sec. IV.

####  II.5.3 Stochastic construction

Our final approach is to stochastically generate samples for
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$, which are
accepted/rejected depending on whether the new samples increases the
population likelihood. This is a form of importance sampling in which an
arbitrary “proposal distribution” is used to generate proposal samples. When a
proposal sample is generated, we add it to a list of previously accepted
points and evaluate  – $\mathcal{L}$ as a Monte Carlo integral,

$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}=\prod_{i=1}^{N}\frac{1}{\xi(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})}\Big{\langle}\mathcal{L}(d_{i}|\theta_{i})\Big{\rangle}_{\theta_{i}\sim\text{\ooalign{\raisebox{-0.60275pt}{--}\cr$\pi$\cr}}(\theta_{i})},$
(35)

where

$\xi(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$M$\cr}})=\Big{\langle}p_{\textrm{det}}(\theta)\Big{\rangle}_{\theta\sim\text{\ooalign{\raisebox{-0.60275pt}{--}\cr$\pi$\cr}}(\theta)}.$
(36)

Here, the angled brackets indicate averaging over the samples. If the addition
of the new sample increases  – $\mathcal{L}$ , we retain the sample in the
list of samples from  – $\pi$ . As the process is repeated, the set of samples
produces an ever-improving representation of  – $\pi$ .

This method can be extended to employ an additional burn-in phase and/or a
thinning phase to ensure more rapid convergence by removing unfavorable
samples that sometimes get accepted early on before the distribution is well-
converged. While this approach converges more slowly than the other two
methods, it does not employ any assumptions about the structure of the
distribution. Thus, this method can be used to validate the conjecture put
forward in Eqs. (33-34), that
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ is a sum of
delta functions.

####  II.5.4 Numerical study

We demonstrate each method using our Gaussian, toy-model distribution
described in the last subsection: true maximum likelihood values
$\widehat{\theta}_{i}$ drawn from zero-mean, unit-variance Gaussian with error
bars drawn from a uniform distribution on the interval (0.01, 1). The observed
maximum likelihood values are then shifted from the true value by an offset
generated from each individual observation’s uncertainty. The results of this
demonstration are compiled in Fig. 2. The three panels of Fig. 2 represent
tests performed with $N=10$, $100$, and $1000$ observations. In each panel,
the black curve represents the true distribution $\pi(\theta)$. The colored
spikes illustrate different numerical solutions for
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$: cyan is the
“combined” approach, which uses the iterative grid to obtain an initial guess
that is subsequently refined using the optimization method. Meanwhile, orange
represents the iterative grid approach by itself. For the grid-based approach
we run 30 iterations of adding peaks with variable but decreasing weights,
before repeating this process an additional ten times. Finally, gray
represents the stochastic approach. For the stochastic method, we generate
3000 samples with 1000 samples for burn-in.

We see that the combined approach better estimates  – $\mathcal{L}$ relative
to the other techniques considered 333A method is “better” if it yields a
larger value of  – $\mathcal{L}$ than another approach.. We observe that, as
$N$ increases, $\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$
increasingly resembles the true Gaussian distribution $\pi(\theta)$ (shown in
Fig. 2 as a black curve). To illustrate this more clearly, we take the
inferred delta function locations from the $N=1000$ “combined” result in Fig.
2 and compute the weighted histogram. This result is directly compared to the
true distribution in Fig. 3, from which we see that indeed the inferred
distribution is (albeit slowly) approaching the true distribution. We
conjecture that, in general,
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ approaches the
true distribution in the infinite-data limit:

$\displaystyle\lim_{N\rightarrow\infty}\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)\rightarrow\pi_{\textrm{true}}(\theta).$
(39)

Figure 3:  Comparison between a binned representation of  – $\pi$ as computed
for the toy model data set with $N=1000$ observations and the true underlying
population distribution. This representation more clearly shows that  – $\pi$
is approaching the true distribution in the limit of many observations.

####  II.5.5 Computational challenges

Before continuing, we discuss two computational challenges. First, we note
that the examples illustrative above are all one-dimensional. The discussion
above generalizes to $\geq 2$ dimensions; we conjecture that
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ is still a sum
of delta functions in $\geq 2$ dimensions. However, it becomes increasingly
challenging to determine the location and height of these peaks in higher
dimensions. It seems plausible that some version of the stochastic method
becomes relatively more efficient when there are many dimensions, though,
additional work is required to test this hypothesis. Furthermore, by
increasing the dimensionality of the problem, constructing continuous
representations of the individual event likelihoods and the detection
probability, $p_{\textrm{det}}(\theta)$, becomes increasing difficult. Recent
developments in using Gaussian mixture models to produce continuous
representations of these distributions might alleviate these concerns Talbot
and Thrane (2022); Golomb and Talbot (2022b). Second, even if we stay in one
dimension, the computational cost of calculating
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}},\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}$
grows with $N$. This is evident from Fig. 2: with more data comes more delta
functions (each with a location and a height), which means the solution is
characterized by an increasingly large list of numbers.

##  III Model criticism with  – $\mathcal{L}$

In this section, we show how the  – $\mathcal{L}$ formalism can be used to
determine if a model $M$ is an adequate description of data. The first step is
to generate synthetic datasets based on the posterior distribution for the
model hyper-parameters $p(\Lambda|d)$. For each data set, we calculate the
maximum population likelihood  – $\mathcal{L}$ (Eq. (9)) as well as the
maximum likelihood for $M$, which we denote

$\displaystyle\mathcal{L}_{\textrm{max}}(M)=\max_{\Lambda\sim
p(\Lambda|d)}{\cal L}(d|\Lambda,M)$ (40)

where $\mathcal{L}(d|\Lambda,M)$ is the population likelihood defined in Eq.
(3). In this way we can estimate

$\displaystyle
p(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M)),$
(43)

the joint distribution for  – $\mathcal{L}$ and
$\mathcal{L}_{\textrm{max}}(M)$ given model $M$. By comparing the measured
values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$,
to this distribution of expected values, one can see if the dataset is typical
of what one would expect given $M$. If the measured values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
are atypical, one can conclude that $M$ is misspecified. Moreover, one may
determine the nature of the misspecification by noting the location of the
observed value of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
relative to the typical values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$.
This is best illustrated with an example.

We first generate 5000 simulated datasets from our toy-model with $N=100$
simulated observations in each dataset, which consists of a Gaussian
distribution with mean $\mu=0$ and width $\sigma=1$. Next, we consider five
different underlying population distributions. Each observed population is
assigned a color:

  * •

Black: model is correctly specified $(\mu=0,\sigma=1)$.

  * •

Purple: model is too wide $(\mu=0,\sigma=0.6)$.

  * •

Blue: model is too narrow $(\mu=0,\sigma=1.4)$.

  * •

Salmon: model is shifted to one side $(\mu=1,\sigma=1)$.

  * •

Yellow: model is too wide and shifted to one side $(\mu=0.8,\sigma=0.6)$.

For each population we generate ten different realizations and compute their
values of
$\\{(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))\\}$,
again with $N=100$ observations. This demonstrates the scenario where a
particular choice in model is made, and the true population distribution may
differ.

The results are shown in Fig. 4. The vertical axis is
$\ln\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}$ while the
horizontal axis is $\ln\mathcal{L}_{\textrm{max}}(M)$. The gray region in the
bottom-right corner is forbidden since
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}\geq\mathcal{L}_{\textrm{max}}(M)$
by definition. The grey contours show the one, two, and three-sigma contours
for our toy-model with $\mu=0$ and $\sigma=1$. The colored dots, meanwhile,
show ten random realizations of each population analyzed with the toy-model.
When the toy-model is sufficiently misspecified with respect to the true
population, the observed population’s values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
typically do not reside within the three-sigma confidence interval of the
expected distribution if the toy-model was well-specified. Moreover, the
different colors of dots all cluster in different regions. For example, when
the toy-model is too broad (purple), the dots cluster above-right of the
correctly specified distribution. When the model is shifted away from the
correct peak (salmon), the dots cluster to the left of the correctly specified
distribution. Finally, if the model is well-specified (black) then the
observed population clusters within the expected distribution for the toy-
model. By studying where the observed values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
fall on this diagram, one can gain some insight into the way in which the
model is misspecified.

Figure 4:  An illustration of model criticism with the  – $\mathcal{L}$
formalism. In the left-hand panel, we plot
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
for five different populations (each with ten different realizations),
analyzed a toy-model with a mean of $\mu=0$ and standard deviation $\sigma=1$.
Each population is represented by a different color. The gray contours show
the 1, 2, and 3-sigma credible intervals for the expected distribution of
$p(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
from the toy-model. By comparing the measured values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
from an observed population to the expected distribution, one may determine if
the dataset is typical of what one would expect given the model. If the
measured values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
fall outside these intervals, one may conclude that the toy-model is
misspecified (does not accurately model the data). Moreover, the location of a
point on this plot relative to the expected distribution, conveys information
about the way in which a model is misspecified. The right-hand panel shows the
toy-model (grey), the true population distribution for the starred and labeled
datapoint (a-d), and the respective  – $\pi$ for the observed data
(turquoise). This demonstrates that shifts away from the expected distribution
(left-hand panel; grey) in
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
can be visually identifiable to the reconstruction of  – $\pi$ .

In practice, it may be computationally challenging to create plots like Fig. 4
for population studies in gravitational-wave astronomy. While it is easy to
create mock datasets, it is time-consuming to calculate individual event
likelihoods for one dataset, let alone thousands. There may be workarounds. We
discuss this possibility in greater detail below.

##  IV Application to gravitational-wave astronomy

In this section, we apply the  – $\mathcal{L}$ formalism to results from
gravitational-wave astronomy to stress-test models for the population of
merging binary black holes. We analyze data from the second gravitational-wave
transient catalog (GWTC-3) Abbott  _et al._ (2021a); LIGO Scientific
Collaboration and Virgo Collaboration and KAGRA Collaboration (2021a), which
includes 69 confidently detected binary black hole mergers with false alarm
rates $<1\,\mathrm{yr^{-1}}$. To ensure similarity to the GWTC-3 LVK
population analysis Abbott  _et al._ (2021c); LIGO Scientific Collaboration
and Virgo Collaboration and KAGRA Collaboration (2021b), we utilize the same
individual event posterior samples — constructed from equally weighted samples
generated from effective-one-body (SEOBNRv3 (Pan  _et al._ , 2014; Taracchini
_et al._ , 2014), SEOBNRv4PHM (Bohé  _et al._ , 2017; Ossokine  _et al._ ,
2020)) and phenomenological (IMRPhenomPv2 (Hannam  _et al._ , 2014),
IMRPhenomXPHM (Pratten  _et al._ , 2021)) waveform results (see Abbott  _et
al._ (2021c) for more details). We calculate the hyperparmeter distributions
and $\mathcal{L}_{\textrm{max}}(M)$ using GWPopulation Talbot  _et al._
(2019), which employs Bilby Ashton  _et al._ (2019); Romero-Shaw  _et al._
(2020) and Dynesty (Speagle, 2020). We utilize the combined injection set from
Ref. LIGO Scientific Collaboration and Virgo Collaboration and KAGRA
Collaboration (2021c) to compute the estimated detectable fraction of binary
black-hole mergers over the first three observing runs. To construct the one-
dimensional individual event likelihoods, we marginalize over the remaining
parameters assuming the maximum a posteriori hyperparameters for the
population distribution.

###  IV.1 Model inspiration through visual inspection

One straightforward application of the  – $\mathcal{L}$ formalism is to
visually compare the reconstructed population distribution (obtained using a
phenomenological model) with
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$. By comparing
these two distributions, it is possible to see which features in the
phenomenological model reconstruction are due to prior assumptions, which
features are due to real trends in the data, and which features might be
missing from the phenomenological model. Formally, we compare
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ to the
population predictive distribution (PPD)

$\displaystyle\text{PPD}(\theta|d,M)=\int\differential\Lambda\,p(\Lambda|d)\pi(\theta|\Lambda,M),$
(44)

which describes the astrophysical distribution of $\theta$ given a
phenomenological model $M$ with hyper-parameters $\Lambda$. Due to the
complexity involved in computing  – $\pi$ , we marginalize all parameters not
involved in the desired  – $\pi$ calculation. For these “nuisance” parameters,
we chose the distributions associated with the maximum a posteriori
hyperparameters from the LVK’s GWTC-3 population analysis with the Power
Law+Peak-Default-Power Law model (Abbott  _et al._ , 2021c).

In Fig. 5, we present
$\text{\ooalign{\raisebox{-0.86108pt}{--}\cr$\pi$\cr}}(\theta)$ with the PPDs
from the LVK analysis of GWTC-3 Abbott  _et al._ (2021c); LIGO Scientific
Collaboration and Virgo Collaboration and KAGRA Collaboration (2021b) for
source-frame primary mass $m_{1}$ (top), the effective inspiral spin parameter
$\chi_{\textrm{eff}}$ (middle), and redshift $z$ (bottom). Each row contains
two sub-panels; the small upper panel shows the maximum-likelihood estimate
for each gravitational-wave event and the $90\%$ confidence interval while the
larger lower panel compares  – $\pi$ with the PPD. The PPD is plotted as a
thick band to show the 90% credibility region at each value of $\theta$.

We first turn our attention to the primary mass distribution in the top row.
There are $M=10$ delta function peaks, implying an informativeness of
$\mathcal{I}=0.15$ (see Eq. (32)). The gray band is the Power Law + Peak model
from Talbot and Thrane (2018) while the orange band is a (more flexible) semi-
parametric power-law-spline model denoted Spline from Edelman  _et al._
(2022a). The agreement between  – $\pi$ and the two PPDs is striking, with
cyan spikes closely matching several of the features in both models including
the turn-over at low masses near $\approx 12M_{\odot}$ and the bump at
$30M_{\odot}$. Furthermore, we see that  – $\pi$ also recovers some of the
finer detail features found only by the Spline model. In particular, the shift
in the low-mass peak and the dips in posterior support at $\sim 16\,M_{\odot}$
and $\sim 25\,M_{\odot}$ are present in the structure of  – $\pi$ . Based on
our visual inspection, it appears that current models are capturing much if
not all of the structure present in  – $\pi$ .

Figure 5:  Population predictive distributions ($90\%$ credibility) and  –
$\pi$ for primary the black-hole mass ($m_{1}$; top), effective inspiral
parameter ($\chi_{\textrm{eff}}$; middle), and redshift ($z$; bottom)
distributions. For the redshift, we divide by the evolution of the comoving
volume and time delay as a function of redshift to plot the merger rate,
$\mathcal{R}(z)$. Comparison of the different models with  – $\pi$ highlights
which features are present in the data and which are due to assumptions in the
model.

Turning our attention to the middle row, we study the distribution of
effective inspiral spin parameter Damour (2001),

$\displaystyle\chi_{\textrm{eff}}\equiv\frac{\chi_{1}\cos\theta_{1}+q\chi_{2}\cos\theta_{2}}{1+q},$
(45)

which measures the mass-weighted black hole spin projected along the orbital
angular momentum 444In Eq. (45), $q\equiv m_{2}/m_{1}$ is mass ratio,
$\chi_{1,2}$ are the dimensionless black hole spins, and $\theta_{1,2}$ are
the spin vector tilt angles relative to the orbital angular momentum.. This
time, only $n=4$ delta function spikes are required to fit the data
($\mathcal{I}=0.06$), showing how much harder it is to measure
$\chi_{\textrm{eff}}$ than $m_{1}$. In gray, we plot the PPD for the Default
model from Refs. Abbott  _et al._ (2021b, c), which draws on work from Refs.
Talbot and Thrane (2017); Wysocki  _et al._ (2019). In orange we plot the PPD
for the Extended model from Refs. Galaudage  _et al._ (2021); Tong  _et al._
(2022), which only analyse 68 binary black-hole events in the population due
to data quality concerns regarding one event (Payne  _et al._ , 2022). To plot
the Extended model results, which incorporates a delta function at
$\chi_{\textrm{eff}}=0$, we plot the $90\%$ interval for the delta function
height, $\delta$, multiplied by the same scale factor as  – $\pi$ . The
continuous contribution to the Extended model is then scaled by the ratio of
the PPD evaluated at only the non-zero $\chi_{\textrm{eff}}$ – $\pi$ delta
functions to the previously computed scaling.

The data-driven  – $\pi$ includes a delta function at
$\chi_{\textrm{eff}}\approx 0$ and three smaller peaks in the
$\chi_{\textrm{eff}}>0$ region, but no peaks with $\chi_{\textrm{eff}}<0$. The
lack of support for $\chi_{\textrm{eff}}<0$ is in contrast to Refs Abbott  _et
al._ (2021b, c), which find support for a sub-population of binary black holes
with $\chi_{\textrm{eff}}<0$. The strong delta function at
$\chi_{\textrm{eff}}=0$ lends support to the argument put forward in Refs.
Miller  _et al._ (2020); Roulet  _et al._ (2021); Galaudage  _et al._ (2021)
that the data can be well-modeled with a sub-population of “non-spinning”
$\chi_{\textrm{eff}}=0$ binaries, even if there is not strong statistical
support for the existence of such a peak Callister  _et al._ (2022); Mould
_et al._ (2022); Tong  _et al._ (2022). However, our visual comparison
suggests that the Extended model may over-predict the abundance of binaries
with $\chi_{\textrm{eff}}\approx 0.3$.

Finally, turning our attention to the bottom row of Fig. 5, we consider the
case of redshift. For this parameter, $n=6$ ($\mathcal{I}=0.09$). Here we plot
the merger rate as a function of redshift, $\mathcal{R}(z)$ by dividing the
posterior predictive distribution by the PPD by the evolution of the comoving
volume and time delay with respect to redshift. The merger rate is more
commonly utilized for interpreting the redshift evolution. The  – $\pi$
distribution fits a decrease in the merger rate at a redshift of $z\sim 0.13$.
While we caution that  – $\pi$ is purely data-informed, and such a feature
might diminish with additional observations, the Power Law model utilized in
Refs. Abbott  _et al._ (2021b, c) lacks the flexibility to resolve such a
feature (see also Edelman  _et al._ (2022b)). As more gravitational-wave
observations are made, finer structure may emerge in the redshift evolution of
the binary merger rate. These differences between the parametric
reconstructions and  – $\pi$ might present the first hints of such structure.
We suggest that future redshift models include additional flexibility to study
the possibility of a deficit of mergers in the nearby Universe.

###  IV.2 Upper bounds on population model likelihoods

In Table 1 we report the difference in natural log likelihood comparing the
various population models to the maximum population likelihood  –
$\mathcal{L}$ :

$\displaystyle\ln\text{\ooalign{\raisebox{2.15277pt}{--}\cr$\mathcal{B}$\cr}}\equiv\ln\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}}-\ln\mathcal{L}_{\textrm{max}}(M).$
(50)

The $\ln\text{\ooalign{\raisebox{2.15277pt}{--}\cr$\mathcal{B}$\cr}}$ values
in Table 1 measure the fit of population models relative to the best possible
fit. Motivated by the typical threshold for model selection in terms of Bayes
factors (Jeffreys, 1961), a value of
$\ln\text{\ooalign{\raisebox{2.15277pt}{--}\cr$\mathcal{B}$\cr}}\lesssim 8$
indicates that the population model is very close to the maximum population
likelihood Thrane and Talbot (2019), which would imply that the fit cannot be
dramatically improved. A large value of
$\ln\text{\ooalign{\raisebox{2.15277pt}{--}\cr$\mathcal{B}$\cr}}$ by itself
does not imply that a model is “wrong” or unsuitable to describe the data, but
it does quantify the extent to which an alternative model can in-principle
improve over the current offerings.

Returning to Table 1, the Power Law + Peak model for $m_{1}$ shows the most
potential room for improvement. This may be due to structure identified using
the Spline model, which is missing from the less flexible Power Law + Peak.
However, the $m_{1}$ measurements are also the most informative in Table 1
(with the largest value of $\mathcal{I}$). With more information, it is
probably easier to concoct an a posteriori model with a large population
likelihood that explains various features in the distribution of $m_{1}$
through over-fitting. The Default and Extended spin models both exhibit
$\ln\text{\ooalign{\raisebox{2.15277pt}{--}\cr$\mathcal{B}$\cr}}<8$, which
implies that neither model can be unequivocally ruled out, though, the
Extended model provides a somewhat better fit with a natural log likelihood
difference of $4.17$. We also note that the $\chi_{\textrm{eff}}$ and $z$
observations are noticeably less informative, and simultaneously the
associated values of $\mathcal{L}_{\textrm{max}}(M)$ are closer to  –
$\mathcal{L}$ . This might indicate that, while there are features present in
– $\pi$ that are present in the data, they are not statistically significant.

Parameter | $\mathcal{I}$ | Model | $\ln\text{\ooalign{\raisebox{2.15277pt}{--}\cr$\mathcal{B}$\cr}}$
---|---|---|---
$m_{1}$ | 0.15 | Power Law + Peak | $14.89$
|  | Spline | 6.66
$\chi_{\textrm{eff}}$ | 0.06 | Default | $7.70$
|  | Extended | 3.53
$z$ | 0.09 | Power Law | $8.93$
Table 1:  The performance of different population models relative to the  –
$M$ . The quantity  – $\mathcal{B}$ (Eq. (50)) is a measure of the population
likelihood of each model relative the maximum possible population likelihood
– $\mathcal{L}$ . The “informativeness” $\mathcal{I}$ (Eq. (32)) is a measure
of the information available about the distribution of each parameter.

###  IV.3 Model criticism in gravitational-wave astronomy

It would be interesting to make a version of the left-hand panel of Fig. 4
using the population models from gravitational-wave astronomy discussed in the
previous subsection. Unfortunately, this would be quite challenging
computationally if we performed the calculation in the most straightforward
way possible. The first (easy) step would be to generate $N\approx 69$ events,
drawn at random from a hyper-posterior sample for one of the population models
and taking into account selection effects. The second (hard) step would be to
run Bayesian parameter on each event to generate posterior samples. The third
(medium-hard) step would be run a population study, which reads in the
posterior samples in order to calculate
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M)$.
In order to reject a model as misspecified with 99% confidence, we would have
to repeat this procedure 99 times.

While that is probably not realistic in the near future, it may be possible to
represent the likelihood functions of simulated events using a Fisher matrix
approximation, which would speed up the calculation significantly. However,
verifying that this approximation produces adequately estimates for
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M)$
could remain a challenge. Another possibility worthy of investigation is the
idea that the distribution of
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M)$
might have some quasi-universal properties. If it can be shown that a large
class of problems produce a similarly-shaped distribution of
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M)$,
perhaps a relatively small number of simulations can be used to work out the
shape of
$p(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$.
We leave this for future work.

##  V Conclusion

The  – $\mathcal{L}$ formalism provides a useful lens through which to view
population studies in gravitational-wave astronomy. It provides an upper bound
on the Bayesian evidence for population models  – $\mathcal{L}$ . The
associated pseudo-prior distribution  – $\pi$ is a sum of delta functions. The
– $\pi$ distribution can be used to see which features in a reconstructed
distribution are model-dependent, and which are genuinely present in the data.
The  – $\pi$ distribution can also draw attention to features in the data that
are not fit by current models, providing a tool for the design of new models.
Finally, the  – $\mathcal{L}$ formalism can be used to determine if a model is
misspecified, by comparing the values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},{\cal Z}_{M})$
to the expected distribution of these quanties given the model $M$. This
comparison can be made quantitatively with a $p$-value. And, by comparing the
measured values of
$(\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M))$
to the distribution expected given the model, it is possible to see the way in
which the model is misspecified. Constructing a distribution of
$\text{\ooalign{\raisebox{0.86108pt}{--}\cr$\mathcal{L}$\cr}},\mathcal{L}_{\textrm{max}}(M)$
may be computationally prohibitive in gravitational-wave astronomy, though,
future work is required to investigate simplifying assumptions that might
bring down the cost.

While we have introduced the  – $\mathcal{L}$ formalism within the context of
gravitational-wave astronomy, the framework is general, and we expect it can
be applied to a broad range of problems in astronomy and beyond where one
seeks to infer the distribution of parameters $\theta$ with potentially
unreliable hierarchical models.

## Acknowledgements

We thank Katerina Chatziioannou, Colm Talbot, Isaac Legred, Isobel Romero-
Shaw, and Paul Lasky for insightful discussions about the  – $\mathcal{L}$
formalism. We thank Rory Smith for input on early discussions regarding using
– $\mathcal{L}$ for model mispecification tests. We are grateful to Jacob
Golomb for discussions focused on computing a continuous representation of the
detection probability for gravitational-wave astronomy. We thank Tom Callister
for comments on an early version of the manuscript.

This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation. This research
has made use of data, software and/or web tools obtained from the
Gravitational Wave Open Science Center (https://www.gw-openscience.org), a
service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo
Collaboration. Virgo is funded by the French Centre National de Recherche
Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare
(INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian
institutes. The authors are grateful for computational resources provided by
the LIGO Laboratory and supported by National Science Foundation Grants
PHY-0757058 and PHY-0823459. This paper carries LIGO Document Number
#P2200309. E.T. is supported through Australian Research Council (ARC) Centre
of Excellence CE170100004.

## References

  * Thrane and Talbot (2019) E. Thrane and C. Talbot, “An introduction to Bayesian inference in gravitational-wave astronomy: parameter estimation, model selection, and hierarchical models,” Pub. Astron. Soc. Aust. 36, E010 (2019).
  * Vitale  _et al._ (2022) S. Vitale, D. Gerosa, W.M. Farr, and S.R. Taylor,  _Inferring the Properties of a Population of Compact Binaries in Presence of Selection Effects_ , Handbook of Gravitational Wave Astronomy (Springer, Singapore, 2022).
  * Mandel  _et al._ (2019) Ilya Mandel, Will M. Farr, and Jonathan R Gair, “Extracting distribution parameters from multiple uncertain observations with selection biases,” Mon. Not. R. Ast. Soc. 486, 1086 (2019).
  * Abbott  _et al._ (2021a) R. Abbott _et al._ , “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run,” (2021a), arxiv/2111.03606.
  * Note (1) We adopt the threshold utilized in (Abbott  _et al._ , 2021c) of a false-alarm-rate $<1\tmspace+{.1667em}\textrm{yr}^{-1}$.
  * Aasi  _et al._ (2015) J. Aasi _et al._ , “Advanced LIGO,” Class. Quantum Grav. 32, 074001 (2015).
  * Acernese  _et al._ (2015) F. Acernese _et al._ , “Advanced virgo: a second-generation interferometric gravitational wave detector,” Class. Quantum Grav. 32, 024001 (2015).
  * Akutsu  _et al._ (2021) T. Akutsu _et al._ , “Overview of KAGRA: Detector design and construction history,” Prog. Theor. Exp. Phys. 2021, 05A101 (2021).
  * Olsen  _et al._ (2022) S. Olsen, T. Venumadhav, J. Mushkin, J. Roulet, B. Zackay, and M. Zaldarriaga, “New binary black hole mergers in the ligo–virgo o3a data,” Phys. Rev. D 106, 043009 (2022).
  * Nitz  _et al._ (2021) A. H. Nitz, C. D. Capano, S. Kumar, Y. Wang, S. Kastha, M. Schäfer, R. Dhurkunde, and M. Cabero, “3-OGC: Catalog of gravitational waves from compact-binary mergers,” Astrophys. J. 922, 76 (2021).
  * Zackay  _et al._ (2021) B. Zackay, L. Dai, T. Venumadhav, J. Roulet, and M. Zaldarriaga, “Detecting gravitational waves with disparate detector responses: Two new binary black hole mergers,” Phys. Rev. D 104 (2021).
  * Venumadhav  _et al._ (2020) T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, “New binary black hole mergers in the second observing run of advanced LIGO and advanced virgo,” Phys. Rev. D 101, 083030 (2020).
  * Zackay  _et al._ (2019) B. Zackay, T. Venumadhav, L. Dai, J. Roulet, and M. Zaldarriaga, “Highly spinning and aligned binary black hole merger in the advanced LIGO first observing run,” Phys. Rev. D 100, 023007 (2019).
  * Abbott  _et al._ (2019) R. Abbott _et al._ , “Binary black hole population properties inferred from the first and second observing runs of advanced LIGO and advanced virgo,” Astrophys. J. 882, L24 (2019).
  * Abbott  _et al._ (2021b) R. Abbott _et al._ , “Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog,” Astrophys. J. 913, L7 (2021b).
  * Abbott  _et al._ (2021c) R. Abbott _et al._ , “The population of merging compact binaries inferred using gravitational waves through GWTC-3,” (2021c), arxiv/2111.03634.
  * Roulet  _et al._ (2021) J. Roulet, H. S. Chia, S. Olsen, L. Dai, T. Venumadhav, B. Zackay, and M. Zaldarriaga, “Distribution of effective spins and masses of binary black holes from the LIGO and virgo o1–o3a observing runs,” Phys. Rev. D 104, 083010 (2021).
  * Farr  _et al._ (2017) W. M. Farr, S. Stevenson, M. Coleman Miller, I. Mandel, B. Farr, and A. Vecchio, “Distinguishing spin-aligned and isotropic black hole populations with gravitational waves,” Nature 548, 426 (2017).
  * Talbot and Thrane (2018) C. Talbot and E. Thrane, “Measuring the binary black hole mass spectrum with an astrophysically motivated parameterization,” Astrophys. J. 856, 173 (2018).
  * Talbot and Thrane (2017) C. Talbot and E. Thrane, “Determining the population properties of spinning black holes,” Phys. Rev. D 96, 023012 (2017).
  * Callister  _et al._ (2021) T. A. Callister, C. Haster, K. K. Y. Ng, S. Vitale, and W. M. Farr, “Who ordered that? unequal-mass binary black hole mergers have larger effective spins,” Astrophys. J. Lett. 922, L5 (2021).
  * Fishbach  _et al._ (2022) M. Fishbach, C. Kimball, and V. Kalogera, “Limits on hierarchical black hole mergers from the most negative $\chi_{\text{eff}}$ systems,” Astrophys. J. Lett. 935, L26 (2022).
  * Biscoveanu  _et al._ (2021) S. Biscoveanu, M. Isi, S. Vitale, and V. Varma, “New spin on LIGO-virgo binary black holes,” Phys. Rev. D 126, 171103 (2021).
  * Biscoveanu  _et al._ (2022) S. Biscoveanu, T. A. Callister, C. Haster, K. K. Y. Ng, S. Vitale, and W. M. Farr, “The binary black hole spin distribution likely broadens with redshift,” Astrophys. J. Lett. 932, L19 (2022).
  * Vitale  _et al._ (2017) S. Vitale, R. Lynch, R. Sturani, and P. Graff, “Use of gravitational waves to probe the formation channels of compact binaries,” Class. Quantum Grav. 34, 03LT01 (2017).
  * Stevenson  _et al._ (2017) S. Stevenson, C. P. L. Berry, and I. Mandel, “Hierarchical analysis of gravitational-wave measurements of binary black hole spin–orbit misalignments,” Mon. Not. R. Ast. Soc. 471, 2801 (2017).
  * Miller  _et al._ (2020) S. Miller, T. A. Callister, and W. M. Farr, “The low effective spin of binary black holes and implications for individual gravitational-wave events,” Astrophys. J. 895, 128 (2020).
  * Galaudage  _et al._ (2021) S. Galaudage, C. Talbot, T. Nagar, D. Jain, E. Thrane, and I. Mandel, “Building better spin models for merging binary black holes: Evidence for nonspinning and rapidly spinning nearly aligned subpopulations,” Astrophys. J. Lett. 921, L15 (2021).
  * Fishbach  _et al._ (2018) M. Fishbach, D. E. Holz, and W. M. Farr, “Does the black hole merger rate evolve with redshift?” Astrophys. J. 863, L41 (2018).
  * Edelman  _et al._ (2022a) B. Edelman, Z. Doctor, J. Godfrey, and B. Farr, “Ain’t no mountain high enough: Semiparametric modeling of LIGO–virgo’s binary black hole mass distribution,” The Astrophysical Journal 924, 101 (2022a).
  * Edelman  _et al._ (2022b) B. Edelman, B. Farr, and Z. Doctor, “Cover Your Basis: Comprehensive Data-Driven Characterization of the Binary Black Hole Population,” (2022b), in prep.
  * Golomb and Talbot (2022a) J. Golomb and C. Talbot, “Searching for structure in the binary black hole spin distribution,” (2022a), in prep.
  * Note (2) Here, we paraphrase the aphorism attributed to statistician, George Box: “all models are wrong, but some are useful.”.
  * Romero-Shaw  _et al._ (2022) I. M. Romero-Shaw, E. Thrane, and P. D. Lasky, “When models fail: an introduction to posterior predictive checks and model misspecification in gravitational-wave astronomy,” Pub. Astron. Soc. Aust. 39, E025 (2022).
  * Gelman  _et al._ (2013) A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin,  _Bayesian Data Analysis, Third Edition_ , Chapman & Hall/CRC Texts in Statistical Science (Taylor & Francis, 2013).
  * Essick  _et al._ (2022) R. Essick, A. Farah, S. Galaudage, C. Talbot, M. Fishbach, E. Thrane, and D. E. Holz, “Probing Extremal Gravitational-Wave Events with Coarse-Grained Likelihoods,” Astrophys. J. 926, 34 (2022).
  * Callister  _et al._ (2022) T. A. Callister, S. J. Miller, K. Chatziioannou, and W. M. Farr, “No evidence that the majority of black holes in binaries have zero spin,” (2022), arxiv/2205.08574.
  * Tong  _et al._ (2022) H. Tong, S. Galaudage, and E. Thrane, “The population properties of spinning black holes using Gravitational-wave Transient Catalog 3,” (2022), arxiv/2209.02206.
  * Jeffreys (1961) H. Jeffreys,  _Theory of Probability_ , 3rd ed. (Oxford, 1961).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods 17, 261–272 (2020).
  * Conn  _et al._ (2000) A. R. Conn, N. I. M. Gould, and P. L. Toint,  _Trust Region Methods_ (Society for Industrial and Applied Mathematics, 2000) https://epubs.siam.org/doi/pdf/10.1137/1.9780898719857 .
  * Note (3) A method is “better” if it yields a larger value of  – $\mathcal{L}$ than another approach.
  * Talbot and Thrane (2022) C. Talbot and E. Thrane, “Flexible and Accurate Evaluation of Gravitational-wave Malmquist Bias with Machine Learning,” Astrophys. J. 927, 76 (2022), arXiv:2012.01317 [gr-qc] .
  * Golomb and Talbot (2022b) J. Golomb and C. Talbot, “Hierarchical Inference of Binary Neutron Star Mass Distribution and Equation of State with Gravitational Waves,” Astrophys. J. 926, 79 (2022b), arXiv:2106.15745 [astro-ph.HE] .
  * LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration (2021a) LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration, “Gwtc-3: Compact binary coalescences observed by ligo and virgo during the second part of the third observing run — parameter estimation data release,”  (2021a).
  * LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration (2021b) LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration, “The population of merging compact binaries inferred using gravitational waves through gwtc-3 \- data release,”  (2021b).
  * Pan  _et al._ (2014) Y. Pan, A. Buonanno, A. Taracchini, L. E. Kidder, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, “Inspiral-merger-ringdown waveforms of spinning, precessing black-hole binaries in the effective-one-body formalism,” Phys. Rev. D 89, 084006 (2014), arXiv:1307.6232 [gr-qc] .
  * Taracchini  _et al._ (2014) A. Taracchini _et al._ , “Effective-one-body model for black-hole binaries with generic mass ratios and spins,” Phys. Rev. D 89, 061502 (2014), arXiv:1311.2544 [gr-qc] .
  * Bohé  _et al._ (2017) A. Bohé _et al._ , “Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors,” Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Ossokine  _et al._ (2020) S. Ossokine _et al._ , “Multipolar Effective-One-Body Waveforms for Precessing Binary Black Holes: Construction and Validation,” Phys. Rev. D 102, 044055 (2020), arXiv:2004.09442 [gr-qc] .
  * Hannam  _et al._ (2014) M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, “Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms,” Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Pratten  _et al._ (2021) G. Pratten _et al._ , “Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes,” Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Talbot  _et al._ (2019) C. Talbot, R. J. E. Smith, E. Thrane, and G. B. Poole, “Parallelized Inference for Gravitational-Wave Astronomy ,” Phys. Rev. D 100, 043030 (2019).
  * Ashton  _et al._ (2019) G. Ashton _et al._ , “Bilby: A user-friendly Bayesian inference library for gravitational-wave astronomy,” Astrophys. J. Supp. 241, 27 (2019).
  * Romero-Shaw  _et al._ (2020) I. M. Romero-Shaw _et al._ , “Bayesian inference for compact binary coalescences with BILBY: Validation and application to the first LIGO–Virgo gravitational-wave transient catalogue,” Mon. Not. R. Ast. Soc. 499, 3295 (2020).
  * Speagle (2020) J. S. Speagle, “dynesty: a dynamic nested sampling package for estimating bayesian posteriors and evidences,” Monthly Notices of the Royal Astronomical Society 493, 3132–3158 (2020).
  * LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration (2021c) LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration, “Gwtc-3: Compact binary coalescences observed by ligo and virgo during the second part of the third observing run — o1+o2+o3 search sensitivity estimates,”  (2021c).
  * Damour (2001) T. Damour, Phys. Rev. D 64, 124013 (2001).
  * Note (4) In Eq. (45), $q\equiv m_{2}/m_{1}$ is mass ratio, $\chi_{1,2}$ are the dimensionless black hole spins, and $\theta_{1,2}$ are the spin vector tilt angles relative to the orbital angular momentum.
  * Wysocki  _et al._ (2019) D. Wysocki, J. Lange, and R. O’Shaughnessy, “Reconstructing phenomenological distributions of compact binaries via gravitational wave observations,” Phys. Rev. D 100, 043012 (2019).
  * Payne  _et al._ (2022) E. Payne, S. Hourihane, J. Golomb, R. Udall, D. Davis, and K. Chatziioannou, “The curious case of GW200129: interplay between spin-precession inference and data-quality issues,” (2022), arXiv:2206.11932 [gr-qc] .
  * Mould  _et al._ (2022) M. Mould, D. Gerosa, F. S. Broekgaarden, and N. Steinle, “Which black hole formed first? mass-ratio reversal in massive binary stars from gravitational-wave data,” (2022), arxiv/2205.12329.