# Binary vision: The merging black hole binary mass distribution via iterative
density estimation

Jam Sadiq Instituto Galego de Física de Altas Enerxías, Universidade de
Santiago de Compostela, Santiago de Compostela, Galicia, Spain  SISSA, Via
Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste  Thomas Dent
Instituto Galego de Física de Altas Enerxías, Universidade de Santiago de
Compostela, Santiago de Compostela, Galicia, Spain  Mark Gieles ICREA, Pg.
Luís Companys 23, E08010 Barcelona, Spain  Institut de Ciències del Cosmos
(ICCUB), Universitat de Barcelona (IEEC-UB), Martí Franquès 1, E08028
Barcelona, Spain

###### Abstract

Binary black hole (BBH) systems detected via gravitational-wave (GW) emission
are a recently opened astrophysical frontier with many unknowns and
uncertainties. Accurate reconstruction of the binary distribution with as few
assumptions as possible is desirable for inference on formation channels and
environments. Most population analyses have, though, assumed a power law in
binary mass ratio $q$, and/or assumed a universal $q$ distribution regardless
of primary mass. Kernel density estimation (KDE)-based methods allow us to
dispense with such assumptions and directly estimate the joint binary mass
distribution. We deploy a self-consistent iterative method to estimate this
full BBH mass distribution, finding local maxima in primary mass consistent
with previous investigations and a secondary mass distribution with a partly
independent structure, inconsistent with both power laws and with a constant
function of $q$. We find a weaker preference for near-equal mass binaries than
in most previous investigations; instead, the secondary mass has its own
“spectral lines” at slightly lower values than the primary, and we observe an
anti-correlation between primary and secondary masses around the
$\sim\\!10\,M_{\odot}$ peak.

Compact binaries, stellar-mass black holes, gravitational waves, density
estimation

O0pt O0pt$\mathrlap{\mathchoice{\raisebox{0.0pt}{$\displaystyle{\mathchar
22\relax\mkern-9.0mu}$}}{\raisebox{0.0pt}{$\textstyle{\mathchar
22\relax\mkern-9.0mu}$}}{\raisebox{0.0pt}{$\scriptstyle{\mathchar
22\relax\mkern-9.0mu}$}}{\raisebox{0.0pt}{$\scriptscriptstyle{\mathchar
22\relax\mkern-9.0mu}$}}}b$

##  1 Introduction

Ever since the first GW detection revealed a binary black hole source with
the—previously unsuspected—component masses of around 35 $M_{\odot}$ (Abbott
et al., 2016a, b), LIGO-Virgo-KAGRA observations of compact binaries have
continued to yield surprises, of which the binary mass distribution arguably
contains the most information bearing on formation environments and channels.
In the first three observing runs of Advanced LIGO (Aasi et al., 2015) and
Advanced Virgo (Acernese et al., 2015) the better part of 100 detections of
binary compact object mergers via gravitational wave (GW) emission were made,
as catalogued in the GWTC releases (Abbott et al., 2019, 2021, 2021a, 2021b).
Having a set of detected events it is possible to study population properties
of these compact binaries and eventually draw implications from these
properties on binary astrophysical formation and evolution. Detailed
investigations of the population properties of BBH mergers, the most commonly
detected source type, were undertaken in Abbott et al. (2020a, 2023a),
focusing on several population characteristics including their component
masses and spins and possible dependence on redshift.

Among parameters estimated and studied in connection with the population
properties of binary compact objects, the component masses are obtained with
least uncertainty. Many parameterized and semi- or non-parametric models have
been proposed to study the mass-dependence of the compact binary merger rate
or the mass distribution of the merger population (see Abbott et al., 2023a,
and references therein). In parametric models, Bayesian hierarchical
techniques are used to infer model hyper-parameter posteriors, and thus the
population distribution (e.g. Mandel, 2010; Thrane & Talbot, 2019). On the
other hand, non-parametric models are data driven methods which learn
population properties either without requiring any specific functional form,
or (for semi-parametric models) allowing for generalised deviations from a
given parametric model (Powell et al., 2019; Tiwari & Fairhurst, 2021; Tiwari,
2021, 2022; Veske et al., 2021; Rinaldi & Del Pozzo, 2021; Edelman et al.,
2021, 2023; Callister & Farr, 2023; Ray et al., 2023; Toubiana et al., 2023).

In Sadiq et al. (2022) we introduced a fast and flexible adaptive width kernel
density estimation (awKDE) as a non-parametric estimation method for
population reconstructions of binary black hole distribution from observed
gravitational wave data. A limitation of this method arose from the
measurement uncertainty in each individual event’s parameters. Given the
relatively low signal-to-noise ratios, the observed component masses have
significant uncertainties (e.g. Veitch et al., 2015) that can bias the overall
population distribution if not properly accounted for.

In this work we are proposing a new method to reduce this uncertainty estimate
in our populations distribution by using iterative re-weighting scheme of
samples for each observed events using the awKDE density estimates as a
probability for re-weighting of samples. The idea is similar to the standard
expectation-maximization algorithm (Dempster et al., 1977).

As an application of this new method, we estimate the full 2-dimensional
component mass distribution without assumptions (aside from the use of
Gaussian kernels) on its functional form. While attention has often focused on
the primary component mass or on the more precisely measured chirp mass (see
among others Dominik et al., 2015; Tiwari & Fairhurst, 2021; Tiwari, 2022,
2023; Edelman et al., 2023; Schneider et al., 2023; Farah et al., 2023), less
attention has been paid to the full binary distribution, either considered via
the secondary $m_{2}$ or mass ratio $q\equiv m_{2}/m_{1}$. We expect these
parameters to bear traces of the possible BBH formation channels (Kovetz et
al., 2017), in that for dynamical (cluster) formation the two masses may be
independent variates, up to a factor modelling probability of binary formation
and merger (Fishbach & Holz, 2020; Antonini et al., 2023) that typically
favors near-equal masses (e.g. Rodriguez et al., 2016; O’Leary et al., 2016);
whereas for isolated binary evolution, some nontrivial though probably highly
model-dependent correlation of component masses may arise (e.g. van Son et
al., 2022).

Typically, parametric models have assumed power-law $m_{2}$ dependence at
fixed $m_{1}$ (Kovetz et al., 2017; Fishbach & Holz, 2017; Talbot & Thrane,
2018; Abbott et al., 2019), recovering mildly positive powers indicating a
preference for equal masses. A more detailed study using GWTC-1 events
concluded that the two BHs of a given binary prefer to be of comparable mass
(Fishbach & Holz, 2020). More recent non-parametric or semi-parametric studies
have relaxed these assumptions, either through allowing the power-law index to
vary over chirp mass (Tiwari, 2022), or allowing $p(q)$ to be a free (data-
driven) function (Edelman et al., 2023; Callister & Farr, 2023) though
enforcing the same dependence over all primary masses. Tiwari (2023)
introduced a more flexible approach with $p(q)$ modelled by a truncated
Gaussian whose parameters depend on chirp mass, finding some significant
variation. Ray et al. (2023) measured the full 2-d distribution with a binned
(piecewise-constant) model over $m_{1}$, $m_{2}$ (including possible redshift
dependence), although they did not consider the $q$ distribution. Note that
the mass ratio distribution presents nontrivial technical issues since (at
least for lower mass BBH) typical event measurement errors are both large, and
correlated with the BH orbit-aligned spin components (Cutler & Flanagan, 1994;
Baird et al., 2013). Concerning measurement errors, we expect our iterative
reweighting scheme to yield a significant advance in reconstructing the full
mass distribution.

The remainder of the paper is organized as follows: in section 2 we explain
our method and demonstrate it using simple one-dimensional mock data. In
section 3 we apply our method to detected BBH in GWTC-3, compare the result
with our previous studies and further use our new method in two mass
dimensions. In section 4 we discuss the implications of our results and
consider extensions of the method. We also describe additional mock data tests
and supplementary results in appendices.

##  2 Method

###  2.1 Statistical framework

Our general approach to population inference can be considered as similar to
maximum likelihood, with uncertainties quantified via empirical bootstrap
methods (Efron, 1979). Given a set of observed events, if we neglect
measurement uncertainty in each event’s parameters, our population estimate is
a KDE where the kernel bandwidth for each event is adjusted (Breiman et al.,
1977; Abramson, 1982; Terrell & Scott, 1992; Sain & Scott, 1996) using an
adaptive scheme to maximize the cross-validated likelihood (Sadiq et al.,
2022). Note that a “maximum likelihood” KDE is not well defined, as the
likelihood increases indefinitely in the limit of small bandwidth kernels
centered on the observations; in this limit, the variance of the density
estimate over realizations of the data becomes infinitely large. The bias-
variance tradeoff is then addressed by adaptive kernel bandwidth, with a
choice of hyperparameters—global bandwidth $h$ and sensitivity parameter
$\alpha$, see Wang & Wang (2011)—optimized by grid search using leave-out-one
cross-validation to calculate a figure of merit. We then quantify counting
uncertainties for the underlying inhomogeneous Poisson process using
generalized bootstrap resampling Chamandy et al. (2012).

We noted in Sadiq et al. (2022) that for nontrivial measurement errors this
method, in addition to possible intrinsic biases due to the choice of a
Gaussian KDE, will be biased towards an over-dispersed estimate of the true
distribution. Here we motivate and present our strategy for correcting this
bias. Our motivation is linked to Bayesian hierarchical population inference
(Mandel, 2010), where measurement errors are treated by considering the true
event properties $\vec{\theta_{i}}$ for events labelled $i=1,\ldots,N$ as
nuisance parameters in the likelihood:

$p(\\{d\\}|\vec{\lambda})=\prod_{i=1}^{N}\int
p(d_{i}|\vec{\theta}_{i})p(\vec{\theta}_{i}|\vec{\lambda})\,d\vec{\theta}_{i},$
(1)

where $\\{d\\}$ is the set of data segments $d_{i}$ corresponding to the
events and $\vec{\lambda}$ are population model hyperparameters (here for
simplicity we omit selection effects). Inference is implemented using
parameter estimation (PE) samples which were generated using a standard or
fiducial prior $p_{\rm PE}(\vec{\theta})$, often chosen as uniform over
parameters of interest (see e.g. Veitch et al., 2015; Thrane & Talbot, 2019).
Samples (labelled by $k$) are distributed as the posterior density using this
prior, hence

$\vec{\theta}_{i}^{k}\sim p(\vec{\theta}_{i}|d_{i},p_{\rm
PE}(\vec{\theta}))\propto
p(d_{i}|\vec{\theta}_{i})p(\vec{\theta}_{i}|\vec{\lambda})\cdot\frac{p_{\rm
PE}(\vec{\theta})}{p(\vec{\theta}_{i}|\vec{\lambda})},$ (2)

