# The cosmic variance of testing general relativity with gravitational-wave
catalogs

Costantino Pacilio  0000-0002-8140-4992<EMAIL_ADDRESS>Dipartimento di Fisica “G. Occhialini”, Universitá degli Studi di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  INFN, Sezione di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  Davide Gerosa
0000-0002-0933-3579 Dipartimento di Fisica “G. Occhialini”, Universitá degli
Studi di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  INFN,
Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  School
of Physics and Astronomy & Institute for Gravitational Wave Astronomy,
University of Birmingham,
Birmingham, B15 2TT, United Kingdom  Swetha Bhagwat  0000-0003-4700-5274
School of Physics and Astronomy & Institute for Gravitational Wave Astronomy,
University of Birmingham,
Birmingham, B15 2TT, United Kingdom

###### Abstract

Combining multiple gravitational-wave observations allows for stringent tests
of general relativity, targeting effects that would otherwise be undetectable
using single-event analyses. We show that the finite size of the observed
catalog induces a significant source of variance. If not appropriately
accounted for, general relativity can be excluded with arbitrarily large
credibility even if it is the underlying theory of gravity. This effect is
generic and entirely analogous to the so-called “cosmic variance” of
cosmology: in essence, we only have one catalog that contains all the events.
We show that the cosmic variance holds for arbitrarily large catalogs and
cannot be suppressed by selecting “golden” observations with large signal-to-
noise ratios. We present a mitigation strategy based on bootstrapping (i.e.
resampling with repetition) that allows assigning uncertainties to one’s
credibility on the targeted test. We demonstrate our findings using both toy
models and real gravitational-wave data. In particular, we quantify the impact
of the cosmic variance on the ringdown properties of black holes using the
latest LIGO/Virgo catalog.

_Introduction_ — Gravitational-wave (GW) detections of binary compact objects
allow for new tests of general relativity (GR) in the strong-field regime [1,
2] adding up to those performed with other experimental and astrophysical
probes [3, 4]. Such tests are limited by the intrinsic challenges of modeling
the strong-field dynamics in theories of gravity beyond GR [5, 6, 7, 8], which
prevents a directed, model-dependent search [9]. In this regime, one primarily
relies on testing the null hypothesis that GR is the underlying theory of
gravity [10].

At the individual-event level, tests of GR have been performed since the very
first GW detection of binary black holes (BHs) [2] and more stringent tests
have since then been reported using the increasing number of detections during
the first three LIGO/Virgo observing runs [11, 10, 12]. Combining multiple
events, or stacking, is key to measure effects that are otherwise undetectable
using single sources.

Existing approaches to stacking can be categorized as (i) multiplication of
the individual likelihoods [13, 14], (ii) multiplication of the individual
Bayes factors [15, 16, 17] and (iii) hierarchical inference [18, 19, 20].
Multiplication of the likelihoods assumes that deviations have the same values
across all the events (e.g., constraints on the mass of the graviton) while
multiplication of the Bayes factors assumes that deviations in multiple events
are uncorrelated (e.g., constraints on additional BH hair) [18]. Both
assumptions are unrealistic and Ref. [19] first proposed hierarchical stacking
as a consistent way of combining observations, similarly to that of
hierarchical Bayesian inference used in GW population studies [21, 22]. In
this context, the consistency of the data with GR can be quantified by
standard metrics such as credible levels and Bayes factors.

Care must be exercised when interpreting the results of tests of GR, as they
can lead to incorrect conclusions in the presence of unmodeled physics (e.g.,
environmental effects [23, 24], eccentricity [25, 26]), systematics in the
waveform templates [27, 28], stealth biases [29], and overlapping signals
[30]. In fact, one could also revert the argument and use tests of GR as a
complementary method to identify the presence of systematics [31].

In this letter, we investigate an additional source of uncertainty when
performing catalog tests of GR, namely the variance originating from the
finite size of the catalog itself. We show that, even if the null hypothesis
is correct, it could be excluded with arbitrarily large credibility from the
posterior of the deviation parameters when combining multiple events. The
issue would be mitigated if one were to repeat the experiment multiple times,
as large deviations would only occur in relatively few repetitions. However,
by definition, we are only going to have one catalog that contains all the
observations. In analogy with the known effect in cosmology [32, 33], we refer
to the considerations above as the “cosmic variance” of strong-field gravity
tests, meaning that we only observe one realization of all the possible
catalogs of sources.

Crucially, our key message is that cosmic variance does not invalidate the use
of null tests of GR, but it must be accounted for when interpreting the
results. First, we show that using Bayes factors provides a more conservative
evidence against violations of the null hypothesis than the corresponding
credible intervals might suggest. Second, we design a mitigation strategy by
assigning uncertainties to credible intervals and Bayes factors. Since one
cannot use multiple realizations, we propose bootstrapping as a partial remedy
[34]. In a nutshell, from the original dataset
$\boldsymbol{d}=\\{d_{1},\dots,d_{N}\\}$ one resamples a new dataset with the
same size $\boldsymbol{d}^{\rm boot}=\\{d^{\rm boot}_{1},\dots,d^{\rm
boot}_{N}\\}$ allowing for repetitions. When resampling $\boldsymbol{d}$ with
replacement, there are $\binom{2N-1}{N}$ distinct combinations and the
probability of obtaining the original dataset is as small as $N!/N^{N}$ [35].

This mimics a set of repeated experiments to study the distribution of the
chosen estimators (Bayes factors or credible intervals), which can then be
used to extract summary statistics (e.g. standard deviation, interquantile
range), thus providing uncertainty estimates.

We focus on hierarchical tests of GR as introduced in Ref. [19] as they
represent the most general case. First, we perform numerical experiments to
show that the cosmic variance holds for arbitrarily large catalogs and it
cannot be mitigated by selecting the observations based on their signal-to-
noise ratio (SNR). Then, we demonstrate the impact of the cosmic variance on
real GW data by reproducing and extending a flagship test of GR. In
particular, we consider the so-called pSEOBNR test [36] which targets
deviations in the dominant frequency and damping time of the ringdown portion
of the signal and was recently applied to the GWTC-3 catalog [12]. We show
that, while the hierarchical analysis of the damping time appears to exclude
GR with high credibility, the corresponding Bayes factor prefers GR and the
bootstrapped distributions have significant support in favor of the null
hypothesis.

_Hierarchical stacking_ — We are interested in testing the null hypothesis
(i.e. GR is the true theory) using a deviation parameter $x$, which is scaled
such that it vanishes when the null hypothesis $\mathcal{H}_{0}$ is satisfied

$\mathcal{H}_{0}:\mathcal{H}\wedge\\{x=0\\}\,,$ (1)

where $\mathcal{H}$ is a broader hypothesis. If the null hypothesis
$\mathcal{H}_{0}$ is inconsistent with the data, we expect deviations $x$ to
spread away from $0$ following unknown patterns that are set by the system
parameters and the nature of the deviations. We model the distribution of $x$
as a normal distribution $\mathcal{N}$ with mean $\mu$ and variance
$\sigma^{2}$

$p_{\rm pop}(x|\mu,\sigma^{2})=\mathcal{N}(x|\mu,\sigma^{2})\,.$ (2)

In terms of these hyper-parameters, the null hypothesis maps to

$\mathcal{H}_{0}:\mathcal{H}\wedge\\{\mu=0,\sigma^{2}=0\\}\,.$ (3)

GR tests are performed by applying hierarchical population inference [21] to
reconstruct the distribution of $\\{\mu,\sigma^{2}\\}$ from the observed
events $\boldsymbol{d}=\\{d_{1},\dots,d_{N}\\}$. The posterior is given by

$p(\mu,\sigma^{2}|\boldsymbol{d})\propto\mathcal{L}(\boldsymbol{d}|\mu,\sigma^{2})\pi(\mu,\sigma^{2})\,,$
(4)

where the hierarchical likelihood

$\mathcal{L}(\boldsymbol{d}|\mu,\sigma^{2})=\prod_{i=1}^{N}\int
dx~{}\mathcal{L}(d_{i}|x)p_{\rm pop}(x|\mu,\sigma^{2})$ (5)

can be expressed in terms of the likelihoods $\mathcal{L}(d_{i}|x)$ of the
individual observations and $\pi(\mu,\sigma^{2})$ models the prior. Equation
(5) assumes that the observations are independent of each other, which would
be violated for, e.g., overlapping events. The expressions above do not
include selection effects [21, 22] because in this context we do not wish to
reconstruct the underlying distribution of $x$ but only constrain its value
using the set of observed sources.

While one could consider more informative ansatzes in place of Eq. (2), Refs.
[19, 20] showed that a Gaussian distribution can identify deviations from the
null hypothesis even when these follow more complex patterns.

The consistency with the null hypothesis can be quantified using the quantile

$Q_{0}=\int_{p\geq p(0,0)}p(\mu,\sigma^{2}|\boldsymbol{d})~{}d\mu
d\sigma^{2}\,,$ (6)

which is defined such that $Q_{0}=0$ ($Q_{0}=1$) indicate full consistency
(full inconsistency).

The Bayes factor $\mathcal{B}$ in favor of the null hypothesis
$\mathcal{H}_{0}$ over the broad hypothesis $\mathcal{H}$ can be estimated
using the Savage-Dickey density ratio [37],

$\mathcal{B}=\frac{p(\mu=0,\sigma^{2}=0|\mathbf{d})}{\pi(\mu=0,\sigma^{2}=0)}\equiv\frac{\mathcal{L}(\mathbf{d}|\mu=0,\sigma^{2}=0)}{\mathcal{Z}}$
(7)

where $\mathcal{L}$ is the hyper-likelihood of Eq. (5) and
$\mathcal{Z}=p(\mathbf{d})$ is the evidence of the data under $\mathcal{H}$.
Bayes factors are often interpreted using Jeffreys’ scale [38], where
$\mathcal{B}\geq 10^{2}$ ($\mathcal{B}\leq 10^{-2}$) denotes “decisive”
evidence in favor of (against) the null hypothesis.

From Eq. (7), the Bayes factor scales as $\mathcal{B}\propto\Delta$, where
$\Delta$ is the prior volume: wide priors favor the null hypothesis, vice-
versa tight priors favor the alternative hypothesis. This implies one can
artificially increase the odds for either of the two competing models by
restricting or enlarging the prior volume [39]. In the following, we fix this
ambiguity by restricting the original prior volume to the $(1-p)$ posterior
credible interval along each axis. For concreteness, the fraction of discarded
posterior samples is set to $p=1.973\times 10^{-9}$, which corresponds to a
$6$-$\sigma$ interval if these were Gaussian distributions. We then rescale
$\mathcal{B}$ by the ratio $\Delta_{\rm new}/\Delta_{\rm old}$ of the
restricted and original prior volumes. The rationale behind our choice is that
$\Delta_{\rm new}$ is just as large to encompass the vast majority of the
posterior support and therefore, the resulting Bayes factor constitutes a
somewhat conservative estimate when testing GR. We denote the resulting Bayes
factor as $\mathcal{B}_{\star}$ to distinguish it from the prior-dependent
expression reported in Eq. (7).

_Cosmic variance_ — If the null hypothesis is correct, one would naively
expect that the posterior for $\mu$ and $\sigma^{2}$ would become sharper
around $\mu=\sigma^{2}=0$ as more events are added to the catalog; vice versa,
it would peak away from zero if the null hypothesis is violated in nature.

It is straightforward to check this expectation using a toy model where $x$
has Gaussian likelihoods for all the events

$\mathcal{L}(d_{i}|x)\propto\mathcal{N}(x|\mu_{{\rm obs},i},\sigma^{2}_{{\rm
obs},i})$ (8)

and errors are homoscedastic, i.e., $\sigma_{{\rm obs},i}=\sigma_{\rm
obs}={\rm const.}$ In the limit of large catalogs $N\gg 1$, Eqs. (4) and (5)
reduce to [20]

$p(\mu,\sigma^{2}|\boldsymbol{d})\approx
p(\mu|\boldsymbol{d})p(\sigma^{2}|\boldsymbol{d})\,,$ (9)

with

$p(\mu|\boldsymbol{d})\propto\mathcal{N}\left(\mu\left|{\rm mean}(\mu_{\rm
obs}),\frac{{\rm var}(\mu_{\rm obs})}{N}\right.\right)$ (10)

and

$p(\sigma^{2}|\boldsymbol{d})\propto\mathcal{N}\left(\sigma^{2}\left|{\rm
var}(\mu_{\rm obs})-\sigma_{\rm obs}^{2},\frac{2\,{\rm var}(\mu_{\rm
obs})^{2}}{N}\right.\right)\,\\!.$ (11)

The true value of $x$ under the null hypothesis is $x_{\rm true}=0$, which
implies the $\mu_{{\rm obs},i}$’s are independently sampled from a normal
distribution

$\mu_{{\rm obs},i}\sim\mathcal{N}(\mu_{{\rm obs},i}|\mu_{\rm
true}=0,\sigma_{\rm obs}^{2})\,,$ (12)

where the variance $\sigma_{\rm obs}^{2}$ accounts for the scatter due to
noise in the detector consistently with the assumption of normal likelihoods
[40]. The central limit theorem implies

$\displaystyle{\rm mean}(\mu_{\rm obs})$
$\displaystyle\sim\mathcal{N}\left({\rm mean}(\mu_{\rm
obs})\left|0,\frac{\sigma_{\rm obs}^{2}}{N}\right.\right)\,,$ (13)
$\displaystyle{\rm var}(\mu_{\rm obs})$
$\displaystyle\sim\mathcal{N}\left({\rm var}(\mu_{\rm obs})\left|\sigma_{\rm
obs}^{2},\frac{2~{}\sigma_{\rm obs}^{2}}{N}\right.\right)\,.$ (14)

Plugging Eq. (14) into Eq. (10) shows that the $\mu$ posterior has variance
$\mathchar 21016\relax\,\sigma_{\rm obs}^{2}/N$ around ${\rm mean}(\mu_{\rm
obs})$ at leading order in $N^{-1}$. By direct comparison with Eq. (13), it
follows that $p(\mu|\boldsymbol{d})$ does not necessarily converge towards
$\mu=0$ as $N$ increases, but there is a chance that the particular draw of
$\\{\mu_{{\rm obs},i}\\}_{i=1}^{N}$ from (12) shifts its peak away from the
true value $\mu_{\rm true}=\sigma^{2}_{\rm true}=0$. A similar conclusion
applies to the recovery of $\sigma^{2}$ by direct comparison of Eq. (11) and
(14). This toy model illustrates how the cosmic variance is due to the
variance associated with the finite size $N$ of the catalog.

In writing Eqs. (10) and (11) we neglected the effects of the prior
$\pi(\mu,\sigma^{2})$, i.e. we have assumed it is uniform and unbounded.
However, the condition that $\sigma^{2}>0$ induces boundary effects in
$p(\sigma^{2}|\boldsymbol{d})$. It follows that the posterior
$p(\mu,\sigma^{2}|\boldsymbol{d})$ lacks a frequentist coverage of credible
intervals, that is, it is not true that $Q_{0}>p$ in a fraction $(1-p)$ of
similar experiments. The consequent difficulty of the statistical
interpretation of $Q_{0}$ was already raised in Ref. [39].

We further highlight the impact of the cosmic variance by considering $1000$
catalogs of $N=10^{4}$ events each, with Gaussian likelihoods for $x$ as per
Eq. (8) and three different choices for the stochastic uncertainties. First,
(i) we consider the case of homoscedastic likelihoods with $\sigma_{\rm
obs}=0.1$. Then, (ii) we assume heteroscedastic Gaussian likelihoods with
$\sigma_{{\rm obs},i}=1/{\rm SNR}_{i}$, where ${\rm SNR}\in[10,1000]$ is a
random variable distributed according to $p({\rm SNR})\propto{\rm SNR}^{-4}$
to mimic the density of SNRs expected from realistic GW detections [41].
Finally, (iii) we isolate from the same heteroscedastic catalogs only the
events with ${\rm SNR}\geq 50$, which mimics a scenario where one performs
tests of GR only on a loud subset of the available GW catalog. For each case,
we draw the maximum-likelihood estimators $\mu_{{\rm obs},i}$ from
$\mathcal{N}(\mu_{{\rm obs},i}|\mu_{\rm true}=0,\sigma_{{\rm obs},i}^{2})$ to
capture noise scattering.

Figure 1: Cumulative distribution function of $Q_{0}$ for three sets of
catalog realizations. We consider Gaussian likelihoods with homoscedastic
errors (solid), heteroscedastic errors (dashed), and heteroscedastic with an
SNR cut (dash-dotted). Orange and blue curves are produced by either assuming
the null hypothesis $\mu_{\rm true}=0=\sigma_{\rm true}^{2}=0$ or drawing
$(\mu_{\rm true},\sigma^{2}_{\rm true})$ from their prior, respectively.

We map the $(\mu,\sigma^{2})$ posterior distribution using the dynesty
implementation of nested sampling [42] with 5000 live points and uniform
priors $\mu\in[-0.9,0.9]$ and $\sigma^{2}\in[0,0.9]$. Figure 1 shows the
resulting coverage of $Q_{0}$. To better highlight the peculiarity of the null
hypothesis, we repeat the experiments without assuming that $\mu_{\rm
true}=\sigma_{\rm true}^{2}=0$, but instead sample $(\mu_{\rm
true},\sigma_{\rm true}^{2})$ from their uniform priors at each catalog
realization. As expected, when catalogs are drawn from uniform priors in
$(\mu_{\rm true},\sigma_{\rm true}^{2})$, the quantile $Q_{0}$ has a
frequentist coverage. Instead, when drawing from the null hypothesis, the
recovered values of $\sigma^{2}$ lie close to the edge of the prior, which
produces an excess of low values of $Q_{0}$, thus pushing its cumulative
distribution to the upper-right portion of Figure 1. That said, while $Q_{0}$
lacks a frequentist interpretation, it nonetheless provides an upper-bound
estimate of the false alarm rate because $Q_{0}>p$ in less than a fraction
$(1-p)$ of the catalog realizations. Figure 1 also shows that restricting to
high-SNR events does not reduce the effect of cosmic variance, in agreement
with our interpretation based on the finite size of the catalog.

_Bootstrapping_ — The cosmic variance can be mitigated by assigning
uncertainties to the chosen estimator. We showcase this idea by selecting a
homoscedastic catalog realization with a high null-hypotesis quantile
$Q_{0}=0.98$. The corresponding Bayes factor is
$\log_{10}\mathcal{B}_{\star}=-0.65$, indicating substantial but not decisive
evidence against the null hypothesis. After resampling for 1000 catalogs via
bootstrap, we find that $Q_{0}>0.77$ and
$\log_{10}\mathcal{B}_{\star}=-0.76^{+1.37}_{-2.61}$ at $90\%$ credibility. In
particular, $\log_{10}\mathcal{B}_{\star}>0$ in $23\%$ of the bootstrapped
catalogs, which is a non-negligible fraction and would suggest great care in
claiming that the measurement provides evidence against the null hypothesis.
This toy model shows that our proposed strategy is robust in mitigating false
positives, even for catalogs with large credible quantiles.