hence the integrals may be performed (up to a constant factor) by summing over
samples _re-weighted_ by the ratio of the population distribution to the PE
prior, $p(\vec{\theta}_{i}|\vec{\lambda})/p_{\rm PE}(\vec{\theta})$.

Here, while not making use of this hierarchical likelihood, we remark that PE
samples give a _biased_ estimate of each event’s properties if the true
population distribution $p_{\rm pop}(\vec{\theta})$ (corresponding to
$p(\vec{\theta}|\vec{\lambda})$ in the parameterised case) is not equal to
$p_{\rm PE}(\vec{\theta})$. Then, if we have access to an estimated population
distribution $\hat{p}_{\rm pop}(\vec{\theta})$ that is more accurate than the
PE prior is, we will obtain more accurate estimates of event properties by
drawing samples weighted proportional to $\hat{p}_{\rm
pop}(\vec{\theta})/p_{\rm PE}(\vec{\theta})$, as described in more detail
below,

To summarize, a KDE obtained by drawing from PE samples will be biased because
these samples are themselves biased, due to the PE prior not being equal to
the true population distribution. However, the more accurate an estimate of
the true distribution we are able to obtain, the smaller will be the bias in
event parameters using reweighted PE samples, and ultimately the smaller will
be the bias of the KDE.

###  2.2 Iterative Reweighting

The above discussion suggests an iterative procedure where, beginning with
both biased PE samples and a biased population KDE, one may be improved in
turn using the other, until – ideally – reaching a stationary state, where
both the sample draws and the corresponding population estimates are unbiased
(up to more fundamental limitations of PE and of our KDE). This iterative
strategy is similar to the Expectation-Maximization (EM) algorithm (Dempster
et al., 1977), a popular method to estimate parameters for statistical models
when there are missing or incomplete data.

Our basic algorithm follows these steps:

  1. 1.

For each GW event, draw Poisson distributed (with mean 1) PE samples weighted
by the current estimate of population density $\hat{p}_{\rm pop}$

  2. 2.

Create an awKDE from this sample set, optimizing the global bandwidth (and
sensitivity parameter $\alpha$, if not fixed)

  3. 3.

Update the current population estimate using one or more KDEs and the
selection function, and go to step 1.

In more detail, in step 1 we draw PE samples with probability proportional to
the ratio of $\hat{p}_{\rm pop}(\vec{\theta}_{i}^{k})$ to the PE prior
distribution. Step 2 reproduces our previous awKDE method. Step 3 relates the
KDE of _detected_ events to an estimate of the true population distribution,
hence in general it requires us to compensate for the selection function over
the event parameter space: i.e. we estimate the true distribution by the KDE
of detected events divided by the probability of detection, as detailed in
(Sadiq et al., 2022, section 3.1).

In step 3 we may choose to derive the updated population density $\hat{p}_{\rm
pop}(\vec{\theta})$ from only the most recently calculated KDE: then the
iterative process is a Markov chain,111Although it may be thought of as a
Markov chain Monte Carlo, our method is entirely unrelated to the Metropolis-
Hastings algorithm. and we may characterize it via the autocorrelation of
various scalar quantities computed at each iteration. We use the optimized
global bandwidth $h$ (and adaptive sensitivity parameter $\alpha$, if not
fixed to unity) for this purpose.

After discarding a small number of initial iterations and then accumulating a
number significantly greater than the autocorrelation time, we expect the
collection of iterations to provide unbiased (though not necessarily
independent) estimates of the population distribution. For subsequent
iterations we then use the median of $\hat{p}_{\rm pop}(\vec{\theta}_{i}^{k})$
over a buffer of previous iterations (usually the previous 100) to determine
the sample draw probabilities for the next iteration. This population estimate
should be more precise than one using only a single previous KDE; and in
addition using the buffer estimate the samples for each successive iteration
are essentially independent.

###  2.3 One-dimensional mock data demonstration

Figure 1: awKDE for 60 events from a mock data mixture distribution with 50%
power law ($\alpha=-0.5$) and 50% Gaussian ($\mu=35$, $\sigma=3$) samples. 100
mock parameter samples are generated for each event with measurement error
$\sigma_{m}=5$. Top: awKDE drawing random parameter samples without
reweighting. Bottom: applying iterative reweighting. The solid (dot-dashed)
lines represent the median (symmetric 90% confidence band) from 900 bootstrap
iterations.

We first test this iterative reweighting method on a simple mock dataset. We
generate true event parameters drawing 30 event each from a truncated power
law and a Gaussian distribution respectively; the power law is $p(x)\sim
x^{-0.5}$ and the Gaussian has mean (s.d.) of $\mu=35$ ($\sigma_{p}=3$). We
then add measurement errors to our true parameters with a s.d. $\sigma_{m}=5$,
hence broader than the true Gaussian peak; 100 mock parameter samples with the
same uncertainty are then generated around each “measured” value. First,
applying awKDE as in Sadiq et al. (2022) to random draws from these mock
parameter samples, as expected we find an over-dispersed estimate around the
peak (Fig. 1, top).

The second (bottom) plot of Fig. 1 shows the awKDE applying our iterative
reweighting algorithm. Here the Gaussian height and s.d. are accurately
reconstructed and the true distribution is well within the 90% percentiles of
iteration samples, except at the step-function truncation of the power law
which cannot be accurately represented by a Gaussian KDE.

Figure 2: Autocorrelation of the (log) optimized global bandwidth series for
the iterative Markov chain. The x-axis shows lag, i.e. separation of the
iterations between which autocorrelation is calculated.

We verify that the initial Markov process has accumulated several independent
samples by plotting the autocorrelation of the optimized global bandwidth $h$
vs. lag (separation along the chain) in Fig. 2. (We fix the adaptive
sensitivity parameter $\alpha$ to unity for 1-d data.) The autocorrelation
drops near zero by a lag of $\lesssim\\!10$ iterations, thus a buffer of 100
iterations contains several independent population estimates. The estimate is
more noisy for larger lags as fewer iterations are available.

More detailed tests of iterative reweighting with a Gaussian KDE in two
dimensions in the presence of correlated parameter errors are given in
Appendix A.

##  3 Results from GWTC-3

As in Sadiq et al. (2022), as input to our analysis we use parameter
estimation samples (LVK, 2021a) for the set of BBH events catalogued in GWTC-3
(Abbott et al., 2021a, 2023b) with false alarm rate below 1 per year. For the
sensitivity estimate we employ a fit to search injection (simulated signal)
results (Wysocki et al., 2019; Wysocki, 2020) released with the catalog (LVK,
2021b).

###  3.1 One-dimensional mass distribution

We start by evaluating the effect of the iterative reweighting method on the
1-d primary mass distribution, taking 100 random PE samples for each of the 69
BBH events; here, we assume a power-law distribution for secondary mass
$p(m_{2})\sim q^{1.5}$. We reproduce the awKDE results from Sadiq et al.
(2022) and use this estimate to seed the reweighted iteration algorithm. After
$\sim$1000 reweighting iterations we compute the median and symmetric 90%
interval from the last 900 rate estimates (the first 100 are used to set up a
buffer for population weighting, as above): results are presented in Fig. 3.

Figure 3: Differential rate over BBH primary mass using the iterative weighted
KDE on GWTC-3 PE samples, assuming a power-law secondary distribution. We
overplot our estimates with two of the models, Flexible mixtures (FM), and
Power Law + Spline (PS) models in (Abbott et al., 2023a). Our estimate is
generally consistent with other non-parametric methods, though with a higher
and narrower peak around $35\,M_{\odot}$ and lacking any feature between the
10 $M_{\odot}$ and 35 $M_{\odot}$ peaks.

Our estimate is generally consistent with other non-parametric or semi-
parametric approaches, represented by the Flexible mixtures, and Power Law +
Spline models in Abbott et al. (2023a), and does not show the over-dispersion
apparent in Figure 8 of Sadiq et al. (2022); specifically, we find a slightly
higher and narrower peak around $35\,M_{\odot}$, but no identifiable feature
around $20\,M_{\odot}$ (compare Tiwari, 2022; Toubiana et al., 2023).

###  3.2 Two—dimensional Mass Reconstruction

Next, we apply our reweighting scheme on PE samples for both component masses
and compute the two-dimensional (2-d) merger rate, using the estimated
sensitive volume$\times$time (VT) as a function of the two masses. We will
first discuss various technical aspects of extending the 1-d calculation
without assuming any power-law dependence for $m_{2}$.

#### Binary exchange symmetry

Typically when presenting binary parameter estimates, the convention
$m_{1}>m_{2}$ is applied. However, all aspects of binary formation physics and
event detection and parameter estimation will be invariant under swapping the
component labels, i.e. exchanging $m_{1}\leftrightarrow m_{2}$ (and at the
same time exchanging the spins). Thus, considering the differential merger
rate $\mathcal{R}(m_{1},m_{2})$ as a function over the whole plane, it must
also have a reflection symmetry about the line $m_{1}=m_{2}$. To respect this
symmetry and remove biases resulting from the apparent lack of support at
$m_{2}>m_{1}$, we train and evaluate KDEs on _reflected_ sample sets which
contain both the released PE samples, and copies of them with swapped
components. Note also that a power-law $m_{2}$ distribution implies the rate
is a non-differentiable function at the equal mass line, whereas a KDE by
construction is smooth and differentiable everywhere.

#### Choice of KDE parameters

In previous work (Sadiq et al., 2022), we mainly considered a KDE constructed
over linear mass (or distance) parameters; however, here we choose the
logarithms of component masses. While this choice is not expected to have a
large impact on the results, since the kernel bandwidth is free to locally
adapt in either case, it is technically preferable for a few reasons: we avoid
any possible KDE support at negative masses; there is a generally higher
density of events towards lower masses considering the entire
$3-100\,M_{\odot}$ range; the density of observed events also shows less
overall variation over log coordinates; and when evaluating the KDE on a grid
with equal spacing, fewer points are required to maintain precision for the
low-mass region.

For a 2-d KDE we also have a choice of kernel parameters, i.e. the Gaussian
covariance matrix: given the similar or identical physical interpretation and
range of values between $\ln m_{1}$ and $\ln m_{2}$, we choose a covariance
proportional to the unit matrix, with an overall factor determined by the
local adaptive bandwidth for each event.

#### PE prior