Figure 2: Distribution of Bayes factors $\log_{10}\mathcal{B}_{\star}$ and GR
quantiles $Q_{0}$ over $1000$ bootstrapped realizations of the flagship
pSEOBNR catalog test. Green crosses and red hexagons indicate constraints on
the ringdown frequency $\delta\hat{f}_{220}$ and damping time
$\delta\hat{\tau}_{220}$, respectively. In the marginalized histograms,
colored triangles indicate 90% confidence intervals. Black markers indicate
the values of $\log_{10}\mathcal{B}_{\star}$ and $Q_{0}$ corresponding to the
original catalog; this is just one possible realization among many!

_State-of-the-art application_ — We now apply our findings to a standard test
of strong-field gravity with GWs. We consider the pSEOBNR family [43, 36, 31,
28] of binary BH waveforms, which are obtained by augmenting effective-one-
body templates with free parameters corresponding to fractional deviations in
the quasi-normal modes of the remnant BH. In the spirit of BH spectroscopy
[44, 45], the pSEOBNR scheme has been used in tests of GR by allowing for
deviations $\delta\hat{f}_{220}$ and $\delta\hat{\tau}_{220}$ in the dominant
frequency and damping time respectively [10, 36, 12]. The latest iteration of
these tests [12] uses 10 GW events and indicates a moderate deviation of
$\delta\hat{\tau}_{220}$ from the GR value $\delta\hat{\tau}_{220}=0$. While
insufficient to claim inconsistencies with GR, the authors themselves indicate
this finding deserves further investigation.

In order to illustrate the role of the cosmic variance in the interpretation
of the results, we reproduce the pSEOBNR analysis of Ref. [12] with a
hierarchical stacking of the events. For consistency with Ref. [12], we
recover the $(\mu,\sigma)$ posterior and set uniform priors $\mu\in[-0.9,0.9]$
and $\sigma\in[0,0.9]$, covering a region that is much broader than the
resulting posterior. Quoting median and $90\%$ credibility, we obtain
$\mu=0.02_{-0.04}^{+0.04}$, $\sigma<0.06$ for $\delta\hat{f}_{220}$ and
$\mu=0.13^{+0.13}_{-0.11}$, $\sigma<0.19$ for $\delta\hat{\tau}_{220}$, which
is in agreement with the analysis of Ref. [12]. Using Eq. (6), we quantify the
consistency between GR and the data as $Q_{0}=0.32$ for $\delta\hat{f}_{220}$
and $Q_{0}=0.81$ for $\delta\hat{\tau}_{220}$. As already reported [12], the
$\delta\hat{\tau}_{220}$ quantile indicates a moderate preference against the
null hypothesis. Using Bayes factors, we find
$\log_{10}\mathcal{B}_{\star}=1.49$ for $\delta\hat{f}_{220}$ and
$\log_{10}\mathcal{B}_{\star}=0.70$ for $\delta\hat{\tau}_{220}$. In
particular, we note that in the case of $\delta\hat{\tau}_{220}$, even if
$Q_{0}=0.81$, the Bayes factor favors the null hypothesis.

We assign uncertainties to $Q_{0}$ and $\log_{10}\mathcal{B}_{\star}$ by
generating 1000 bootstrapped catalog realizations. For each of these, we
repeat the hierarchical analysis and extract the corresponding values of
$Q_{0}$ and $\mathcal{B}_{\star}$. The analysis of Ref. [12] uses 10 GW
events, which implies there are $\mathchar 21016\relax\,10^{5}\gg 1000$ [35]
independent realizations and the probability of duplications is consequently
small. Our results are shown in Figure 2. For $\delta\hat{f}_{220}$ we find
that $\log_{10}\mathcal{B}_{\star}=1.45_{-0.83}^{+0.25}$ and $Q_{0}<0.77$ at
90% confidence. For $\delta\hat{\tau}_{220}$, we find
$\log_{10}\mathcal{B}_{\star}=0.62_{-1.19}^{+0.70}$ and $Q_{0}>0.42$.

Our bootstrap procedure returns broad histograms for $Q_{0}$; in particular,
the credible quantile of the null hypothesis for $\delta\hat{\tau}_{220}$ can
be as low as $Q_{0}=0.42$ within the $90\%$ range. Accounting for the cosmic
variance mitigates the significance of the inference performed with the
original observed catalog. Moreover, the distribution of the Bayes factors for
$\delta\hat{\tau}_{220}$ does not signal any substantial evidence against the
null hypothesis at $90\%$ credibility; rather, $83\%$ of the samples have
$\log_{10}\mathcal{B}_{\star}>0$, indicating support for the null hypothesis.

Finally, the correlation between $Q_{0}$ and $\log_{10}\mathcal{B}_{\star}$
shown Figure 2 indicates that Bayes factors provide weaker evidence against
the null hypothesis than the corresponding credible level. In particular,
while there are individual catalog realizations with $Q_{0}\approx 1$, the
corresponding Bayes factors barely meet the threshold for decisive evidence.

_Discussion_ — Combining information from multiple observations is a natural
strategy to strengthen one’s statistical inference on a physical phenomenon.
Testing gravity with GWs is no exception. GR is a fundamental pillar of our
understanding of the Universe and, when the stakes are so high, our confidence
in experimental bounds becomes critical. The interpretation of tests of GR
with GW catalogs depends on both the statistics (e.g. quantiles and Bayes
factors) as well as the techniques (e.g. hierarchical stacking) used to
combine the inferences in favor or against the null hypothesis. Crucially, one
must also quantify the cosmic variance originating from the single realization
of the catalog of GW events at our disposal.

In particular, three key points are worth stressing:

  * (i)

The net effect of the cosmic variance is to soften one’s claim in favor of
violations of GR. After all, the GW catalog we happened to observe might be an
unlucky realization of the null hypothesis.

  * (ii)

The cosmic variance can be quantified by producing multiple mock catalogs. We
propose a data-driven approach that does not rely on assuming a population of
sources but instead resamples the observed catalog with repetition.

  * (iii)

The cosmic variance does not vanish as either the size of the catalog or the
SNR of the events increase.

Points (i) and (ii) are best exemplified on the BH ringdown test we borrowed
from the flagship analysis of Ref. [12]. We show that, while the current
catalog presents a quantile $Q_{0}$ that might be interpreted as a moderate
deviation from GR, this evidence turns out to be insignificant when the
original measurement is considered as a part of a distribution of bootstrapped
estimators. We have illustrated point (iii) with a toy model based on Gaussian
likelihoods.

Our findings lead to the conceptual issue of whether one should test the null
hypothesis using Bayesian model selection in the context of tests of GR. As
pointed out in Ref. [39], reporting the evidence against GR with Bayesian
estimators using free deviation parameters is questionable: results are prior-
dependent and not reparametrization-invariant, while credible intervals lack a
frequentist interpretation. On the other hand, a frequentist approach based on
the $p$-value only assesses the likelihood of the experimental outcome given
the null hypothesis and can be considered more resilient. Unfortunately,
implementing a pure $p$-value test in this context is, in practice, unfeasible
because one would need to know the true population distribution of the events.

Bootstrapping is a possible way out but only provides a partial solution.
Bootstrap samples inevitably inherit the peculiarities (e.g. outliers) of the
specific catalog realization we have observed. We speculate another promising
avenue in this direction is to incorporate population inference into tests of
GR [46] while relying on the notion of “Bayesian $p$-values” [47]. A safer
solution is to settle for weaker but more confident statements. This can be
done trivially by breaking down the catalogs into chunks, using fewer events
to compute the chosen estimator but obtaining multiple estimates.

While we concentrated on tests of GR, the cosmic variance is a generic effect,
both within and outside of GW astronomy (e.g. Ref. [48]). For instance,
astrophysical inferences from GW observations of binary populations [49] and
cosmological models [50] are impacted by the cosmic variance in much the same
fashion. The considerations put forward in this work are relevant to assess
the statistical significance of some of those findings, especially when the
significance itself is deemed to be weak.

_Acknowledgments_ — We thank Andrea Maselli, Riccardo Buscicchio, and Golam
Shaifullah for discussions. C.P. and D.G. are supported by ERC Starting Grant
No. 945155–GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant
No. 2022-Z9X4XS, and the ICSC National Research Centre funded by
NextGenerationEU. D.G. is supported by MSCA Fellowship No.
101064542–StochRewind and Leverhulme Trust Grant No. RPG-2019-350. S.B. is
supported by UKRI Stephen Hawking Fellowship No. EP/W005727. Computational
work was performed at CINECA with allocations through INFN and Bicocca.

## References

  * Baker  _et al._ [2015] T. Baker, D. Psaltis, and C. Skordis, Astrophys. J. 802, 63 (2015), arXiv:1412.3455 [astro-ph.CO] .
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ , Phys. Rev. Lett. 116, 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Berti  _et al._ [2015] E. Berti _et al._ , Class. Quant. Grav. 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Will [2014] C. M. Will, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Witek  _et al._ [2019] H. Witek, L. Gualtieri, P. Pani, and T. P. Sotiriou, Phys. Rev. D 99, 064035 (2019), arXiv:1810.05177 [gr-qc] .
  * Okounkova  _et al._ [2020] M. Okounkova, L. C. Stein, J. Moxon, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 101, 104016 (2020), arXiv:1911.02588 [gr-qc] .
  * Okounkova [2020] M. Okounkova, Phys. Rev. D 102, 084046 (2020), arXiv:2001.03571 [gr-qc] .
  * East and Ripley [2021] W. E. East and J. L. Ripley, Phys. Rev. D 103, 044040 (2021), arXiv:2011.03547 [gr-qc] .
  * Yunes  _et al._ [2016] N. Yunes, K. Yagi, and F. Pretorius, Phys. Rev. D 94, 084002 (2016), arXiv:1603.08955 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ , Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ , Phys. Rev. D 100, 104036 (2019), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ , (2021b), arXiv:2112.06861 [gr-qc] .
  * Del Pozzo  _et al._ [2011] W. Del Pozzo, J. Veitch, and A. Vecchio, Phys. Rev. D 83, 082002 (2011), arXiv:1101.1391 [gr-qc] .
  * Ghosh  _et al._ [2016] A. Ghosh _et al._ , Phys. Rev. D 94, 021101 (2016), arXiv:1602.02453 [gr-qc] .
  * Gossan  _et al._ [2012] S. Gossan, J. Veitch, and B. S. Sathyaprakash, Phys. Rev. D 85, 124056 (2012), arXiv:1111.5819 [gr-qc] .
  * Agathos  _et al._ [2014] M. Agathos, W. Del Pozzo, T. G. F. Li, C. Van Den Broeck, J. Veitch, and S. Vitale, Phys. Rev. D 89, 082001 (2014), arXiv:1311.0420 [gr-qc] .
  * Meidam  _et al._ [2014] J. Meidam, M. Agathos, C. Van Den Broeck, J. Veitch, and B. S. Sathyaprakash, Phys. Rev. D 90, 064009 (2014), arXiv:1406.3201 [gr-qc] .
  * Zimmerman  _et al._ [2019] A. Zimmerman, C.-J. Haster, and K. Chatziioannou, Phys. Rev. D 99, 124044 (2019), arXiv:1903.11008 [astro-ph.IM] .
  * Isi  _et al._ [2019] M. Isi, K. Chatziioannou, and W. M. Farr, Phys. Rev. Lett. 123, 121101 (2019), arXiv:1904.08011 [gr-qc] .
  * Isi  _et al._ [2022] M. Isi, W. M. Farr, and K. Chatziioannou, Phys. Rev. D 106, 024048 (2022), arXiv:2204.10742 [gr-qc] .
  * Mandel  _et al._ [2019] I. Mandel, W. M. Farr, and J. R. Gair, Mon. Not. Roy. Astron. Soc. 486, 1086 (2019), arXiv:1809.02063 [physics.data-an] .
  * Vitale  _et al._ [2020] S. Vitale, D. Gerosa, W. M. Farr, and S. R. Taylor, in  _Handbook of Gravitational Wave Astronomy_ (Springer, 2020) arXiv:2007.05579 [astro-ph.IM] .
  * Barausse  _et al._ [2014] E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D 89, 104059 (2014), arXiv:1404.7149 [gr-qc] .
  * Berti  _et al._ [2022] E. Berti, V. Cardoso, M. H.-Y. Cheung, F. Di Filippo, F. Duque, P. Martens, and S. Mukohyama, Phys. Rev. D 106, 084011 (2022), arXiv:2205.08547 [gr-qc] .
  * Saini  _et al._ [2022] P. Saini, M. Favata, and K. G. Arun, Phys. Rev. D 106, 084031 (2022), arXiv:2203.04634 [gr-qc] .
  * Bhat  _et al._ [2023] S. A. Bhat, P. Saini, M. Favata, and K. G. Arun, Phys. Rev. D 107, 024009 (2023), arXiv:2207.13761 [gr-qc] .
  * Moore  _et al._ [2021] C. J. Moore, E. Finch, R. Buscicchio, and D. Gerosa, iScience 24, 102577 (2021), arXiv:2103.16486 [gr-qc] .
  * Toubiana  _et al._ [2023] A. Toubiana, L. Pompili, A. Buonanno, J. R. Gair, and M. L. Katz, (2023), arXiv:2307.15086 [gr-qc] .
  * Vallisneri and Yunes [2013] M. Vallisneri and N. Yunes, Phys. Rev. D 87, 102002 (2013), arXiv:1301.2627 [gr-qc] .
  * Hu and Veitch [2023] Q. Hu and J. Veitch, Astrophys. J. 945, 103 (2023), arXiv:2210.04769 [gr-qc] .
  * Maggio  _et al._ [2023] E. Maggio, H. O. Silva, A. Buonanno, and A. Ghosh, Phys. Rev. D 108, 024043 (2023), arXiv:2212.09655 [gr-qc] .
  * Kamionkowski and Loeb [1997] M. Kamionkowski and A. Loeb, Phys. Rev. D 56, 4511 (1997), arXiv:astro-ph/9703118 .
  * Portsmouth [2004] J. Portsmouth, Phys. Rev. D 70, 063504 (2004), arXiv:astro-ph/0402173 .
  * Ivezić  _et al._ [2020] Ž. Ivezić, A. J. Connolly, J. T. VanderPlas, and A. Gray,  _Statistics, Data Mining, and Machine Learning in Astronomy_ (Princeton, 2020).
  * [35] N. J. Sloane _et al._ , The On-Line Encyclopedia of Integer Sequences Sequence A001700, https://oeis.org/A001700.
  * Ghosh  _et al._ [2021] A. Ghosh, R. Brito, and A. Buonanno, Phys. Rev. D 103, 124041 (2021), arXiv:2104.01906 [gr-qc] .
  * Dickey [1971] J. M. Dickey, The Annals of Mathematical Statistics 42, 204 (1971).
  * Jeffreys [1998] H. Jeffreys,  _The Theory of Probability_ (Oxford, 1998).
  * Chua and Vallisneri [2020] A. J. K. Chua and M. Vallisneri, (2020), arXiv:2006.08918 [gr-qc] .
  * Vallisneri [2008] M. Vallisneri, Phys. Rev. D 77, 042001 (2008), arXiv:gr-qc/0703086 .
  * Schutz [2011] B. F. Schutz, Class. Quant. Grav. 28, 125023 (2011), arXiv:1102.5421 [astro-ph.IM] .
  * Speagle [2020] J. S. Speagle, Mon. Not. Roy. Astron. Soc. 493, 3132 (2020), arXiv:1904.02180 [astro-ph.IM] .
  * Brito  _et al._ [2018] R. Brito, A. Buonanno, and V. Raymond, Phys. Rev. D 98, 084038 (2018), arXiv:1805.00293 [gr-qc] .
  * Dreyer  _et al._ [2004] O. Dreyer, B. J. Kelly, B. Krishnan, L. S. Finn, D. Garrison, and R. Lopez-Aleman, Class. Quant. Grav. 21, 787 (2004), arXiv:gr-qc/0309007 .
  * Berti  _et al._ [2006] E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D 73, 064030 (2006), arXiv:gr-qc/0512160 .
  * Payne  _et al._ [2023] E. Payne, M. Isi, K. Chatziioannou, and W. M. Farr, (2023), arXiv:2309.04528 [gr-qc] .
  * Vallisneri  _et al._ [2023] M. Vallisneri, P. M. Meyers, K. Chatziioannou, and A. J. K. Chua, (2023), arXiv:2306.05558 [astro-ph.HE] .
  * Roebber  _et al._ [2016] E. Roebber, G. Holder, D. E. Holz, and M. Warren, Astrophys. J. 819, 163 (2016), arXiv:1508.07336 [astro-ph.CO] .
  * Abbott  _et al._ [2023a] R. Abbott _et al._ , Phys. Rev. X 13, 011048 (2023a), arXiv:2111.03634 [astro-ph.HE] .
  * Abbott  _et al._ [2023b] R. Abbott _et al._ , Astrophys. J. 949, 76 (2023b), arXiv:2111.03604 [astro-ph.CO] .