The PE samples released by LVK use a prior uniform in component masses (LVK,
2021a) up to a factor dependent on cosmological redshift; we currently do not
consider reweighting relative to the default PE cosmological model. As the
prior is a density, it transforms with a Jacobian factor when changing
variables to $\ln m_{1},\ln m_{2}$, thus, we must divide the estimated rate
$\mathcal{R}(\ln m_{1},\ln m_{2})$ by a prior $\propto m_{1}m_{2}$ when
obtaining reweighted draw probabilities.

Figure 4: Two dimensional rate estimated from PE samples for 69 BBH events
from GWTC-3. Red $+$ symbols show the median masses for each event. The
contours and color scale shows the two-dimensional differential merger rate
over $\ln m_{1},\ln m_{2}$ from iterative reweighted KDE. Two main maxima are
visible at $m_{1}\,(m_{2})\sim 10\,(8)\,M_{\odot}$ and
$\sim\\!35\,(32)\,M_{\odot}$ with a possible less significant overdensity
around $m_{1}\sim 20\,M_{\odot}$.

With these technical choices, we perform 1500 reweighting iterations in total,
the first 600 using the Markov chain (i.e. the immediately preceding rate
estimate) for sample draw weights, and the remaining 900 using the buffer
median estimate.

Figure 5: Autocorrelations of the KDE (log) global bandwidth and adaptive
parameter $\alpha$ using the initial 600 Markov chain iterations. The
autocorrelation drops close to 0 by a lag of $\sim$30 iterations. The estimate
becomes noisy at high lags as a smaller number of iterations is available.

Fig. 4 shows the rate estimate computed with iterative reweighting for BBH
events in GWTC-3. The autocorrelations of optimal global bandwidth and
sensitivity parameter $\alpha$ for the first 600 iterations are shown in Fig.
5: the correlation drops close to $0$ at a lag of $\sim\,$30 iterations.

The mass distribution shows several interesting features in addition to the
expected peaks (overdensities) around primary masses of $\sim\\!10\,M_{\odot}$
and $\sim\\!35\,M_{\odot}$, with corresponding peaks over secondary mass. For
primary masses $\sim\\!35\,M_{\odot}$ up to $80\,M_{\odot}$, the most likely
secondary mass is $\sim\\!30\,$–$\,35\,M_{\odot}$. Thus, over this range the
two component masses appear almost independently chosen. Around the $m_{1}\sim
10\,M_{\odot}$ peak, there appears to be some _anti-correlation_ of the two
components, i.e. higher $m_{1}$ favors lower $m_{2}$. Between the two peaks
the distribution of mass ratios appears broader than at either one (as hinted
at in Tiwari, 2023), although the apparent trend is based on a small number of
events. We also see a narrow lower density region just above the
$\sim\\!10\,M_{\odot}$ peak (cf. the local minimum at chirp mass
$\sim\\!11\,M_{\odot}$ in Tiwari 2023).

Figure 6: Merger rates over component masses, obtained by numerical
integration of our 2-d KDE rate estimates in Fig. 4: the medians and symmetric
90% confidence regions are shown.

We also integrated the 2-d KDE rate estimate numerically over both $m_{1}$ and
$m_{2}$ to obtain merger rates over component mass. As shown in Fig. 6, we
recover features consistent with the 2-d estimates and with other methods.
Each component mass distribution appears well modelled by a combination of two
Gaussian peaks and (broken) power laws.

Figure 7: Secondary mass $m_{2}$ distributions estimated via iterative
reweighted KDE for various fixed values of primary mass: for each $m_{1}$ (in
$M_{\odot}$) we plot the median and symmetric $90\%$ confidence band. Notice
the top plot x-xais is in log scale unlike the bottom plot with linear scale
in $x$-axis.

To elucidate features in the 2-d distribution, we choose various
representative values of primary mass to plot the distribution of $m_{2}$ in
Fig. 7. The similarity between secondary distributions for $m_{1}\gtrsim
35\,M_{\odot}$ is evident.

Figure 8: Mass ratio $q$ distributions estimated via iterative reweighted KDE
for various fixed values of primary mass: for each $m_{1}$ (in $M_{\odot}$) we
plot the median and symmetric $90\%$ confidence band. For comparison we
overplot a power law dependence $\propto q^{1.5}$.

We may also derive the distribution of mass ratio $q$ from our 2-d rate
estimate. We plot this for various representative values of primary mass in
Fig. 8, and compare to a typical power law $\propto q^{1.5}$. First, we see
that the $q$ distribution varies over primary mass; hence, models where it is
forced to the same form over the whole mass range are likely to have
nontrivial bias. For some primary masses, $p(q)$ is consistent with a
monotonic increasing function such as a positive power, but for others it
clearly decreases over some of the range. Roughly, if $m_{1}$ is close to a
peak then $p(q)$ is consistent with an increasing power law, but for other
values the mass ratio rather shows a maximum at intermediate values, down to
$q\sim 0.5$ for $m_{1}=15$ or $m_{1}=70$.

This behaviour may suggest that the primary and secondary masses are
independently drawn from similar distributions, modulo a $q$-dependent
“pairing factor” (Fishbach & Holz, 2020) which influences the relative
probability of binary merger. However the preference towards $q\sim 1$ seen in
previous work is not confirmed here. Callister & Farr (2023) reached a similar
conclusion although assuming a “universal” $p(q)$ over all primary masses.

#### Results including GW190814

We also estimate the 2-d and 1-d integrated merger rates using our iterative
reweighted KDE method when the outlier event GW190814 (Abbott et al., 2020b),
which has a mass ratio $\mathcal{O}(10)$ and a secondary mass barely above the
likely neutron star maximum mass, is included in the BBH population. Detailed
results are presented in Appendix B: roughly summarizing the trends seen
there, the bulk of the estimated distribution remains little changed by the
addition of the extra event, although the peak in secondary mass below
$10\,M_{\odot}$ is shifted towards lower values and both higher and broader;
this is likely due to a general increase in KDE bandwidth when optimized with
cross-validation. (In parameterized models the estimated mass distribution is
also highly sensitive to inclusion of GW190814 (Abbott et al., 2020a, 2023a).)
It is not clear whether other methods for bandwidth choice would yield more
accurate estimates; higher event statistics in the low $m_{2}$ regime are
clearly desirable.

##  4 Discussion

#### Summary of results

In this work we undertook a detailed investigation of the full 2-d mass
distribution of merging compact binary black holes observed by LIGO-Virgo-
KAGRA up to the O3 run without assuming any specific functional form for the
secondary mass or mass ratio, enabled by a new method of iterative density
estimation to address mass measurement uncertainties. Although we reproduce
the broad features and local maxima seen in other parametric and non-
parametric analyses, we find significantly less preference for near-equal
masses than in most previous works; we also find that the mass ratio
distribution cannot be described by a single function over the whole
population (compare Tiwari, 2023). For a range of primary masses, we find non-
monotonically varying secondary and mass ratio distributions, thus a power-law
dependence is ruled out. Furthermore, we find that for primary masses above
$35\,M_{\odot}$ the secondary mass distribution is nearly independent of
$m_{1}$, with a “preferred partner” mass of $m_{2}\simeq 30-35\,M_{\odot}$.
Conversely, near the low-mass peak $m_{1}\simeq 10\,M_{\odot}$ we observe an
anticorrelation between the two components, i.e. higher $m_{1}$ implies lower
$m_{2}$.

#### Possible astrophysical interpretations

Our new estimate of the joint $m_{1}$–$m_{2}$ distribution may be compared to
model predictions in the literature; because our individual component marginal
distributions are similar to previous findings, we focus here on the mass
ratio. Broadly speaking, we can distinguish model predictions from the
isolated binary and dynamical channels.

The isolated binary channel predicts relatively flat $p(q)$ distributions
(e.g. Belczynski et al., 2020; Olejak et al., 2021) compared to the dynamical
channel (see Fig. 2 in Baibhav et al. 2023 and Fig. 1 in Zevin et al. 2021).
Our estimates show in general flatter $q$ distributions than the GWTC-3
results presented in Abbott et al. (2023a).

Because predictions for $p(q)$ in the isolated binary channel depend strongly
on the adopted parameters (see e.g. Broekgaarden et al., 2022), our results
provide an important step towards constraining astrophysical parameters with
GWs. For example, the steep $p(q)$ found for very small common envelope
efficiency parameter ($\alpha_{\rm CE}\simeq 0.2$, Baibhav et al., 2023) and
the chemically homogeneous evolution model (Mandel & de Mink, 2016) seem
disfavoured, implying that these routes cannot account for the majority of the
observed population.

The stable mass transfer channel is efficient for primary masses near
$\sim\\!10\,M_{\odot}$. van Son et al. (2022) predicts a dearth of near-equal
mass mergers, which is because the binary needs to be relatively unequal in
mass during the second mass transfer phase for the orbit to shrink, but not
too unequal to avoid unstable mass transfer. This is only partly supported by
our Fig. 4 in that the low-mass peak has support from equal mass out to
$q\simeq 0.5$. For some parameter choices their models predict bi-modality in
$p(q)$, with peaks at $q\simeq 0.35$ and $q\simeq 0.75$: our results suggest a
peak at $q\simeq 0.8$ for $m_{1}\simeq 10\,M_{\odot}$ and at $q\simeq 0.45$
for $m_{1}\simeq 15\,M_{\odot}$ (see Fig. 8), suggesting that a more detailed
comparison may yield interesting constraints.

For the dynamical channel, it is interesting to consider whether models now
predict $q$ distributions that are too steeply rising. Rodriguez et al. (2016)
modelled BBH mergers that formed dynamically in globular clusters: they find a
median mass ratio of 0.87, with 68% of sources having mass ratios $q>0.8$. As
shown in Fig. 8 we find comparable support for near-equal mass only at
$m_{1}\sim 10\,M_{\odot}$ or $\sim\\!35\,M_{\odot}$; elsewhere our median $q$
is significantly lower.

Antonini et al. (2023) model BBH merger in globular clusters in comparison to
the GWTC-3 $q$ distribution: their model distributions are flatter and
underestimate the power-law LVK fits by an order of magnitude at $q\simeq 1$.
They find a final $q$ distribution flatter than the $q$ distribution of
dynamically formed BBHs ($p(q)\propto q^{4}$ for metal-poor clusters) because
the BH mass function is not always sufficiently sampled, such that a secondary
BH with a mass similar to the primary is present in a cluster; and due to a
slight bias against equal-mass BBH due to their lower inspiral probability.
The reported $p(q)$ in their Fig. 1 is relatively flat for $q\gtrsim 0.7$,
qualitatively in agreement with our findings for $m_{1}\gtrsim 20\,M_{\odot}$
(Fig. 8, lower panel; their models cannot reproduce observed rates for lower-
mass primaries.) Due to the predicted pair instability gap, all BBH mergers in
their models with $m_{1}\gtrsim 50\,M_{\odot}$ are hierarchical mergers, i.e.
a BBH in which at least one of the components is a BBH merger remnants that
was retained in the cluster (e.g. Antonini & Rasio, 2016; Rodriguez et al.,
2019; Kimball et al., 2021). Mergers with second-generation primaries are
expected to have a mass ratio $q\simeq 0.5$, which is supported by our
distribution for $m_{1}=70\,M_{\odot}$ (Fig. 8).