executed, returning True. In the former case the final $W$ has $|\varphi|+1$
worlds, and the structure represents a looping model: $M$ is built on the
frame $(W,R)$ where $R=\\{(w_{k^{\prime}},w_{k^{\prime\prime}})\mid
k^{\prime}\leq k^{\prime\prime}$ or $k^{\prime}=k^{\prime\prime}=|\varphi|\\}$
(notice that $R$ self-loops on $w_{|\varphi|}$). In the latter case (Shorten
is executed) we have that $|W|\leq|\varphi|$, and the structure represents a
finite model with no loops: $M$ is built on the frame $(W,R)$ where
$R=\\{(w_{i},w_{j})\mid i\leq j\\}$. In both cases, for every $w\in W$ we
define the valuation as $V(w)=L(w)|_{\mathcal{P}}$. We want to prove that
$M,w_{0}\Vdash\varphi$, where $w_{0}$ is the root of the pre-linear
$\operatorname{\mathbf{K4}}$-model $M$. We prove the correctness of the
procedure only for the case of a ‘looping’ pre-linear model. The case of a
non-looping model is simpler and can be proved in the very same way. Suppose,
by way of contradiction, that $M,w_{0}\not\Vdash\varphi_{i}$ for some clause
$\varphi_{i}$. Assume that
$\varphi_{i}={\Box}^{s}(\lambda_{1}\wedge\ldots\wedge\lambda_{n}\rightarrow\lambda)$,
and let $w_{0},w_{1},\ldots$ be the enumeration of the of worlds in $W$.
Hence, we have that there exists $w_{t}\in W$, for some $t\geq s$, such that
$M,w_{t}\Vdash\lambda_{1}\wedge\ldots\wedge\lambda_{n}$ but
$M,w_{t}\not\Vdash\lambda$. Since $w_{t}\in W$, the for cycle of the main
procedure has run enough times to create $w_{t}$, and since
$\operatorname{\mathbf{K4}}$-Saturate has been executed after that $w_{t}$ has
been created, $(\lambda_{1}\wedge\ldots\wedge\lambda_{n}\rightarrow\lambda)$
has been inserted in $H(w_{t})$ (lines 12-16). We first show that, for every
positive literal $\lambda_{j}$ in the head of the clause, we have
$\lambda_{j}\in H(w_{t})\cup L(w_{t})$. If $\lambda_{j}=p$, then, since
$M,w_{t}\Vdash p$ we have that $p\in L(w_{t})$ by the definition of $M$. If
$\lambda_{j}={\Box}^{m}p$, since we know that $M,w_{t}\Vdash\lambda_{j}$, it
must be the case that either $t+m\leq|\varphi|$ and $M,w_{u}\Vdash p$ for each
$u\geq t+m$, or $t+m>|\varphi|$ and $M,w_{|\varphi|}\Vdash p$ (because of the
self-loop in the model). Suppose that we are in the latter case, and let
$N=|\varphi|<t+m$: since $M,w_{N}\Vdash p$, by the definition of $M$ we have
that this implies that $p\in L(w_{N})$. Hence, the last execution of
$\operatorname{\mathbf{K4}}$-saturate must have executed lines 4–5, adding
${\Box}p$ to $H(w_{N})$. But, from lines 4–5, we have that also
${\Box}{\Box}p,\ldots,{\Box}^{m}p$ are in $H(w_{N})$, while from lines 6–8 we
can conclude that ${\Box}^{m}p\in H(w_{t})$. When $t+m\leq|\varphi|$, a
similar argument allow us to conclude that also in this case $\lambda_{j}\in
H(w_{t})$. Since every positive literal in the body of the clause
$(\lambda_{1}\wedge\ldots\wedge\lambda_{n}\rightarrow\lambda)$ belongs to
$H(w_{t})$, we have that some execution of
$\operatorname{\mathbf{K4}}$-Saturate must have put $\lambda\in H(w_{t})$
(lines 17-20). If $\lambda=p$, within the same execution of
$\operatorname{\mathbf{K4}}$-Saturate we would have that $p\in L(w_{t})$
(lines 9-11), and thus that $M,w_{t}\Vdash p$, which contradicts the
hypothesis that $M$ does not satisfy the clause. Conversely, if
$\lambda={\Box}^{m}p$, then, the execution of
$\operatorname{\mathbf{K4}}$-Saturate that takes place either at the moment of
creation of the point $w_{t+m}$ (when $t+m\leq N$), or at the creation of the
point $w_{N}$ (when $t+m>N$) puts $p$ respectively in $L(w_{t+m})$ or
$L(w_{N})$, which, again, contradicts our hypothesis. In all cases we have a
contradiction and thus we have proved that $M,w_{0}\Vdash\varphi$.

Finally, the case of $\operatorname{\mathbf{L}}=\operatorname{\mathbf{S4}}$
can be proved by combining the proofs for $\operatorname{\mathbf{T}}$ and
$\operatorname{\mathbf{K4}}$.

The definition of compatibility of a pre-linear non-looping
$\operatorname{\mathbf{K}}$-model $M$ with respect to a structure $(W,H,L)$
given above applies also to pre-linear reflexive, transitive, and reflexive
and transitive models. The following lemma shows that the saturation
procedures preserve compatibility also for $\operatorname{\mathbf{T}}$,
$\operatorname{\mathbf{K4}}$, and $\operatorname{\mathbf{S4}}$.

###### Lemma 21 .

Let
$\operatorname{\mathbf{L}}\in\\{\operatorname{\mathbf{T}},\operatorname{\mathbf{K4}},\operatorname{\mathbf{S4}}\\}$,
let $M$ be a pre-linear $\operatorname{\mathbf{L}}$-model, and let $(W,H,L)$
be a structure compatible with $M$. Then the structure obtained after the
execution of $\operatorname{\mathbf{L}}$-Saturate($W,H,L$) is compatible with
$M$.

###### Proof 4.9 .

We first show that the claim holds for
$\operatorname{\mathbf{L}}=\operatorname{\mathbf{T}}$. We proceed as in the
proof of Lemma 18 to show that after every iteration of the while loop the
structure $(W,H,L)$ remains compatible with $M$. From now on, we assume that
$M$ has no loops and we identify the worlds of the structure $(W,H,L)$ as
$w_{0},w_{1},\ldots$, and the worlds in $M$ as
$w_{0}^{\prime},w_{1}^{\prime},\ldots$. Let $w_{k}\in W$ and $\psi\in
H(w_{k})$ be respectively the world and the formula selected by the let
statement of line 3. If the corresponding $w_{i}^{\prime}$ is not a world of
$M$, then the property trivially holds. Otherwise, we proceed by following the
pseudocode. Lines 4–6 insert ${\Box}\psi\in H(w_{k-1})$ if the predecessor
$w_{k-1}$ of $w_{k}$ exists, ${\Box}\psi\in Cl(\varphi)$ and $\psi\in
H(w_{k-1})\cup L(w_{k-1})$. By inductive hypothesis, since $\psi$ was already
in $H(w_{k})\cup L(w_{k})$ and in $H(w_{k-1})\cup L(w_{k-1})$, it was already
true that $M,w_{k}^{\prime}\Vdash\psi$ and $M,w_{k-1}^{\prime}\Vdash\psi$.
Since $M$ is a reflexive pre-linear model, we also have that
$M,w_{k-1}^{\prime}\Vdash{\Box}\psi$, as required by the fact that now
${\Box}\psi\in H(w_{k-1})\cup L(w_{k-1})$. The other steps depend on the
structure of $\psi$. If $\psi=p$, then lines 7–9 move $\psi$ from $H(w_{k})$
to $L(w_{k})$. By inductive hypothesis, since $p$ was already in $H(w_{k})\cup
L(w_{k})$, it was already true that $M,w_{k}^{\prime}\Vdash p$. If
$\psi={\Box}\xi$, then line 11 inserts $\xi$ in $H(w_{k})$. Moreover, if there
exists a successor $w_{k+1}$ of $w_{k}$, then lines 12–15 move ${\Box}\xi$
from $H(w_{k})$ to $L(w_{k})$, and insert $\xi\in H(w_{k+1})$. By inductive
hypothesis, $M,w_{k}^{\prime}\Vdash{\Box}\xi$: by reflexivity, this means that
$M,w_{k}^{\prime}\Vdash\xi$ and that, if the successor $w_{k+1}^{\prime}$
exists in $M$, $M,w_{k+1}^{\prime}\Vdash\xi$, as required. The proof for the
cases when $\psi$ is a clause is similar to Lemma 18.

Let us prove now that the claim holds for
$\operatorname{\mathbf{L}}=\operatorname{\mathbf{K4}}$. We proceed as before
to show that after every iteration of the while loop the structure $(W,H,L)$
remains compatible with $M$. Again, let $w_{k}\in W$ and $\psi\in H(w_{k})$ be
respectively the world and the formula selected by the let statement of line
3. If the corresponding $w_{k}^{\prime}$ is not a world of $M$, then the
property trivially holds. Otherwise, we proceed by following the pseudocode by
considering lines 4–5 that inserts ${\Box}\psi\in H(w_{N})$ if ${\Box}\psi\in
Cl(\varphi)$ and $\psi\in H(w_{N})\cup L(w_{N})$. By inductive hypothesis,
since $\psi$ was already in $H(w_{N})\cup L(w_{N})$, it was already true that
$M,w_{N}^{\prime}\Vdash\psi$. Since $w_{N}$ is the last world of a looping
model we also have that $M,w_{N}^{\prime}\Vdash{\Box}\psi$, as required by the
fact that now ${\Box}\psi\in H(w_{N})\cup L(w_{N})$. Then, lines 6–8 insert
${\Box}\psi\in H(w_{k-1})$ if the predecessor $w_{k-1}$ of $w_{k}$ exists,
${\Box}\psi\in Cl(\varphi)$ and ${\Box}\psi\in H(w_{k})\cup L(w_{k})$. By
inductive hypothesis, since $\psi$ and ${\Box}\xi$ were already in
$H(w_{k})\cup L(w_{k})$, it was already true that $M,w_{k}^{\prime}\Vdash\psi$
and $M,w_{k}^{\prime}\Vdash{\Box}\psi$. Since $M$ is a transitive pre-linear
model, we also have that $M,w_{k-1}^{\prime}\Vdash{\Box}\psi$, as required by
the fact that now ${\Box}\psi\in H(w_{k-1})\cup L(w_{k-1})$. The other steps
depend on the structure of $\psi$. If $\psi=p$, then lines 9–11 move $\psi$
from $H(w_{k})$ to $L(w_{k})$. By inductive hypothesis, since $p$ was in
$H(w_{k})\cup L(w_{k})$, it was already true that $M,w_{k}^{\prime}\Vdash p$.
In the case $\psi={\Box}\xi$, if there exists a successor $w_{k+1}$ of
$w_{k}$, then lines 13–16 move ${\Box}\xi$ from $H(w_{k})$ to $L(w_{k})$, and
insert $\xi,{\Box}\xi\in H(w_{k+1})$. By inductive hypothesis,
$M,w_{k}^{\prime}\Vdash{\Box}\xi$, which imply that
$M,w_{k+1}^{\prime}\Vdash\xi$ and $M_{k+1}^{\prime}\Vdash{\Box}\xi$, as we
wanted. The rest of the cases are treated as in Lemma 18.