The $p(q)$ distribution is expected to be slightly flatter for dynamically
formed BBHs in young (open) star clusters, because they have fewer BHs per
cluster and their higher metallicities lead to steeper BH mass functions and
therefore lower companion masses (e.g. Banerjee, 2021). A accurate picture of
the mass ratio distribution is therefore important to understand the relative
contribution of dynamically formed BBHs in young (and metal-rich) and old (and
metal-poor) star clusters; and also more generally for understanding the
relative contributions of isolated and dynamically formed binaries in the
population as a whole (Zevin et al., 2021; Baibhav et al., 2023).

An intriguing apparent feature in our reconstruction, the anticorrelation
between $m_{1}$ and $m_{2}$ in the low-mass ($m\sim 10\,M_{\odot}$) peak,
suggests a connection to isolated binary dynamics, though it would be
premature to link it with a specific mechanism.

#### Technical issues and biases

As noted in the introduction, measurement errors of the binary mass ratio are
correlated with those in (orbit-aligned) spins. Since we have so far not
attempted to reconstruct or estimate the merging binary spin distribution, we
implicitly assume that distribution is equal to the prior used for parameter
estimation (uniform in magnitude and isotropic in direction): this is a
potential source of bias which remains to be addressed by future work. The
distribution of aligned spins has been found to be concentrated near zero,
with a slight preference for positive aligned spin (Miller et al., 2020;
Abbott et al., 2020a, 2023a); hence, the degree of bias may be limited.
Callister et al. (2021) also note the intriguing possibility that the _true_
mass ratio and aligned spin (after allowing for measurement errors) are anti-
correlated.

A converse question concerns inferences on BH spin distributions which either
assume a specific distribution in $q$, or a power law with index as a
hyperparameter: if the $p(q)$ model is significantly inaccurate, are such spin
inferences biased? (Ng et al. 2018 and Miller et al. 2020 contain detailed
discussion of potential biases in aligned spin population estimates.) The
effect may not be large, as most BBH events by necessity have parameter values
close to the observed peaks, for which we find a $q$ distribution which is not
far from power-law.

#### Extensions of the method

As already noted, here we restricted the application of our KDE to the binary
mass distribution; component spins, and distance or redshift are then the next
relevant parameters for population analysis. We expect to encounter a
technical issue in optimizing the Gaussian kernel for a multi-dimensional data
set, where it will not be appropriate (or even meaningful, given the different
units) to impose equal variances over different parameters as we currently do
for (log) $m_{1}$ and $m_{2}$. For more than two dimensions a grid search may
not be practicable; more sophisticated methods may be required in order to
realize the potential of iterative KDE over a full set of population
parameters.

## Acknowledgements

We thank Daniel Wysocki for making available fitted sensitivity estimates for
binary mergers in the O1-O3 data. We also benefited from conversations with
Lieke van Son, Floor Broekgaarden and Fabio Antonini, and with Will Farr,
Thomas Callister, Amanda Farah, Vaibhav Tiwari and others in the LVK Binary
Rates & Populations group. This work has received financial support from Xunta
de Galicia (CIGUS Network of research centers), by European Union ERDF and by
the “María de Maeztu” Units of Excellence program CEX2020-001035-M and the
Spanish Research State Agency. TD and JS are supported by research grant
PID2020-118635GB-I00 from the Spanish Ministerio de Ciencia e Innovación. JS
also acknowledges support from the European Union’s H2020 ERC Consolidator
Grant “GRavity from Astrophysical to Microscopic Scales” (GRAMS-815673) and
the EU Horizon 2020 Research and Innovation Programme under the Marie
Sklodowska-Curie Grant Agreement No. 101007855. MG acknowledges support from
the Ministry of Science and Innovation (EUR2020-112157, PID2021-125485NB-C22,
CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033) and from AGAUR
(SGR-2021-01069).

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 research has made use of data or software obtained from the
Gravitational Wave Open Science Center (gwosc.org), a service of the LIGO
Scientific Collaboration, the Virgo Collaboration, and KAGRA. This material is
based upon work supported by NSF’s LIGO Laboratory which is a major facility
fully funded by the National Science Foundation, as well as the Science and
Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-
Society (MPS), and the State of Niedersachsen/Germany for support of the
construction of Advanced LIGO and construction and operation of the GEO600
detector. Additional support for Advanced LIGO was provided by the Australian
Research Council. Virgo is funded, through the European Gravitational
Observatory (EGO), by the French Centre National de Recherche Scientifique
(CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch
Nikhef, with contributions by institutions from Belgium, Germany, Greece,
Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported
by Ministry of Education, Culture, Sports, Science and Technology (MEXT),
Japan Society for the Promotion of Science (JSPS) in Japan; National Research
Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia
Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.

## References

  * Aasi et al. (2015) Aasi, J., et al. 2015, Class. Quant. Grav., 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2016a) Abbott, B. P., et al. 2016a, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2016b) —. 2016b, Astrophys. J. Lett., 818, L22, doi: 10.3847/2041-8205/818/2/L22
  * Abbott et al. (2019) —. 2019, Phys. Rev. X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2019) —. 2019, Astrophys. J. Lett., 882, L24, doi: 10.3847/2041-8213/ab3800
  * Abbott et al. (2020a) Abbott, R., et al. 2020a.  https://arxiv.org/abs/2010.14533
  * Abbott et al. (2020b) —. 2020b, Astrophys. J. Lett., 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abbott et al. (2021) —. 2021, Phys. Rev. X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2021a) —. 2021a.  https://arxiv.org/abs/2111.03606
  * Abbott et al. (2021b) —. 2021b.  https://arxiv.org/abs/2108.01045
  * Abbott et al. (2023a) —. 2023a, Phys. Rev. X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * Abbott et al. (2023b) —. 2023b.  https://arxiv.org/abs/2302.03676
  * Abramson (1982) Abramson, I. S. 1982, The Annals of Statistics, 10, 1217 , doi: 10.1214/aos/1176345986
  * Acernese et al. (2015) Acernese, F., et al. 2015, Class. Quant. Grav., 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Antonini et al. (2023) Antonini, F., Gieles, M., Dosopoulou, F., & Chattopadhyay, D. 2023, MNRAS, 522, 466, doi: 10.1093/mnras/stad972
  * Antonini & Rasio (2016) Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187, doi: 10.3847/0004-637X/831/2/187
  * Baibhav et al. (2023) Baibhav, V., Doctor, Z., & Kalogera, V. 2023, Astrophys. J., 946, 50, doi: 10.3847/1538-4357/acbf4c
  * Baird et al. (2013) Baird, E., Fairhurst, S., Hannam, M., & Murphy, P. 2013, Phys. Rev. D, 87, 024035, doi: 10.1103/PhysRevD.87.024035
  * Banerjee (2021) Banerjee, S. 2021, MNRAS, 503, 3371, doi: 10.1093/mnras/stab591
  * Belczynski et al. (2020) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104, doi: 10.1051/0004-6361/201936528
  * Breiman et al. (1977) Breiman, L., Meisel, W., & Purcell, E. 1977, Technometrics, 19, 135.  http://www.jstor.org/stable/1268623
  * Broekgaarden et al. (2022) Broekgaarden, F. S., Berger, E., Stevenson, S., et al. 2022, MNRAS, 516, 5737, doi: 10.1093/mnras/stac1677
  * Callister & Farr (2023) Callister, T. A., & Farr, W. M. 2023.  https://arxiv.org/abs/2302.07289
  * Callister et al. (2021) Callister, T. A., Haster, C.-J., Ng, K. K. Y., Vitale, S., & Farr, W. M. 2021, Astrophys. J. Lett., 922, L5, doi: 10.3847/2041-8213/ac2ccc
  * Chamandy et al. (2012) Chamandy, N., Muralidharan, O., Najmi, A., & Naidu, S. 2012, Estimating Uncertainty for Massive Data Streams, https://research.google/pubs/pub43157/
  * Cutler & Flanagan (1994) Cutler, C., & Flanagan, E. E. 1994, Phys. Rev. D, 49, 2658, doi: 10.1103/PhysRevD.49.2658
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M., & Rubin, D. B. 1977, Journal of the Royal Statistical Society. Series B (Methodological), 39, 1.  http://www.jstor.org/stable/2984875
  * Dominik et al. (2015) Dominik, M., Berti, E., O’Shaughnessy, R., et al. 2015, Astrophys. J., 806, 263, doi: 10.1088/0004-637X/806/2/263
  * Edelman et al. (2021) Edelman, B., Doctor, Z., Godfrey, J., & Farr, B. 2021.  https://arxiv.org/abs/2109.06137
  * Edelman et al. (2023) Edelman, B., Farr, B., & Doctor, Z. 2023, Astrophys. J., 946, 16, doi: 10.3847/1538-4357/acb5ed
  * Efron (1979) Efron, B. 1979, The Annals of Statistics, 7, 1 .  https://doi.org/10.1214/aos/1176344552
  * Farah et al. (2023) Farah, A. M., Edelman, B., Zevin, M., et al. 2023.  https://arxiv.org/abs/2301.00834
  * Fishbach & Holz (2017) Fishbach, M., & Holz, D. E. 2017, Astrophys. J. Lett., 851, L25, doi: 10.3847/2041-8213/aa9bf6
  * Fishbach & Holz (2020) —. 2020, Astrophys. J. Lett., 891, L27, doi: 10.3847/2041-8213/ab7247
  * Kimball et al. (2021) Kimball, C., et al. 2021, Astrophys. J. Lett., 915, L35, doi: 10.3847/2041-8213/ac0aef
  * Kovetz et al. (2017) Kovetz, E. D., Cholis, I., Breysse, P. C., & Kamionkowski, M. 2017, Phys. Rev. D, 95, 103010, doi: 10.1103/PhysRevD.95.103010
  * LVK (2021a) LVK. 2021a, GWTC-3: Parameter estimation data release, Zenodo, doi: 10.5281/zenodo.5546663
  * LVK (2021b) —. 2021b, GWTC-3: O1+O2+O3 Search Sensitivity Estimates, Zenodo, doi: 10.5281/zenodo.5636816
  * Mandel (2010) Mandel, I. 2010, Phys. Rev. D, 81, 084029, doi: 10.1103/PhysRevD.81.084029
  * Mandel & de Mink (2016) Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634, doi: 10.1093/mnras/stw379
  * Miller et al. (2020) Miller, S., Callister, T. A., & Farr, W. 2020, Astrophys. J., 895, 128, doi: 10.3847/1538-4357/ab80c0
  * Ng et al. (2018) Ng, K. K. Y., Vitale, S., Zimmerman, A., et al. 2018, Phys. Rev. D, 98, 083007, doi: 10.1103/PhysRevD.98.083007
  * O’Leary et al. (2016) O’Leary, R. M., Meiron, Y., & Kocsis, B. 2016, Astrophys. J. Lett., 824, L12, doi: 10.3847/2041-8205/824/1/L12
  * Olejak et al. (2021) Olejak, A., Belczynski, K., & Ivanova, N. 2021, A&A, 651, A100, doi: 10.1051/0004-6361/202140520
  * Powell et al. (2019) Powell, J., Stevenson, S., Mandel, I., & Tino, P. 2019, Mon. Not. Roy. Astron. Soc., 488, 3810, doi: 10.1093/mnras/stz1938
  * Ray et al. (2023) Ray, A., Magaña Hernandez, I., Mohite, S., Creighton, J., & Kapadia, S. 2023\.  https://arxiv.org/abs/2304.08046
  * Rinaldi & Del Pozzo (2021) Rinaldi, S., & Del Pozzo, W. 2021.  https://arxiv.org/abs/2109.05960
  * Rodriguez et al. (2016) Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, Phys. Rev. D, 93, 084029, doi: 10.1103/PhysRevD.93.084029
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 100, 043027, doi: 10.1103/PhysRevD.100.043027
  * Sadiq et al. (2022) Sadiq, J., Dent, T., & Wysocki, D. 2022, Phys. Rev. D, 105, 123014, doi: 10.1103/PhysRevD.105.123014
  * Sain & Scott (1996) Sain, S. R., & Scott, D. W. 1996, J. Am. Stat. Assoc., 91, 1525.  https://www.tandfonline.com/doi/abs/10.1080/01621459.1996.10476720
  * Schneider et al. (2023) Schneider, F. R. N., Podsiadlowski, P., & Laplace, E. 2023, Astrophys. J. Lett., 950, L9, doi: 10.3847/2041-8213/acd77a
  * Talbot & Thrane (2018) Talbot, C., & Thrane, E. 2018, Astrophys. J., 856, 173, doi: 10.3847/1538-4357/aab34c
  * Terrell & Scott (1992) Terrell, G. R., & Scott, D. W. 1992, The Annals of Statistics, 20, 1236 .  https://doi.org/10.1214/aos/1176348768
  * Thrane & Talbot (2019) Thrane, E., & Talbot, C. 2019, Publ. Astron. Soc. Austral., 36, e010, doi: 10.1017/pasa.2019.2
  * Tiwari (2021) Tiwari, V. 2021, Class. Quant. Grav., 38, 155007, doi: 10.1088/1361-6382/ac0b54
  * Tiwari (2022) —. 2022, Astrophys. J., 928, 155, doi: 10.3847/1538-4357/ac589a
  * Tiwari (2023) —. 2023.  https://arxiv.org/abs/2304.03498
  * Tiwari & Fairhurst (2021) Tiwari, V., & Fairhurst, S. 2021, Astrophys. J. Lett., 913, L19, doi: 10.3847/2041-8213/abfbe7
  * Toubiana et al. (2023) Toubiana, A., Katz, M. L., & Gair, J. R. 2023.  https://arxiv.org/abs/2305.08909
  * van Son et al. (2022) van Son, L. A. C., de Mink, S. E., Renzo, M., et al. 2022, ApJ, 940, 184, doi: 10.3847/1538-4357/ac9b0a
  * Veitch et al. (2015) Veitch, J., et al. 2015, Phys. Rev. D, 91, 042003, doi: 10.1103/PhysRevD.91.042003
  * Veske et al. (2021) Veske, D., Bartos, I., Márka, Z., & Márka, S. 2021.  https://arxiv.org/abs/2105.13983
  * Wang & Wang (2011) Wang, B., & Wang, X. 2011.  https://arxiv.org/abs/0709.1616
  * Wysocki (2020) Wysocki, D. 2020, Calibrating semi-analytic VT’s to injections in O3a, https://dcc.ligo.org/LIGO-T2000432/public
  * Wysocki et al. (2019) Wysocki, D., Lange, J., & O’Shaughnessy, R. 2019, Physical Review D, 100.  http://dx.doi.org/10.1103/PhysRevD.100.043012
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152, doi: 10.3847/1538-4357/abe40e
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, Astrophys. J., 910, 152, doi: 10.3847/1538-4357/abe40e

##  Appendix A Two-dimensional mock data test of iterative reweighting

Figure 9:  Left: The KDE estimate from data with error for an anticorrelation
is added in it. The true distribution is a mixture of a uniform distribution
and a bivariate Gaussian with equal and uncorrelated variances in the two
dimensions. We then simulate a measurement error which is strongly
anticorrelated between the two dimensions to generate mock data parameter
samples (see main text for details). We see the contours around (50, 50) is
wider (elliptical) and clearly show an anti-correlation. Right: The result
after reweighting iteration are applied. We see iterative reweighting reduce
the artificial correlation and we do not see elliptical regions around the
peaks and results are close to our true distribution.

We here demonstrate iterative reweighted KDE over 2D mock data sets. We start
with 60 mock events, 50% from a two-dimensional uniform distribution and 50%
from a bi-variate normal (Gaussian) distribution. The normal distribution mean
is $\mu=(50,50)$ (in arbitrary units) and we take the two dimensions to be
uncorrelated each with a variance of $9$. We added error in our mock data
values with an anti-correlation between the two dimensions using covariance
$\Sigma_{p}=\bigl{(}\begin{smallmatrix}32.5&-31.5\\\
-31.5&32.5\end{smallmatrix}\bigr{)}$, corresponding to an error ellipse with
8:1 axes at $45^{\circ}$, to mimic the anticorrelation between $m_{1}$ and
$m_{2}$ along a contour of constant chirp mass. 100 mock parameter samples
with the same error covariance are then generated around each measured value.
We construct KDE on our errored data points and apply our reweighting scheme.
The main idea we want to check if our scheme can reduce this artificial
correlation we added in true distribution. The left panel of Fig. 9 is the KDE
estimate on errored data before applying iterative reweighting. We can clearly
see an anti-correlation we added into our true distribution around the peaks
$(50,50)$. The right panel of Fig 9 is the results after we apply our
iterative reweighting scheme and use the last 900 iterations to get median
curves. We can see that anti-correlation is removed and our method recovered
the true distribution. This example shows that our method can reduce effect of
artificial correlation in the real data. We also compute 1-d results from
these 2-d iterative reweighted average results and compared this corresponding
true distribution in 1-d that we expect from true 2-d distribution. Note that
our 2-d distribution can be a union of two 1-d uniform distribution with a
Gaussian peak $\mu=50$. Fig. 10 shows the true distribution in 1d (black
dashed) where the red curves are results of numerical integration of iterative
reweighted KDE estimates for the last 900 iterations. The median (solid) red
curve matched very closely to the true distribution and recover the peak
height and width of Gaussian peak in true distribution.

Figure 10: Integrated 1-d component distribution derived from 2-d mock data
results. The black dashed curves is the true 1-d distribution composed of
uniform and Gaussian components with a Gaussian mean $\mu=50$. The red curved
are x coordinate from numerical integration of two dimensional mock data
results (as shown on right of Fig. 9) after iterative reweighting is done.

##  Appendix B Results including GW190814

Figure 11: Left: Two dimensional rate estimated with random PE samples for ln
$m_{1}$ and ln $m_{2}$ of 70 events (69BBH) and including GW190814 from
GWTC-3. The red plus symbols show the median values $m_{1}$ and $m_{2}$ of
samples of each event. Right: The corresponding one dimensional component mass
distributions from numerical integration of two dimensional results.

In addition to our main results from significant BBH events in GWTC-3, we
examine the effect of including the outlier event GW190814 (Abbott et al.,
2020b) with low secondary mass $m_{2}=2.59\,M_{\odot}$ $(q=0.112)$. We
slightly extend the range of $m_{2}$ over which the KDE is evaluated in order
to include the additional event. Fig. 11 summarizes our results, which are
also briefly discussed in the main text.

  *[BBH]: Binary black hole
  *[GW]: gravitational-wave


# Blockchain for Access Control in e-Health Scenarios

João Pedro Dias  INESC TEC and
Faculty of Engineering
University of Porto
Porto, Portugal
<EMAIL_ADDRESS>Luís Reis  INESC TEC
FEUP Campus
Rua Dr. Roberto Frias
Porto, Portugal
<EMAIL_ADDRESS>Hugo Sereno Ferreira  INESC TEC and
Faculty of Engineering
University of Porto
Porto, Portugal
<EMAIL_ADDRESS>Ângelo Martins  INESC TEC
FEUP Campus
Rua Dr. Roberto Frias
Porto, Portugal
<EMAIL_ADDRESS>
###### Abstract

Access control is a crucial part of a system’s security, restricting what
actions users can perform on resources. Therefore, access control is a core
component when dealing with e-Health data and resources, discriminating which
is available for a certain party. We consider that current systems that
attempt to assure the share of policies between facilities are prone to
system’s and network’s faults and do not assure the integrity of policies
lifecycle. By approaching this problem with the use of a distributed ledger,
namely a consortium blockchain, where the operations are stored as
transactions, we ensure that the different facilities have knowledge about all
the parties that can act over the e-Health resources while maintaining
integrity, auditability, authenticity, and scalability.

###### Index Terms:

e-Health, Access Control, Blockchain, Distributed Ledger Technology, Security,
Distributed Systems

##  1 INTRODUCTION

Healthcare is experiencing an explosion of data partially due to the
widespread of health data collection systems such as wearables (e.g. fitness
trackers) [1], health tracking applications (e.g. diet tracking) [2] and
ambient assisted living systems such as CAALYX [3]. By now, it is estimated
that medical data will grow at a rate of 48% per year, reaching 2.3 zettabytes
by the year of 2020 [4].

Lots of new smart objects are empowering the creation of cyber-physical smart
pervasive systems with application in multiple domains, including healthcare.
These smart objects are under the umbrella of the Internet of Things paradigm,
that foresees the advance towards new smart and inter-connected systems by the
means of ubiquitous computing [5].