Once we have proved that $\operatorname{\mathbf{L}}$-Saturate preserves
compatibility, proving that HornBoxSat is complete becomes straightforward.

###### Lemma 22 .

Let
$\operatorname{\mathbf{L}}\in\\{\operatorname{\mathbf{T}},\operatorname{\mathbf{K4}},\operatorname{\mathbf{S4}}\\}$
and let $\varphi$ be a $\operatorname{\mathbf{L^{Horn,\Box}}}$-formula. If
$\varphi$ is satisfiable then HornBoxSat($\operatorname{\mathbf{L}},\varphi$)
returns True.

###### Proof 4.10 .

We can proceed as in the proof of Lemma 19. We assume that $\varphi$ is
satisfiable, and thanks to Theorem 15, we know that there exist some pre-
linear models of size bounded either by $md(\varphi)+1$ (for
$\operatorname{\mathbf{T}}$) or by $|\varphi|$ (for
$\operatorname{\mathbf{K4}}$ and $\operatorname{\mathbf{S4}}$). Since our
procedure starts iterating from a structure $(W,H,L)$ compatible with $M$, and
since $\operatorname{\mathbf{L}}$-Saturate preserves compatibility (Lemma 21),
we can prove, as in Lemma 19, that
HornBoxSat($\operatorname{\mathbf{L}},\varphi$) builds the longest model $M$
of $\varphi$, and therefore returns True.

The following theorem summarises the above results and shows that the
procedure is correct, complete and polynomial for all the considered
languages.

###### Theorem 23 .

Let
$\operatorname{\mathbf{L}}\in\\{\operatorname{\mathbf{K}},\operatorname{\mathbf{T}},\operatorname{\mathbf{K4}},\operatorname{\mathbf{S4}}\\}$
and let $\varphi$ be a $\operatorname{\mathbf{L^{Horn,\Box}}}$-formula. Then,
the procedure HornBoxSat($\mathbf{L},\varphi$) is correct, complete and
terminates within $O(poly(|\varphi|))$ time.

###### Proof 4.11 .

Correctness and completeness of HornBoxSat follows from Lemma 17 and 19 when
$\operatorname{\mathbf{L}}=\operatorname{\mathbf{K}}$ and from Lemma 20 and 22
when
$\operatorname{\mathbf{L}}\in\\{\operatorname{\mathbf{T}},\operatorname{\mathbf{K4}},\operatorname{\mathbf{S4}}\\}$.

As for the time complexity of the procedure, we observe that the most external
cycle runs $O(|\varphi|)$ times in the worst case. The dimension of the
candidate models grows from 1 to $O(|\varphi|)$ points, and, then, down to 1
again; in the $H$ component of each point there are, at most, $O(|\varphi|)$
formulas on which Saturate has effect. Therefore, the total time spent is
$O(|\varphi|)\cdot 2\cdot\Sigma_{i=1}^{O(|\varphi|)}i=O(|\varphi|^{3})$, which
is polynomial in $|\varphi|$.

We conclude this section by studying the particular case of the fragment
$\operatorname{\mathbf{K^{core,\Box}}}$. The results obtained so far allow us
only to fix the complexity of the satisfiability problem in these cases
between NLogSpace and P, the lower bound being a consequence of the NLogSpace-
completeness of the satisfiability problem for binary propositional clauses
(i.e., the 2SAT problem [Pap03]). We want to prove that the satisfiability
problem for $\operatorname{\mathbf{K^{core,\Box}}}$, the satisfiability
problem is, in fact, NLogSpace-complete; to this end, we extend the algorithm
for 2SAT. The standard algorithm to solve a 2SAT instance exploits the
properties of the _implication graph_ of the formula. An implication graph is
a directed graph in which there is one vertex per propositional variable or
negated variable, and an edge connecting one vertex to another whenever the
corresponding variables are related by an implication in the formula. A
propositional binary formula is unsatisfiable if and only if there exists a
path in the implication graph that starts from a variable $p$, leads to its
negation $\neg p$ and then goes back to $p$. The NLogSpace algorithm simply
guesses the variable $p$ and search nondeterministically for the chain of
implications that corresponds to such a contradictory path. An alternative,
indirect way to obtain this result would be to explore the characteristics of
SLD-resolution [Ngu03] under the Horn box restriction of
$\operatorname{\mathbf{K}}$.

In the implication graph of a $\operatorname{\mathbf{K^{core,\Box}}}$-formula
the vertexes are pairs where the first component is a positive literal
$\lambda$ or its negation $\neg\lambda$, and the second component is a natural
number representing the depth of the world where $\lambda$
(resp.,$\neg\lambda$) holds, measured as the distance from the root of the
model; it is worth to recall that $\operatorname{\mathbf{K^{core,\Box}}}$, as
$\operatorname{\mathbf{K^{Horn,\Box}}}$, enjoys the pre-linear model property,
so that the notion of distance from the root makes perfect sense. The edges in
the graph correspond either to implications in the formula (as in the
propositional case), or to jumps from a world at some depth $d$ to the world
at depth $d+1$ or $d-1$. Thus, given
$\varphi\in\operatorname{\mathbf{K^{Horn,\Box}}}$, where
$\varphi_{1},\ldots,\varphi_{l}$, and given a natural number $D\leq
md(\varphi)$, we define the _$D$ -implication graph of $\varphi$_ as a
directed graph $G_{\varphi}^{D}=(V,E)$ where $V=\\{(\lambda,d)\mid\lambda\in
Cl(\varphi)$ or $\neg\lambda\in Cl(\varphi)$ and $0\leq d<D\\}$, and:  _(i)_
$((\lambda,d),(\lambda^{\prime},d)),((\neg\lambda^{\prime},d),(\neg\lambda,d))\in
E$if and only if $\varphi_{i}={\Box}^{d}(\neg\lambda\lor\lambda^{\prime})$ for
some $i$;  _(ii)_ $(({\Box}\lambda,d),(\lambda,d+1))\in E$if and only if
$d<D-1$;  _(iii)_ $((\lambda,d),({\Box}\lambda,d-1))\in E$if and only if
$d>0$;  _(iv)_ $((\neg{\Box}\lambda,d),(\neg\lambda,d+1))\in E$if and only if
$d<D-1$; and  _(v)_ $((\neg\lambda,d),(\neg{\Box}\lambda,d-1))\in E$if and
only if $d>0$.

Notice that, in this case, we treat
$\operatorname{\mathbf{K^{core,\Box}}}$-clauses as disjunctions, instead of
implications; in this way, clauses that contain a single positive (resp.,
negative) literal $\lambda$ (resp., $\neg\lambda$) are dealt with as
disjunctions of the type $(\lambda\vee\bot)$ (resp., $(\neg\lambda\vee\bot)$).
Moreover, for technical reasons, we add to each formula the clauses
${\Box}^{s}(\bot\rightarrow\top)$ for each $s\leq md(\varphi)$. The main
property of the $D$-implication graph is that the formula on which built is
unsatisfiable on (pre-linear) models of length $D$ if and only if there exists
a cycle $\pi$ that starts and ends in some vertex $(\lambda,d)$ and visits the
vertex $(\neg\lambda,d)$; in this case, $\pi$ is called contradictory cycle.

###### Lemma 24 .

Let $\varphi$ be a formula of $\operatorname{\mathbf{K^{Horn,\Box}}}$. Then,
if there exists a contradictory cycle $\pi$ in $G_{\varphi}^{D}$, then
$\varphi$ is unsatisfiable on a linear model with precisely $D$ distinct
worlds.

###### Proof 4.12 .

Let $\varphi$ be a formula of $\operatorname{\mathbf{K^{Horn,\Box}}}$, and
suppose that there exists a vertex $(\lambda,d)$ and a contradictory cycle
$\pi=(\lambda_{1},d_{1}),(\lambda_{2},d_{2}),\ldots,(\lambda_{|\pi|},d_{|\pi|})$
in the implication graph $G_{\varphi}^{D}$ such that
$(\lambda,d)=(\lambda_{1},d_{1})$. Suppose by contradiction that $\varphi$ is
satisfiable on a pre-linear model $M$ with precisely $D$ distinct worlds such
that $M,w_{0}\Vdash\varphi$. Two cases may arise: either
$M,w_{d}\Vdash\lambda$ or $M,w_{d}\Vdash\neg\lambda$. We consider here only
the former case (the latter one can be proved in the very same way), and we
show that for every vertex $(\lambda_{i},d_{i})\in\pi$ it is the case that
$M,w_{d_{i}}\Vdash\lambda_{i}$. We proceed by induction on $i$. When $i=1$
then $(\lambda_{1},d_{1})=(\lambda,d)$, and since $M,w_{d}\Vdash\lambda$ the
property trivially holds. Now, suppose that the property is true for the
vertex $(\lambda_{i},d_{i})$, and consider the $i+1$-th vertex
$(\lambda_{i+1},d_{i+1})$; three cases may arise:

  * •