The explosion of data being collected and, a posteriori, analyzed by different
entities, leads to the debut of data security and privacy issues. By one hand,
these issues are taken into account because such smart devices may be
connected to the Internet at some point for accessing its collected data
anytime and anywhere [5]. The data being collected from those devices may be
part of the Personal Health Records (PHR) and this is typically owned by the
patients and may be or not, shared with third-entities [6]. On the other hand,
Electronic Medical Records (EMR) and Electronic Health Records (EHR) store
individual information that is required by the healthcare professionals and
may be shared among different institutions and facilities [6].

Hence, there is the need to control the accesses to this data resources by
third-entities. Access control is concerned about determining the allowed
activities of certain users, mediating every attempt by a user to access a
resource in the system [7]. Dealing with the user access control to health
data, personal or medical, holden by different parties, which may be required
to be accessed by third-parties with different goals (e.g. insurance companies
versus doctors), is not an easy task. This is especially problematic when we
are still moving towards a unified and interoperable electronic health
(e-Health) systems [8].

In this paper, we suggest an approach to the problem of access control in
large scale and distributed systems, as it is observed in e-Health scenarios
where different entities and users should be able to access data with
different permission levels and granularities. The Data Keepers should be able
to manage the accesses to their data by the means of adding, changing or
revoking permissions. Such system should be also capable of defining fine-
grained permissions both, at the user level and, at the resource level.

The system must also be fault tolerant, which means that it must not be
dependable on a centralized architecture. Upon these considerations:
consistency, integrity and authenticity of the operations among nodes should
be assured. The system must be also immutable providing an accurate audit
trail.

This paper is structured as follows: firstly, it is given an overview over the
related work in the scope of permission management in e-Health systems,
focusing also blockchain approaches for access control. Afterwards is it given
a description of the purposed solution architecture. Then we address some core
details of the proof-of-concept implementation. Finally, some final remarks
are presented, summing up the contributions as well as pointing out further
developments.

##  2 BACKGROUND & STATE OF THE ART

###  2.1 BLOCKCHAIN

In our solution, we take into consideration a Distributed Ledger Technology
(DLT), specifically, blockchain. A distributed ledger (also known as shared
ledger) consists of a consensus of replicated, shared and synchronized digital
data distributed along a set of nodes, working as a distributed database,
generally geographically dispersed [9]. It is important to note that, despite
all blockchains being distributed ledgers, not all distributed ledgers are
blockchains.

The blockchain is a specific type of distributed ledger conceptualized by
Satoshi Nakamoto and used as a core component of the digital currency Bitcoin
[10]. Data in a blockchain should be tamper-proof, specifically accomplished
by the use of cryptography, by the means of digital signatures and digital
fingerprints (hashing) [9, 10]. Also, consensus must be assured among peers
considering scenarios where some of the peers are providing erroneous data, by
partially or completed computer/network failures or, even, by malicious intent
when some party tries to subvert the ledger [10].

A blockchain consists of a chain of blocks that contains information about
transactions. Each one of these transactions, is digitally signed by the
entity emitting them. Transactions are combined into a block, that is
committed to the chain, establishing the blockchain. Each block contains the
hash of the previous block, being this propagated along the chain until the
first block, created when the blockchain was firstly created, designed genesis
block [10].

We can then consider that a blockchain works as a state transaction system
(state machine), where there is a state that corresponds to the snapshot of
the chain (the result of all transaction until now) and, after adding a new
block of transactions to the chain, we got a new snapshot that corresponds to
a new state of the system, as result of the new transactions [11].

In order to a validate a block, it is necessary a proof-of-work. This
mechanism is used in order to get a consensus in the peer-to-peer network
[10]. In Bitcoin an HashCash proof-of-work is used, being the work effort
called mining. The mining consists of finding a nonce (by the means of brute-
force) that satisfies the condition of generating a digest with the required
number of leading zeroes. This proof-of-work guarantees consensus in a network
following the principle that the nodes will always accept the longest
available chain [10]. This also implies that older blocks - those further back
in the blockchain - are more secure than newer ones.

There are alternatives to proof-of-work. In the proof-of-stake, as it is being
considered to be used in Ethereum [11], the creator of the next block to be
pushed in the chain is chosen in a deterministic way based on the wealth of
the node [9]. Another one, as used in the Sawtooth Lake [12], uses a Proof of
Elapsed Time (PoET), which is a lottery-based consensus protocol that takes
advantage of the trusted execution environments provided by Intel’s Software
Guard Extensions.

Notwithstanding the common use of blockchain for trade currencies, like
Bitcoin, there exists an array of other applications for the technology. This
is possible because, as blockchain is used to store coin transactions but it
can be used to store any other domain transactions. Furthermore, it can be
used as a general-purpose database distributed system, therefore make it
useful in a large variety of situations [13].

Lastly, blockchains can be considered of three main kinds, as stated by
Buterin [14], namely: public, fully-private and consortium. Public blockchains
(e.g. Bitcoin), is a type of blockchain in which anyone can read, send
transactions to and expect to see them included if they are valid, and,
further, anyone in the world can participate in the consensus process. Fully-
private blockchains consists of blockchains where write permissions are kept
centralized to one organization (even if spread among facilities), existing a
closed group of known participants (e.g. a supply chain) [13]. Finally,
consortium blockchains, are partly private in such way that the consensus
process is controlled by a number of pre-selected set of nodes. In this type
of blockchain, the right to query the blockchain may be public, or restricted
to the participants (e.g. governmental institutions and partners).

###  2.2 ACCESS CONTROL

The problem of access control has already been covered in the literature. We
can observe different ways of controlling and managing accesses in different
situations in our everyday technological systems. There exists a problem of
defining permission rules, typically known as policies, alongside with the
problems related with inconsistency, especially in e-Health systems [15].

One of the more common approaches is the use of Access Control Lists (ACL),
commonly used in modern operating systems. ACL consists of a list associated
with an object that specifies all the subjects that can access it, along with
the access level (or rights) [7].

Other systems use Access Control Matrix, in which, each row represents a
subject, each column an object and each entry is the set of access rights for
that subject to that object [7].

Specifically in healthcare, Role-Based Access Control (RBAC) and Attribute-
Based Access Control (ABAC) have been applied [16]. RBAC defines the user’s
access right basing itself on his/her roles and the role-specific privileges
associated with them. The ABAC system extends the RBAC role-based features to
attributes, such as properties of the resource, entities, and the environment
[17]. Policies in ABAC can be expressed resorting to the eXtensible Access
Control Markup Language (XACML), defined by the OASIS consortium [18]. The
XACML standard also includes a reference architecture for designing and
implement access control systems, defining the system components and usage
flow, that can be used in multiple application domains.

Another approach to access control is the Entity-Based Access Control (EBAC)
[19], which allows the definition of more expressive access control policies.
This is accomplished by supporting both, the comparison of attribute values,
as well as traversing relationships along arbitrary entities. Moreover,
Bogaerts et al. presents Auctoritas as an authorization system that specifies
a practical policy language and evaluation engine for the EBAC system [19].

###  2.3 BLOCKCHAIN APPLIED TO ACCESS CONTROL

With the widespread use of DLT, and, more specifically, blockchain, novel
approaches to resolve the access control problematic appeared using this
paradigm, including in e-Health scenarios. Maesa et al. [20] proposes a
blockchain-based access control, implementing ABAC on top of the blockchain
technology, following the XACML reference architecture. This approach
validates itself through a reference implementation on top of Bitcoin. However
this solution does not encompass the particularities of using such in e-Health
context, namely, the possibility of having different authorities and/or
entities as resource owners.

In the application of blockchain for e-Health, Yue et al. [21] proposes the
use of a Healthcare Data Gateway (HGD) to enable the patient to own, control
and share their data while maintaining data privacy. This solution also
encompasses that all the patient e-Health record is stored in a blockchain.
Although the novelty of such approach, it implies a disruptive change on the
already-existent systems of storing and retrieving e-Health data, what would
require a considerable effort to implement which may call into question its
current applicability. Also, there are cases where it is needed to access data
without the explicit agreement from the patient itself (e.g. due to the
patient inability to allow the access or by some governmental requirements)
and this solution does not provide the ability to do such (e.g. some family
member allow the data access). Also, considering the growth of e-Health data,
storing this data on the blockchain itself will result in a rapid growth on
its size, exceeding publicly available hard drive capacity, requiring special
hardware to full nodes and, further, could lead to the centralization of the
blockchain [22].

##  3 ILLUSTRATIVE SCENARIO

As of today, there exists an increment of the e-Health data being produced at
different sources. Such data can come from more traditional origins like
medical exams or medical reports, but, with the advent of the Internet-of-
Things, things like wearables (e.g. fitness trackers) and ambient assisted
living systems, even more data is being produced and consumed by the
individual and/or by 3rd-parties. However, it is hard to keep track on the
localization of such data as well as when and who is consuming that data. As
such, there must exist an access control system transversal to the e-Health
domain that allows one to keep data ownership, managing easily the access
policies in place. Such system must be clear to the user, as it can easily
authorize, deny or revoke permissions on-the-fly (e.g. by the means of a
mobile application and push notifications).

Figure 1: An example of an access to e-Health data, including its different
participants. In detail, a third-party entity can request access to a given
document, and, in order to do so, the access must be allowed by a request
handler system. Such handler should then verify the request by checking an
access control policies repository and, further, it must grant or deny the
access. If some previously unknown request appears, Data Keepers must issue a
new access control policy.

An illustrative scenario for observing how an access control system is
required to work and, how it would impact the e-Health workflow a case study
was designed. A draft of such scene is given on Figure 1. Here we can define
the two major actors of such scenario, Data Keepwers and Third-Parties.

Data Keepers is the set of individuals or entities that have ownership over a
certain data-entry. They can be, but are not limited to, the data creators,
institution responsible, or the data subject(s). The possibility of defining
different sets of keepers is a requirement since, depending on the e-Health
record or source, different restrictions may apply. Namely, when dealing with
PHR data the only keeper is the individual or the legal responsible [23].
However, EMR data is controlled by a set of clinicians and staff within one
health care organization, and, further, EHR data can even be shared between
more than one health organization [23].

Third-parties are any organization, entity or individual that have interest in
accessing the data (e.g. insurance companies or research institutes).
Generally, these consumers have time-limited access to a small portion of data
of a dataset or to a specific individual.

In detail, we can sum up the interaction of the different parties and the
system as it follows (using the Figure 1 as reference):

  1. 1.