$d_{i+1}=d_{i}$: in this case we have that the clause
${\Box}^{d_{i}}(\neg\lambda_{i}\lor\lambda_{i+1})$ is a conjunct of $\varphi$.
Since $M,w_{0}\Vdash{\Box}^{d_{i}}(\neg\lambda_{i}\lor\lambda_{i+1})$ and
$w_{i}$ is at distance $d_{i}$ form $w_{0}$, we have that
$M,w_{d_{i}}\Vdash\neg\lambda_{i}\lor\lambda_{i+1}$. Since by inductive
hypothesis we have that $M,w_{d_{i}}\Vdash\lambda_{i}$, then, as we wanted,
$M,w_{d_{i+1}}\Vdash\lambda_{i+1}$.

  * •

$d_{i+1}=d_{i}+1$: in this case we have that either
$\lambda_{i}={\Box}\lambda$ or $\lambda_{i}=\neg{\Box}\lambda$. In the first
case, the inductive hypothesis is that $M,w_{d_{i}}\Vdash{\Box}\lambda$, which
is to say that $M,w_{d_{i}+1}\Vdash\lambda$, that is,
$M,w_{d_{i}+1}\Vdash\lambda_{i+1}$, as we wanted. In the second case the
inductive hypothesis is that $M,w_{d_{i}}\Vdash\neg{\Box}\lambda$, which is to
say that $M,w_{d_{i}+1}\Vdash\neg\lambda$, that is,
$M,w_{d_{i}+1}\Vdash\lambda_{i+1}$, as we wanted.

  * •

$d_{i+1}=d_{i}-1$: in this case we have that either
$\lambda_{i+1}={\Box}\lambda$ or $\lambda_{i+1}=\neg{\Box}\lambda$. In the
first case, it must be that $\lambda_{i}=\lambda$; therefore, by inductive
hypothesis, $M,w_{d_{i}}\Vdash\lambda$, which is to say that
$M,w_{d_{i}-1}\Vdash{\Box}\lambda$ (since $M$ is a pre-linear model), that is,
$M,w_{d_{i}-1}\Vdash\lambda_{i+1}$, as we wanted. In the second case it must
be that $\lambda_{i}=\neg\lambda$; therefore, by inductive hypothesis,
$M,w_{d_{i}}\Vdash\neg\lambda$, which is to say that
$M,w_{d_{i}-1}\Vdash\neg{\Box}\lambda$, that is,
$M,w_{d_{i}-1}\Vdash\lambda_{i+1}$, as we wanted.

Since both $(\lambda,d)$ and $(\neg\lambda,d)$ belong to $\pi$, by the above
property we have that $M,w_{d}\Vdash\lambda$ and $M,w_{d}\Vdash\neg\lambda$,
and a contradiction is found. Hence, $\varphi$ cannot be satisfied in a pre-
linear model with precisely $D$ distinct worlds.

###### Lemma 25 .

Let $\varphi$ be a formula of $\operatorname{\mathbf{K^{Horn,\Box}}}$. Then,
if there is no contradictory cycle $\pi$ in $G_{\varphi}^{D}$, then $\varphi$
is satisfiable on a linear model with precisely $D$ distinct worlds.

###### Proof 4.13 .

Let $\varphi$ be a formula of $\operatorname{\mathbf{K^{Horn,\Box}}}$, and
suppose that no contradictory cycle is in $G_{\varphi}^{D}$. We show how to
build a pre-linear model for $\varphi$ with exactly $D$ distinct worlds. Let
$W=\\{w_{0},\ldots,w_{D-1}\\}$, and let $(w_{i},w_{i+1})\in R$ for each $i\leq
D-1$. Now, we have to show how to build a consistent evaluation for each
world. As a preliminary step, we show the following property on the
$D$-implication graph for $\varphi$: if
$(\lambda,d)\rightsquigarrow(\lambda^{\prime},d^{\prime})$ (where
$\rightsquigarrow$ means that there is a path from $(\lambda,d)$ to
$(\lambda^{\prime},d^{\prime})$), then
$(\neg\lambda^{\prime},d^{\prime})\rightsquigarrow(\neg\lambda,d)$. We can
prove this by induction on the length of the path. The base case is $|\pi|=2$
(the case in which $|\pi|=1$ is trivial), where three cases may arise:

  * •

$d=d^{\prime}$. In this case, we have that the clause
${\Box}^{d}(\neg\lambda\lor\lambda^{\prime})$ is a conjunct of $\varphi$; but
the rules for the construction of the implication graph include the edge
$(\neg\lambda^{\prime},d),(\neg\lambda,d)$ as well, and, thus, we have the
claim.

  * •

$d^{\prime}=d+1$. In this case we have that either
$\lambda={\Box}\lambda^{\prime}$, or $\lambda=\neg{\Box}\overline{\lambda}$
and $\lambda^{\prime}=\neg\overline{\lambda}$. In the first case, the rules
for the construction of the implication graph include the edge
$(\neg\lambda^{\prime},d+1),(\neg{\Box}\lambda^{\prime},d)$, implying the
claim; in the second case the rules for the construction of the implication
graph include the edge
$(\overline{\lambda},d+1),({\Box}\overline{\lambda},d)$, implying, again, the
claim.

  * •

$d^{\prime}=d-1$. This case is completely symmetric to the above one.

The inductive case can be proved by using precisely the same rules; the
property holds up to a certain length $i$, and, then, a single application of
the above cases leads to prove that it holds for the length $i+1$. Now, we can
build a function $V^{*}:W\rightarrow 2^{Cl(\varphi)}$, as follows. Starting
from the world $w_{0}$, we choose $\lambda\in Cl(\varphi)$, and we set
$\lambda\in V(w_{0})$ if and only if there is no path from $(\lambda,0)$ to
$(\neg\lambda,0)$ in $G_{\varphi}^{D}$. Then, for each $(\lambda^{\prime},d)$
reachable from $(\lambda,0)$ in $G_{\varphi}^{D}$, we set $\lambda^{\prime}\in
V(w_{d})$. We now repeat the same process for each $\lambda^{\prime}\in
Cl(\varphi)$ for which a decision for $w_{0}$ has not been taken yet. Finally,
we repeat all three steps for each $w_{d}\in W$. At this point, we can
synthesize an evaluation $V:W\rightarrow 2^{\mathcal{P}}$ from $V^{*}$ by
simply projecting $V^{*}$ over the propositional letters. Since the above
property implies that such a labeling is contradiction-free, it is easy to see
that $M,w_{0}\Vdash\varphi$.

The above lemmas can be then immediately used to devise a non-deterministic
algorithm for $\operatorname{\mathbf{K^{core,\Box}}}$: we simply guess the
length of the candidate model $D\leq md(\varphi)$, and we check that no
contradictory cycle exists, which proves that $\varphi$ has a pre-linear model
of dimension $D$. Then, the following result holds.

###### Theorem 26 .

The satisfiability problem for $\operatorname{\mathbf{K^{core,\Box}}}$ is
NLogSpace-complete.

###### Proof 4.14 .

The fact that the satisfiability problem for
$\operatorname{\mathbf{L^{core,\Box}}}$ is in NLogSpace is a consequence of
the above argument and the fact that checking the non-existence of a cycle in
a direct graph is NLogSpace. To avoid the explicit construction of the entire
graph, that would require a polynomial amount of space, the NLogSpace
algorithm simulates the construction and checks for the non-existence of
contradictory cycles on-the-fly. In [Pap03] NLogSpace-hardness of 2SAT is
proved by reducing the unreachability problem on acyclic graphs to 2SAT. Since
the formula used for the reduction includes only clauses of the form $(\neg
p\lor q)$, $p$, and $\neg p$, it is indeed a formula of the core fragment of
propositional logic. Hence, NLogSpace-hardness of
$\operatorname{\mathbf{L^{core,\Box}}}$ immediately follows.

The complexity results of this of this section can be seen in Fig. 4.

$\operatorname{\mathbf{L^{core}}}$$\operatorname{\mathbf{L^{Bool}}}$$\operatorname{\mathbf{L^{core,\Box}}}$$\operatorname{\mathbf{L^{Krom}}}$$\operatorname{\mathbf{L^{Horn,\Box}}}$$\operatorname{\mathbf{L^{Horn}}}$$\operatorname{\mathbf{L^{Krom,\Diamond}}}$$\operatorname{\mathbf{L^{Horn,\Diamond}}}$$\operatorname{\mathbf{L^{Krom,\Box}}}$$\operatorname{\mathbf{L^{core,\Diamond}}}$$\in\textsc{P}$
[Ngu04], $\mathbf{L}=\operatorname{\mathbf{S5}}$ PSpace-complete [CL94],
---
$\mathbf{L}=\operatorname{\mathbf{K}},\operatorname{\mathbf{T}},\operatorname{\mathbf{K4}},\operatorname{\mathbf{S4}}$
NP-complete [Ngu04], $\mathbf{L}=\operatorname{\mathbf{S5}}$ NLogSpace-hard,
$\in\textsc{P}$,
---
$\mathbf{L}=\operatorname{\mathbf{T}},\operatorname{\mathbf{K4}},\operatorname{\mathbf{S4}}$,
NLogSpace-complete, $\mathbf{L}=\operatorname{\mathbf{K}}$
P-complete,
---
$\mathbf{L}=\operatorname{\mathbf{K}},\operatorname{\mathbf{T}},\operatorname{\mathbf{K4}},\operatorname{\mathbf{S4}}$
Figure 4. Complexity results

##  5\. Conclusions

In this paper, we considered the expressive power and complexity of sub-
propositional fragments of the basic modal logics of all directed graphs
($\operatorname{\mathbf{K}}$) and its axiomatic reflexive
($\operatorname{\mathbf{T}}$), transitive ($\operatorname{\mathbf{K4}}$), and
reflexive and transitive ($\operatorname{\mathbf{S4}}$) extensions. The
relative expressive power and the complexity results of the sub-propositional
fragments of modal logic studied in this paper is depicted in Figures 3 and 4,
respectively. In most cases relative expressivity coincides with syntactical
containment, with the notable exception of the Krom fragments, that are
expressively equivalent, but not weakly expressively equivalent. Because of
our very general approach for comparing the expressive power of languages,
most of our result can be transferred to other sub-propositional modal logics
such as the fragments of $\mathbf{LTL}$ without Since and Until studied in
[AKRZ13] and the sub-propositional fragments of $\mathbf{HS}$ [AKRZ15,
BKMV+17, BMS14]. To the best of our knowledge, this is the first work where
sub-Krom and sub-Horn fragments of these modal logics have been considered.
Concerning complexity, starting from the known result that every modal logic
between $\operatorname{\mathbf{K}}$ and $\operatorname{\mathbf{S4}}$ is still
PSpace-complete even under the Horn restriction [Ngu04] (in opposition to the
case of $\operatorname{\mathbf{S5}}$, which goes from NP-complete to
P-complete), we proved that eliminating the use of diamonds in Horn formulas
reduces the complexity of the satisfiability problem in all considered cases;
some of our results can be also proved by using the machinery in [Ngu00,
Ngu16], which, on the other hand, does not immediately provide an
implementable satisfiability-checking procedure. We proved that the
satisfiability problem becomes P-complete in the Horn fragment without
diamonds, between NLogSpace and P in the core fragment without diamonds (and
NLogSpace-complete in the special case of
$\operatorname{\mathbf{K^{core,\Box}}}$), and we devised a fast
satisfiability-checking algorithms for all considered cases.

While this study has been inspired from similar results in the temporal case,
such as for fragments of $\mathbf{LTL}$ without Since and Until and fragment
of $\mathbf{HS}$, several natural questions concerning sub-Horn temporal
logics remain open, and our results can certainly offer a solid starting point
to their tackling. To mention one interesting case, consider $\mathbf{HS}$
under the Horn restriction without diamonds: it is known that it is still
undecidable in the discrete case [BKMV+17], but such a proof makes extensive
use of many different temporal operators; what would happen by considering
simple fragments such as $\mathrm{A}$ or $\mathrm{D}$ (i.e., the fragments
with a single modal operator corresponding to, respectively, Allen’s relation
meets and Allen’s relation during) is unclear. Our technique, based on the
intrinsic inability of $\operatorname{\mathbf{S4}}$ to force the non-linearity
of a model under the considered restriction, seems, at least in principle,
applicable to such a case. Other interesting cases include the behaviour of
$\mathbf{LTL}$ with Since/Until under the restrictions considered in this
paper.

## Acknowledgements

The authors acknowledge the support from the Italian INdAM – GNCS Project 2018
“ _Formal methods for verification and synthesis of discrete and hybrid
systems_ ” (D. Bresolin and G. Sciavicco), and the Spanish Project
TIN15-70266-C2-P-1 (E. Muñoz-Velasco).

## References

  * [AKRZ13] A. Artale, R. Kontchakov, V. Ryzhikov, and M. Zakharyaschev.  The complexity of clausal fragments of LTL.  In Proc. of the 19th International Conference on Logic for Programming, Artificial Intelligence, and Reasoning (LPAR), volume 8312 of LNCS, pages 35–52. Springer, 2013.
  * [AKRZ14] A. Artale, R. Kontchakov, V. Ryzhikov, and M. Zakharyaschev.  A cookbook for temporal conceptual data modelling with description logics.  ACM Transactions on Computational Logic, 15(3):1–50, 2014.
  * [AKRZ15] A. Artale, R. Kontchakov, V. Ryzhikov, and M. Zakharyaschev.  Tractable interval temporal propositional and description logics.  In Proc. of the 29th AAAI Conference on Artificial Intelligence (AAAI), pages 1417–1423, 2015.
  * [APT79] B. Aspvall, M. F. Plass, and R. E. Tarjan.  A linear time algorithm for testing the truth of certain quantified boolean formulas.  Information Processing Letters, 8(3):121–123, 1979.
  * [APvB07] M. Aiello, I. Pratt-Hartmann, and J. van Benthem, editors.  Handbook of Spatial Logics.  Springer, 2007.
  * [BCM+03] F. Baader, D. Calvanese, D.L. McGuinness, D. Nardi, and P.F. Patel-Schneider, editors.  The Description Logic Handbook: Theory, Implementation, and Applications.  Cambridge University Press, 2003.
  * [BFG+90] H. Barringer, M. Fisher, D. Gabbay, G. Gough, and R. Owens.  Metatem: A framework for programming in temporal logic.  In Proc. of the Workshop Stepwise Refinement of Distributed Systems Models, Formalisms, Correctness. Springer, 1990.
  * [BJP+12] R. Bloem, B. Jobstmann, N. Piterman, A. Pnueli, and Y. Sa’ar.  Synthesis of reactive(1) designs.  Journal of Computer and System Sciences, 78(3):911–938, 2012.
  * [BKMV+17] D. Bresolin, A. Kurucz, E. Muñoz-Velasco, V. Ryzhikov, G. Sciavicco, and M. Zakharyaschev.  Horn fragments of the Halpern-Shoham interval temporal logic.  ACM Transactions on Computational Logic, 18(3):22:1 – 22:39, 2017\.
  * [BMS14] D. Bresolin, E. Muñoz-Velasco, and G. Sciavicco.  Sub-propositional fragments of the interval temporal logic of Allen’s relations.  In Proc. of the 14th European Conference on Logics in Artificial Intelligence (JELIA), volume 8761 of LNCS, pages 122–136. Springer, 2014\.
  * [BvBW06] P. Blackburn, J.F. van Benthem, and F. Wolter.  Handbook of modal logic, volume 3.  Elsevier, 2006.
  * [CL93] C.C. Chen and I.P. Lin.  The computational complexity of satisfiability of temporal Horn formulas in propositional linear-time temporal logic.  Information Processing Letters, 45(3):131–136, 1993.
  * [CL94] C.C. Chen and I.P. Lin.  The computational complexity of the satisfiability of modal Horn clauses for modal propositional logics.  Theoretical Computer Science Comput, 129(1):95–121, 1994.
  * [CN10] S. Cook and P. Nguyen, editors.  Logical foundations of proof complexity.  Cambridge University Press, 2010.
  * [DEL92] F. Debart, P. Enjalbert, and M. Lescot.  Multimodal logic programming using equational and order-sorted logic.  TCS, 105(1):141–166, 1992.
  * [FP87] L. Fariñas Del Cerro and M. Penttonen.  A note on the complexity of the satisfiability of modal Horn clauses.  Journal of Logic Programming, 4(1):1–10, 1987.
  * [GG13] D.M. Gabbay and F. Guenthner, editors.  Handbook of Philosophical Logic.  Springer Publishing Company, Incorporated, 2013.
  * [GHR94] D. Gabbay, I. Hodkinson, and M. Reynolds.  Temporal Logic: mathematical foundations and computational aspects.  Oxford University Press, 1994.
  * [Hor51] A. Horn.  On sentences which are true of direct unions of algebras.  Journal of Symbolic Logic, 16(1):14–21, 1951.
  * [HS91] J.Y. Halpern and Y. Shoham.  A propositional modal logic of time intervals.  Journal of the ACM, 38:279–292, 1991.
  * [Kro70] M.R. Krom.  The decision problem for formulas in prenex conjunctive normal form with binary disjunction.  Journal of Symbolic Logic, 35(2):14–21, 1970.
  * [Lad77] R.E. Ladner.  The computational complexity of provability in systems of modal propositional logic.  SIAM Journal of Computing, 6(3):467–480, 1977.
  * [Ngu00] L. A. Nguyen.  Constructing the least models for positive modal logic programs.  Fundamenta Informaticae, 42(1):29–60, 2000.
  * [Ngu03] L.A. Nguyen.  A fixpoint semantics and an SLD-resolution calculus for modal logic programs.  Fundamenta Informaticae, 55(1):63–100, 2003.
  * [Ngu04] L.A. Nguyen.  On the complexity of fragments of modal logics.  Advances in Modal Logic, 5:318–330, 2004.
  * [Ngu11] L.A. Nguyen.  Results on Modal Reasoning with Applications to Modal Deductive Databases.  PhD thesis, Warsaw University, 2011.
  * [Ngu16] L.A. Nguyen.  A tractable rule language in the modal and description logic that combines CPDL with regular grammar logic.  Fundam. Inform., 143(1–2):113–139, 2016.
  * [Pap03] C.H. Papadimitriou.  Computational complexity.  John Wiley and Sons Ltd., 2003.
  * [vHLP07] F. van Harmelen, V. Lifschitz, and B. Porter, editors.  Handbook of Knowledge Representation.  Elsevier Science, 2007.