There is a 3rd Party who have interest in viewing a specific e-Health record
(in Figure 1 as example we use the document id #123). This party asks to an
intermediate if it can access the record mentioned;

  2. 2.

The intermediate checks the access against an Access Policy repository where
it matches with already existent policies. If there is no Access Policy within
a particular third-entity and a document, an extra two steps occurs: (1) A
notification with the details of the access request is sent to the set of
document keepers; (2) The access keepers must allow or deny the access
request. If there is more than one keeper on the document, a consensus must be
reached. This consensus depends on the proprieties of the document and it can
require the approval of one keeper, majority or all.

  3. 3.

After checking the request against the Access Policies repository, the
intermediate grants or denies access to the document by the 3rd-party;

  4. 4.

If access is allowed, the 3rd Party can now view the specific e-Health record
(in Figure 1 the document id #123). If not, the request is denied and the
entity cannot access the document in any way.

Taking into account the interaction flow we can say that the scenario shares
the three main components of common access control systems (specified in XACML
standard), namely, a Policy Decision Point (PDP), a Policy Administration
Point (PAP) and a Policy Enforcement Point (PEP) [18].

Here, the PEP is the intermediate negotiator that intercepts the users’
requests and enforces the PDP’s decision. Further, the PDP is responsible for
evaluating the requests against an Access Policies repository. However, in our
case, there is not a traditional PAP, because there is no central authority
managing the policies, being each one of the document keepers managers of
their own documents, working as a distributed PAP system.

In this paper, we leverage the use of blockchain as a repository for Access
Policies, and, furthermore, as a way of enabling the existence of a
distributed PAP component.

##  4 PROPOSED APPROACH

Our approach consists of using Blockchain technology as a way to accomplish a
more reliable and user-empowered solution for access control management in an
e-Health environment. Such approach allows us to define fine-grained access
control while maintaining the consensus in a distributed system, authenticity,
immutability and auditability.

A proof-of-concept of the approach hereby described and detailed was
implemented in order to verify its feasibility.

In our solution, we use an approach similar to the Access Control Matrix,
which allows the establishment of a correspondence between a subject, an
object and a set of rights. However, this information is not stored as is, due
to the inherent proprieties of the use of blockchain. As a transaction-based
state machine, we store transactions corresponding to a pre-defined set of the
state machine on the Access Policies repository.

###  4.1 Access Control Model

Upfront for defining our model, it is needed to define all the entities and
relationships enrolled in it. Such model is given in Figure 2, and five
classes can be identified in our approach.

Figure 2: UML diagram specifying the classes of the system and their
relationships.

First off, we have Entities that can be 3rd Party’s or Data Keepers. Further,
we have Policies and Records. Each Policy refers to a relation between one and
only one 3rd Party and an e-Health Record, with the respective level of
access, PermissionLevel. In turn, each Record can have one or more Data
Keepers that have partial or total ownership over it.

Figure 3: State machine diagram representing the lifecycle of an access
request by a 3rd Party. The composite state AUTH_CHECK represents the
individual authorization requests needed by the Data Keepers of the record
being queried. Figure 4: State machine diagram representing the individual
authorization lifecycle.

Moreover, the information related to the model is preserved by the means of
storing transactions (since this model is compatible with the use of
blockchain). The transactions contain information about with 3 different state
machines, that have dependencies between them, being always related to the
class model defined.

Access Policy State Machine, as represented in Figure 3, is the state machine
related to the main logic of creating Access Policies.

Record Life-cycle State Machine is responsible for the operations over an
e-Health record. The record lifecycle begins with its creation, CREATE, then
it can suffer diverse updates, UPDATE, until it is removed, REMOVE.

Individual Authorization State Machine, as represented in Figure 4, describes
the lifecycle for each user access grants over a given REQUIRE, which lead to
a number of instantiations equal to the number of Data Keepers required. Each
individual instantiation evokes a REQUIRE_ACTION, then the Data Keeper can
allow (AUTH_GRANT) or deny (AUTH_DENY) the access. Eventually, before reaching
the final state, the Data Keeper can revoke (AUTH_REVOKE) a previously granted
authorization.

As stated, the main logic of attributing Access Policies for 3rd Party access
an e-Health record is controlled by the Access Policy State Machine, detailed
in Figure 3. This state machine jumps the init state when an access request,
REQUEST from a 3rd Party enters the PEP. Then, the REQUEST is verified,
VERIFY, against the already existent information on the blockchain (by the
means of a snapshot operation). If, and only if, the information about this
particular access is present in the snapshot, the request can be granted,
AUTH_GRANT or denied AUTH_DENY.

Additionally, if there is no information about an access request, the Access
Policy must be required, REQUIRE, by the means of checking the number of
permissions required form the textttData Keepers in order to get a consensus
(specific to the record), using for that the Individual Authorization State
Machine. This means that there is no central authority authorizing requests
from 3rd Parties, and it is required that some set of Data Keepers allow the
access. While this process is running, the state machine enters into a waiting
state, WAITING_AUTH_CHECK.

At last, there can be a point in the future when it is needed to revoke a
previously granted access, AUTH_REVOKE. The final state is, by this, reached
by the existence of an AUTH_REVOKE or an AUTH_DENY.

###  4.2 Block Model

As our Access Control follows the model presented in Subsection 4.1, based
upon transactions, we leverage the use of blockchain as a repository for these
transactions. This way we never have a current state of the permissions
written in the repository, but, notwithstanding, we can take a snapshot to the
blockchain at any moment and, as result, check the Access Policies in place.

The basic structure of a block in our chain is the following:

  * •

index: Corresponds to the index of the current block in the blockchain.

  * •

timestamp: Timestamp corresponding to when the block was generated.

  * •

previousHash: Hash of the previous block in the chain.

  * •

digitalSign: Digital signature of the current block data.

  * •

data: Content of the block. Corresponds to a set of transactions describing
the access control policies, records information and individual
authorizations.

  * •

nonce: Value that is set so that the hash of the block will contain a run of
leading zeros. This value is calculated iteratively until the required number
of zero bits in the hash is found. This requires time and resources, making it
so that the correct nonce value constitutes proof-of-work.

  * •

hash: A SHA256 hash corresponding to the block data. This hash must have a
leading a priori defined sequence being this leading sequence what defines the
effort of the proof-of-work. In Bitcoin, this leading sequence corresponds to
a certain number of zeros in the beginning of the hash.

Additionally, focusing on our approach, the data field should be detailed, as
it is the field that serves as transaction information storage. This data
field includes in it three sub-fields, namely:

  * •

records: Transaction information relative to transactions of the state machine
presented in Figure LABEL:fig:doc-sm, about creation, update or deletion of
e-Health records of any kind.

  * •

policies: Transaction information relative to transactions of the state
machine presented in Figure 3, about creation and revocation of access
policies.

  * •

individualAuths: Transactions about individual authorization by each one of a
Record Data Keeper in relation to each Policy.

The use of hashes allows us to maintain integrity along the immutable chain of
transactions without a central authority, since any change in the data would
result in a different hash, invalidating the next blocks in the blockchain.
Additionally, as result, we can also achieve accountability and auditability.
Authenticity is assured by the assign of a key-pair to each entity with access
to the blockchain, identifying who write each block in the chain.

###  4.3 Architectural Design

Figure 5: Overview of the system architecture and interactions. The e-Health
data access requests treated seamlessly by the handler. Additionally, the
resources (e-Health data) are, typically, associated with one facility or
repository, but all the available resources are known by all nodes
(facilities).

At the system architectural level, our approach uses blockchain, being
distributed by default, working as a peer-to-peer network connecting the
different nodes, corresponding to the diverse facilities or organizations that
can store, create or/and change e-Health data. Nodes in the network
synchronize between them by following a set of rules:

  * •

When a new block is generated by a node, this block is broadcasted to the
network;

  * •

When a node connects to a new peer in the network it queries for the latest
block;

  * •

If a node finds a block that has a higher index than the last known block, it
either adds the block to its current blockchain (in the case of the difference
is equal to one node) or queries another node to get the full blockchain.

The system implements a consortium blockchain[14], which means that the
blockchain is partly private, in the way that instead of allowing any person
to participate in the process of transaction verification or, allowing only
one entity to have full control, a few selected nodes are predetermined,
providing the same benefits associated with private blockchain. The Data
Keepers and 3rd Party can then interact with the system by the means of using
public available Application Programming Interfaces (API’s) or applications
designed to do such.

Furthermore, despite the use of XACML standard for access control systems and
blockchain as storage, there are some architectural choices due to the
e-Health domain restrictions and more complex use cases. An overview can be
observed in Figure 5.

Focusing on the e-Health Data, as of today, this data is not aggregated in one
storage, being spread by multiple institutions and organizations. As such,
every time that any 3rd Party requests access to a specific record there is
the need of locating this information, and then, proceed to check if the
request is already approved or if there is the need to create a new access
policy. So, as an improvement, information about the creation of new records
must be kept and spread along all the organizations and facilities in such way
that a request to a resource can be handled by any member of the private
blockchain.

Additionally, aiming attention to the Data Keepers, it was noticed that there
is a set of situations where the ownership of e-Health data records is not
explicit to only one entity but shared among more than one entity or
individual, as is the case of EMR versus PHR. Taking this into account, we set
up a mechanism of consensus when creating new access policies. Each e-Health
record has a level of agreement that must be achieve before allowing a 3rd
Party the access to a Resource, being this level associated with the Resource
itself (Figure 2). Then, there is a number of executions of the individual
authorization state machine (Figure 4) corresponding to the number of that
Record Data Keepers. Reaching the minimum number of individual authorization
(that can be either AUTH_GRANT or AUTH_DENY), the access request state machine
(Figure 3) will create an access policy accordantly with the consensus reached
(that can be, once more, either AUTH_GRANT or AUTH_DENY).

Figure 6: Sequential view on how some 3rd Party can access or request access
to an e-Health resource, detailing the communication between the inner modules
of the architecture. It is also visible the process of creation of new access
control policies.

From the functional architecture viewpoint, we can sum up the system
interaction as stated in the sequence diagram of the Figure 6. This diagram
describes the process of a 3rd Party requesting access to an e-Health Record,
owned partially or totally by one or more Data Keepers. Further, the diagram
describes both the case of checking against an already existent access policy
of the 3rd Party and the Record or the process of creating a new Policy by
checking the necessary number of Data Keepers.

###  4.4 Security Threat Analysis

In our approach, there are a number of security questions that must be taken
into consideration since we have to consider that mistakes can be made by, for
example, careless operators. Alongside with threats coming from human
mistakes, there is the need of considering also faults introduced by
intentional system manipulation coming from individuals with malicious intent.
Note that we assume that system’s operators do not have physical access to the
machines where the Access Control system is deployed, communicating only by an
existent GUI/CLI.

#### Why existing Access Control system storage do not suffice?

Currently deployed systems use a centralized or almost-centralized (quasi-
decentralized) solution to store and access control policies, independently on
how these policies are described (e.g. ACL, RBAC or EBAC). Such systems are
vulnerable to network or machine failures since in case of failure all the
facilities are affected by the impossibility of validating access control
policy rules over their resources. The use of a DLT as a purely-distributed
database for storing the policies allows the correct functioning of all non-
falling facilities and automatically synchronizes when the normal function of
the system is restored. Additionally, even if a node is compromised, the
threat is confined to the node specifically, not compromising the whole
network.

#### Why do Access Control systems have difficulties ensuring the validity of
policies lifecyle?

Access Control systems typical record the addition, revoking or change in a
policy by the means of logging. Such approach is vulnerable to modifications
by a malicious third-party because there is no way of assuring the integrity
of these logs. The use of a blockchain as a way of storing policies as
transactions, assures us that older policies (and operations over them) have
integrity and are frequently recomputed from a set of trusted commands.
Additionally, we can ensure that, as time passes (and the blockchain grows),
the blocks are more secure as you go further back in the chain.

#### What does it take to a third-entity with malicious intent to add a rogue
policy?

As the access to resources present in the e-Health system must be validated by
the means of a policy in the Access Control system, successful violation by a
malicious party is unconditionally dependent on the injection of a rogue
permission into the system. Such attempt can be mainly made by two attack
vectors: (i) by attempting to mine a new block in the blockchain, though our
approach makes that unlikely due to the usage of a consortium blockchain,
where only allowed entities can add blocks; hence, there is the need of such
malicious entity to take over more than 51% of the entire blockchain network
to accomplish such task; and (ii) by compromising an Access Control system
machine or any third-application that communicates with it, allowing a party
to submit one or more rogue policies to the system. This last type of attack
cannot be prevented by the blockchain, and thus relies on good security
practices on edge applications and devices.

#### What are the privacy implications of having a transparent ACL system?

As the data stored in the blockchain is available to all nodes that belong to
it, some privacy concerns can arise. However, since the approach consists of a
consortium blockchain, all the nodes are, a priori, well-known and
trustworthy, mitigating the privacy concerns of a public blockchain.

##  5 IMPLEMENTATION DETAILS

During the implementation of the proof-of-concept, some decisions were made;
we describe these details here with the intent of helping the reader to re-
implement a similar proof-of-concept or production-ready system.

The whole system, from the transactions logic to the writing to the
blockchain, was implemented using Javascript due to the simplicity of the
language and availability of libraries. Such libraries, as for example, the
built-in crypto module provides us with a mechanism to digitally sign the
blocks payload using public-key cryptography (using RSA-SHA256) and to
calculate the hashes of each block using SHA256 algorithm.

From our perspective and due to the huge disparity among data types, three
main classes were created: Document, Entity and Transaction.

The information of the State Machines mentioned in section 4.1 is given by the
Transactions data present in each moment and not by the states of the machine
itself. Thus, the current state of the system given by the process of
snapshot. A snapshot is like a picture taken from the blockchain, consisting
of applying all previous transactions until now in order to get the state of
all operations. As example, for evaluate a access request requires the
execution of all transactions about that specific 3rd Party, Record and
associated Data Keepers.

The implemented blockchain uses proof-of-work that bases itself on a brute-
force mechanism of hash calculation. This mechanism works as the nonce is
iteratively incremented until the resulting SHA256 hash matches the a priori
defined number of leading zeroes — this is similar to the Bitcoin system and
establishes the “effort”.

However, we can easily tweak the “effort” to better suit our use-case.
Although we implemented the proof-of-work mechanism, due to its simplicity,
there are alternatives as Proof of Elapsed Time (PoET) or proof-of-stake
mechanisms, that are not as resource consuming as the implemented one, which
can be better to validate the transactions in the context of a consortium
blockchain applied to e-Health scenarios.

One of the details that required some attention was the writing of blocks to
the blockchain. As previously mentioned, each transaction has a unique
identifier (ID) to allow the identification of the same transaction over its
different states. This helps to prevent different states of the same
transaction from being written to the same block.

##  6 CONCLUSIONS

In this paper, we present an approach to solving the problem of managing
access control in the e-Health ecosystem. Access Control is a special complex
task in e-Health since resources and data are distributed among different
facilities and institutions. Further, this is even more problematic because in
some cases, e-Health resources are not owned or managed by a single entity or
individual. As a way of overcoming this complexity, we propose an approach
that leverages the use of blockchain for store transactional information about
e-Health records and access control policies.

For purposes of supporting the plausibility of our solution, a proof-of-
concept was designed and implemented. This proof-of-concept allowed us to make
some, even if preliminary, tests and validations over our approach from a
functional and applicational perspective.

Overall, we determine that our approach is viable, giving diverse advantages
when comparing to the in-place systems. This advantages although, not limited
to, includes the integrity, transparency, and authenticity of the access
control policies in the system, being, this information distributed and
synchronized by all the institutions and organizations that make part of the
consortium.

## ACKNOWLEDGEMENTS

This work was supported by Project “NanoSTIMA: Macro-to-Nano Human Sensing:
Towards Integrated Multimodal Health Monitoring and
Analytics/NORTE-01-0145-FEDER-000016” financed by the North Portugal Regional
Operational Programme (NORTE 2020), under the PORTUGAL 2020 Partnership
Agreement, and through the European Regional Development Fund (ERDF).

## References

  * [1] P. Lukowicz, T. Kirstein, and G. Troster, “Wearable systems for health care applications,” _Methods of Information in Medicine-Methodik der Information in der Medizin_ , vol. 43, no. 3, pp. 232–238, 2004.
  * [2] K. Patrick, W. G. Griswold, F. Raab, and S. S. Intille, “Health and the mobile phone,” _American journal of preventive medicine_ , vol. 35, no. 2, p. 177, 2008.
  * [3] M. N. K. Boulos, A. Rocha, A. Martins, M. E. Vicente, A. Bolz, R. Feld, I. Tchoudovski, M. Braecklein, J. Nelson, G. Ó. Laighin _et al._ , “Caalyx: a new generation of location-based services in healthcare,” _International journal of health geographics_ , vol. 6, no. 1, p. 9, 2007\.
  * [4] IDC, “The digital universe: Driving data growth in healthcare,” EMC Corporation and International Data Corporation, Report, 2014.
  * [5] L. Tan and N. Wang, “Future internet: The internet of things,” _2010 3rd International Conference on Advanced Computer Theory and Engineering(ICACTE)_ , vol. 5, pp. V5–376–V5–380, Aug 2010.
  * [6] P. C. Tang, J. S. Ash, D. W. Bates, J. M. Overhage, and D. Z. Sands, “Personal health records: definitions, benefits, and strategies for overcoming barriers to adoption,” _Journal of the American Medical Informatics Association_ , vol. 13, no. 2, pp. 121–126, 2006.
  * [7] V. C. Hu, D. Ferraiolo, and D. R. Kuhn, _Assessment of access control systems_. US Department of Commerce, National Institute of Standards and Technology, 2006.
  * [8] W. O. Nijeweme-d’Hollosy, L. van Velsen, M. Huygens, and H. Hermens, “Requirements for and barriers towards interoperable ehealth technology in primary care,” _IEEE Internet Computing_ , vol. 19, no. 4, pp. 10–19, July 2015.
  * [9] Deloitte., “Bitcoin, Blockchain & distributed ledgers: Caught between promise and reality,” Centre for the Edge, Australia, Tech. Rep., 2015.
  * [10] S. Nakamoto, “Bitcoin: A Peer-to-Peer Electronic Cash System,” p. 9, 2008.
  * [11] F. Ethereum, “Ethereum: A next-generation smart contract and decentralized application platform,” 2014.
  * [12] I. Corporation, “Sawtooth lake latest documentation,” 2015, (Accessed on 06/02/2017).
  * [13] S. Underwood, “Blockchain beyond bitcoin,” _Commun. ACM_ , vol. 59, no. 11, pp. 15–17, Oct. 2016. [Online]. Available: http://doi.acm.org/10.1145/2994581
  * [14] V. Buterin, “On public and private blockchains,” August 2015, (Accessed on 06/06/2017).
  * [15] R. Krawiec and M. White, “Blockchain: Opportunities for health care,” August 2016, (Accessed on 06/05/2017). [Online]. Available: https://www2.deloitte.com/us/en/pages/public-sector/articles/blockchain-opportunities-for-health-care.html
  * [16] M. Li, S. Yu, K. Ren, and W. Lou, “Securing personal health records in cloud computing: Patient-centric and fine-grained data access control in multi-owner settings,” _Lecture Notes of the Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering_ , vol. 50 LNICST, pp. 89–106, 2010.
  * [17] V. C. Hu, D. Ferraiolo, R. Kuhn, A. R. Friedman, A. J. Lang, M. M. Cogdell, A. Schnitzer, K. Sandlin, R. Miller, K. Scarfone _et al._ , “Guide to attribute based access control (abac) definition and considerations (draft),” _NIST special publication_ , vol. 800, no. 162, 2013.
  * [18] S. Godik and T. Moses, “Oasis extensible access control markup language (xacml),” _OASIS Committee Secification cs-xacml-specification-1.0_ , 2002\.
  * [19] J. Bogaerts, M. Decat, B. Lagaisse, and W. Joosen, “Entity-based access control: Supporting more expressive access control policies,” in _Proceedings of the 31st Annual Computer Security Applications Conference_ , ser. ACSAC 2015. New York, NY, USA: ACM, 2015, pp. 291–300.
  * [20] D. Di Francesco Maesa, P. Mori, and L. Ricci, “Blockchain based access control,” _Distributed Applications and Interoperable Systems: 17th IFIP WG 6.1 International Conference, DAIS 2017, Held as Part of the 12th International Federated Conference on Distributed Computing Techniques, DisCoTec 2017, Neuchâtel, Switzerland, June 19–22, 2017, Proceedings_ , pp. 206–220, 2017.
  * [21] X. Yue, H. Wang, D. Jin, M. Li, and W. Jiang, “Healthcare Data Gateways: Found Healthcare Intelligence on Blockchain with Novel Privacy Risk Control,” _Journal of Medical Systems_ , vol. 40, no. 10, 2016.
  * [22] A. Chepurnoy and D. Meshkov, “On space-scarce economy in blockchain systems,” _IACR Cryptology ePrint Archive_ , vol. 2017, p. 644, 2017. [Online]. Available: http://eprint.iacr.org/2017/644
  * [23] NAHIT, “Report to the office of the national coordinator for health information technology on defining key health information technology terms,” Office of the National Coordinator for Health Information Technology, Tech. Rep., April 2008.