# The Causal Marginal Polytope for Bounding Treatment Effects

Jakob Zeitler  University College London  Ricardo Silva  University College
London

###### Abstract

Due to unmeasured confounding, it is often not possible to identify causal
effects from a postulated model. Nevertheless, we can ask for _partial
identification_ , which usually boils down to finding upper and lower bounds
of a causal quantity of interest derived from all solutions compatible with
the encoded structural assumptions. One appealing way to derive such bounds is
by casting it in terms of a constrained optimization method that searches over
all causal models compatible with evidence, as introduced in the classic work
of Balke and Pearl (1994) for discrete data. Although by construction this
guarantees tight bounds, it poses a formidable computational challenge. To
cope with this issue, alternatives include algorithms that are not guaranteed
to be tight, or by introducing restrictions on the class of models. In this
paper, we introduce a novel alternative: inspired by ideas coming from belief
propagation, we enforce compatibility between _marginals_ of a causal model
and data, without constructing a _global_ causal model. We call this
collection of locally consistent marginals the _causal marginal polytope_. As
global independence constraints disappear when considering small dimensional
tractable marginals, this also leads to a rethinking of how to elicit and
express causal knowledge. We provide an explicit algorithm and implementation
of this idea, and assess its practicality with numerical experiments.

##  1 Introduction

In causal inference, we are concerned about identifying and estimating
properties of a distribution under an _intervention_. Using Pearl’s notation
[Pearl, 2009], a quantity such as $\mathbb{E}[Y\leavevmode\nobreak\
|\leavevmode\nobreak\ do(X=x)]$ describes the expected value of a random
variable $Y$ when $X$ is _controlled_ to take the value $x$. It differs from
the regression relationship $\mathbb{E}[Y\leavevmode\nobreak\
|\leavevmode\nobreak\ X=x]$, in which $X$ is a random variable and $X=x$ is a
particular realization of this random quantity. Those two function will be
equal only under some causal structures [Pearl, 2009, Spirtes et al., 2000],
which in general must come from background knowledge or causal discovery
assumptions that are not testable without experimental data. When hidden
common causes between treatment $X$ and outcome $Y$ cannot be blocked by
observed variables, there might be more than one, and perhaps uncountably
many, values for $\mathbb{E}[Y\leavevmode\nobreak\ |\leavevmode\nobreak\
do(X=x)]$ which are compatible with the observed data.

The problem of _partial identification_ amounts to providing the set of causal
quantities compatible with observable evidence. It will be a non-trivial set
under some background knowledge that rules or constraint the relationship
among variables. For instance, in the _instrumental variable_ setup we have
access to an auxiliary variable $Z$ that is associated with $X$ and causally
disconnected from $Y$ under an intervention on $X$. Under such a structure,
Balke and Pearl [1994] provided an early example on how to provide such a set
by finding an upper bound and a lower bound on quantities such as the average
treatment effect (ATE), $\mathbb{E}[Y\leavevmode\nobreak\
|\leavevmode\nobreak\ do(X=1)]-\mathbb{E}[Y\leavevmode\nobreak\
|\leavevmode\nobreak\ do(X=0)]$ for a binary $X$.

The approach gets tight bounds by parameterizing the causal model without
ruling out any possible structural equation that can causally link discrete
variables, and by solving the problem as a linear program. The number of
decision variables, however, quickly becomes unmanageable both in the number
of variables and number of categories each (discrete) variable can take,
motivating approaches which attempt to derive a subset of constraints [e.g.
Evans, 2012, Finkelstein and Shptiser, 2020].

In this paper, we contribute with a different direction that complements the
existing methods for partial identification in discrete distributions. In a
nutshell, it relaxes the problem of bounding causal effects by removing
constraints. The relaxation, however, follows a pattern that taps in the
familiar formulation of Balke and Pearl [1994] without having to invent new
sophisticated approaches for deriving such constraints. Instead, it takes a
leaf from generalized belief propagation Yedidia et al. [2005], concepts which
are already well-established in the uncertainty in artificial intelligence
literature. In particular, we show how to cast the problem in terms of small-
dimensional marginals of a causal model, enforcing coherence on overlapping
marginals. This forms what we call _the causal marginal polytope_. This
requires also thinking about the implications on how to express causal
knowledge, as conditional independence constraints disappear as we marginalize
a graphical model.

In Section 2 we describe the problem background, and provide further related
work. In Section 3, we describe the main construction of the causal marginal
polytope for discrete data, assuming first knowledge of true distributions
under particular observational and interventional regimes, and then discussing
a basic overview on how to use the method with finite samples. Section 4
illustrates the method with numerical examples

##  2 Model Setup and Related Work

We start reviewing the parameterization of discrete distributions according to
the setup of Balke and Pearl [1994]. This is contrasted with further
extensions and alternative approaches.

###  2.1 Parameterizing discrete causal models

A _structural causal model_ (SCM) [Pearl, 2009] describes a system of observed
variables $\mathbf{V}\equiv\\{V_{1},V_{2},\dots,V_{p}\\}$ as causally
connected by a set of _structural equations_ $\\{f_{1},f_{2},\dots,f_{p}\\}$
and _background_ variables $\mathbf{U}\equiv\\{U_{1},U_{2},\dots,U_{p}\\}$
following a distribution. That is, for each $V_{i}\in\mathbf{V}$, we have a
function $f_{i}$ such that

$V_{i}=f_{i}(V_{pa_{i}},U_{i}),$

where $V_{pa_{i}}\subset\mathbf{V}$ are called the (observable) parents of
$V_{i}$, in the sense that we can construct a directed acyclic graph (DAG)
among the elements of $\mathbf{V}$ with the corresponding set of parents
following the usual definitions from graph theory.

We assume there is a joint distribution $F_{\mathbf{U}}$ of $\mathbf{U}$,
keeping in mind that in principle each $U_{i}$ could be infinite-dimensional,
potentially forming a collection of stochastic processes. However, when all
variables in $\mathbf{V}$ are discrete, it is possible to assume without loss
of generality that each $U_{i}$ can be mapped to a particular value of a
scalar categorical variable $R_{i}$ with finitely many values, as originally
introduced by Balke and Pearl [1994]. In particular, considering the case
where all variables in $\mathbf{V}$ are binary, any structural causal model
can be represented as

$V_{i}=f_{R_{i}}(V_{pa_{i}})),$

where $R_{i}\in\\{1,2,\dots,2^{2^{|pa(V_{i})|}}\\}$ follows a categorical
distribution. The finite function space $\\{f_{r_{i}}\\}$ spans all functions
in $\\{0,1\\}^{|pa(V_{i})|}\rightarrow\\{0,1\\}$. The joint distribution of
$R_{1},R_{2},\dots,R_{p}$ can incorporate probabilistic constraints such as
marginal independence. It is commonly the case that we represent a model of
marginal independence using a _bidirected graphical model_ , where each pair
of vertices is either connected by a bidirected edge or has no edge. The
global Markov structure of such models describes that any two sets of vertices
$R_{S_{a}}$ and $R_{S_{b}}$, for which no bidirected edge exists across them,
should have a corresponding marginal independence structure, that is,
$p(r_{S_{a}},r_{S_{b}})=p(r_{S_{a}})p(r_{S_{b}})$. When combined with the DAG
structure implied by the structural equations, the result is an acyclic
directed mixed graph model (ADMG), with Markovian properties described by
Richardson [2003]. In this case, the ADMG with vertex set $\mathbf{V}$
includes bidirected edge $V_{i}\leftrightarrow V_{j}$ if and only if the
bidirected graph model for $\mathbf{R}$ contains the edge
$R_{i}\leftrightarrow R_{j}$.

The lack of a bidirected edge is commonly used to encode lack of unmeasured
confounding under particular conditioning sets. For instance, if the ADMG
implies that some set $V_{B}$ satisfies the backdoor criterion [Pearl, 2009]
for a pair $(V_{i},V_{j})$, then $p(v_{j}\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{i}),v_{B})=p(v_{j}\leavevmode\nobreak\
|\leavevmode\nobreak\ v_{i},v_{B})$.

Any interventional distribution can be calculated given $\mathcal{F}$, the
function spaces of structural equations for all observed variables, and the
distribution of the background distribution $p(r)$. For instance, if we divide
$\mathbf{V}$ into two sets $\mathbf{V}_{A}$ and $\mathbf{V}_{B}$, we have:

$\begin{array}[]{c}\displaystyle p(v_{A}^{\star}\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{B}^{\star}))=\\\
\displaystyle\sum_{r}p(r)\prod_{V_{i}\in\mathbf{V}_{A}}\mathbbm{1}(f_{r_{i}}(v_{pa_{i}})=v_{i}^{\star})\prod_{V_{j}\in\mathbf{V}_{B}}\mathbbm{1}(v_{j}=v_{j}^{\star}),\end{array}$
(1)

where $\mathbbm{1}(\cdot)$ is the indicator function. Information propagates
from parents to children given $r$ as the starting point, that is, the value
of vector $r$ deterministically implies $v$.

With $\mathcal{F}$ usually fixed, the only free parameter is the joint
distribution $p(r)$. Provided a collection $\mathcal{P}$ of $k$ distributions
in $\mathbf{V}$,

$\mathcal{P}\equiv\\{p(v_{\backslash D_{1}}\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{D_{1}})),\dots,p(v_{\backslash
D_{k}}\leavevmode\nobreak\ |\leavevmode\nobreak\ do(v_{D_{k}}))\\},$

where $\mathbf{V}_{D_{k}}\subset\mathbf{V}$ and $\mathbf{V}_{\backslash
D_{k}}\equiv\mathbf{V}\backslash\mathbf{V}_{D_{k}}$, we say that $p(r)$ is
_partially identifiable_ from $\mathcal{P}$ if it is a strict subset of all
distributions $\mathcal{R}$ that are allowed a priori, such as the space of
all distributions which are Markov with respect to some bidirected graph.

###  2.2 Related work

Restrictions in $\mathcal{F}$ and in $\mathcal{R}$ lead to restrictions in
$\mathcal{P}$. This can lead to model testability, meaning falsification of
the restrictions in $\mathcal{F}$ and $\mathcal{R}$. Early work includes
falsifying discrete instrumental variable models [Pearl, 1995]. Most
restrictions consist of independence constraints corresponding to lack of
directed or bidirected edges in the corresponding ADMG of the model, although
particular functional constraints such as monotonicity are also sometimes
used.

Many methods exist for deriving the constraints in $\mathcal{P}$ that are
implied by $\mathcal{F}$ and $\mathcal{R}$. Some common techniques include
deriving symbolically the polytope describing distributions in $\mathcal{P}$
using algebraic methods [Balke, 1995, Dawid, 2003], but the computational cost
of such methods is prohibitive for all but the smallest models. Alternatives
include non-exhaustive methods that infer a subset of constraints by a clever
combination of graph theory and algebra [e.g. Evans, 2012, Finkelstein and
Shptiser, 2020], or exhaustive and still expensive approaches that exploit
independence constraints (lack of edges) for classes of constraints such as
equalities [Tian and Pearl, 2002, Evans, 2018] and inequalities [Navascues and
Wolfe, 2020]. Deriving such constraints allows bounding of causal quantities,
such as the average treatment effect (ATE).

Numerical approaches express the problem directly by optimizing a quantity of
interest such as the ATE (or even a constant, if all we want is to test
feasibility) subject to equating the right-hand side of all instances of Eq.
(1) to (an estimate of) its left-hand side. When marginal independence
constraints are to be imposed in $p(r)$, parameterizations such as Drton and
Richardson [2008], or alternative ones such as Duarte et al. [2021], require
nonlinearities. As discussed by Duarte et al. [2021], however, such
parameterizations are polynomial programs still amenable to some guarantees of
optimality, and allow for the use of efficient implementations of mathematical
programming packages.

It still is the case that, even in sparse graphs, the number of decision
variables in the mathematical program are in general far too large. One
alternative is to constrain the function space $\mathcal{F}$. In particular,
one can postulate that we have a latent variable model with, say, one latent
variable per clique in the corresponding bidirected graph of background
variables $\mathbf{U}$, along with non-deterministic conditional probabilities
$p(v_{i}\leavevmode\nobreak\ |\leavevmode\nobreak\ v_{pa_{i}},u)$. The model
is unidentifiable, but a Bayesian posterior is still well-defined on its
parameter space [Chickering and Pearl, 1996]. The posterior, however, is on
the unidentifiable causal parameters and hence cannot be interpreted as
providing uncertainty on the bounds [Silva and Evans, 2016]. Although by
construction the posterior will always be within the bounds that an
optimization method can obtain (in the limit of infinite data), barring
algorithmic challenges such as the bad mixing of a Markov chain Monte Carlo
approach, it will in general be overconstrained, sometimes much narrower than
the result of a bounding procedure. How overconstrained will depend on the
prior distribution and, even if it correctly captures the estimand of
interest, it will be unwarranted if the prior does not follow from acceptable
background knowledge about the domain. The shape of the e.g. ATE posterior
within the bounds is completely independent of the data (the likelihood being
flat in that region), and will not be meaningful unless the prior is
informative [Silva and Evans, 2016].

Bounding methods for continuous variables will in general require constraints
in (the possibly infinite-dimensional) $\mathcal{F}$. See Gunsilius [2021,
2020], Kilbertus et al. [2020], Hu et al. [2021], Xia et al. [2021] for
examples. All of these methods rely heavily on Monte Carlo-powered
optimization, special algorithmic tricks to stabilize them, and particular
choices of function space. As we are specializing our approach to discrete
data, they are out of scope.

##  3 The Causal Marginal Polytope

Our approach is meant to complement ways of parameterizing $\mathcal{F}$ and
$\mathcal{R}$ in a way to scale up mathematical programs for bounding causal
effects. In particular, we will focus on the case where _there is no sparsity_
: we will not consider graphs with missing edges to avoid a discussion of ways
of parameterizing $p(r)$ [e.g., two alternatives are discussed by Drton and
Richardson, 2008, Duarte et al., 2021]. Instead, we will provide a way of
getting a fully linear program by constraining the “strength” of particular
edges as motivated by Ramsahai [2012] and Silva and Evans [2016]. We will
assume all variables are binary to simplify presentation. We also assume for
now that we know the probability distributions $\mathcal{P}$ from an oracle,
while in Section 3.3 we discuss what to do to estimate them from data.

###  3.1 Model Space and Data Binding

Let $\mathcal{M}$ be a collection of (conditional) _marginal causal models_
over subsets of set $\mathbf{V}$, variables which are assumed to follow a SCM
$\mathcal{M}_{global}$. We define a SCM $\mathcal{M}_{j}\in\mathcal{M}$ over a
subset $\mathbf{V}_{M_{j}}\subset\mathbf{V}$. Moreover, we also formulate the
framework by allowing conditioning on realizations of another subset
$\mathbf{V}_{M^{\prime}_{j}}\subset\mathbf{V}$, which contains only non-
descendants of $\mathbf{V}_{M_{j}}$ in the causal graph of
$\mathcal{M}_{global}$ (to be clear, in our numerical experimentation in
Section 4, we will not make use of conditioning sets $\mathbf{V}_{M_{j}}$ in
order to focus on other comparative aspects of the model space).

Model $\mathcal{M}_{j}$ is assumed to have a function space
$\mathcal{F}_{\mathcal{M}_{j}}$, which follows the same ancestral
relationships as those implied by $\mathcal{M}_{global}$, and a background
probability mass function $p_{\mathcal{M}_{j}}(r_{M_{j}}\leavevmode\nobreak\
|\leavevmode\nobreak\ v_{M^{\prime}_{j}})$. The background pmf is
parameterized by $\theta_{j.v_{M^{\prime}_{j}}}$, as a contingency table for
each instantiation $v_{M^{\prime}_{j}}$ that we include in $\mathcal{M}_{j}$.
Not all possible instances of the sample space of
$\mathbf{V}_{M^{\prime}_{j}}$ need to be defined by $\mathcal{M}_{j}$: a
reduced number of values $v_{M^{\prime}_{j}}$ can be selected according to
computational bottlenecks at the expense of introducing further relaxations to
the mathematical program. We call the set of possible values
$v_{M^{\prime}_{j}}$ defined by $\mathcal{M}_{j}$ its _conditional scope_.

In the examples of Section 4, the function space
$\mathcal{F}_{\mathcal{M}_{j}}$ is exhaustive with respect to
$\mathcal{M}_{j}$, i.e., it contains all possible corresponding binary
structural equations, as in Balke and Pearl [1994], among the causally-ordered
elements of $\mathbf{V}_{M_{j}}$. Any variable in
$\mathbf{V}_{M^{\prime}_{j}}$ in $\mathcal{M}_{j}$ will play a role of
“background” variable similar to the original definition of $\mathbf{U}$. This
means that the (linear) mapping from
$(\mathcal{F}_{\mathcal{M}_{j}},\theta_{j.v_{M^{\prime}_{j}}})$ to
$p_{\mathcal{M}}(v_{\backslash D_{M_{j}}}\leavevmode\nobreak\
|\leavevmode\nobreak\ v_{M^{\prime}_{j}},do(v_{D_{M_{j}}}))$,
$\mathbf{V}_{D_{M_{j}}}\subset\mathbf{V}_{M_{j}}$, $\mathbf{V}_{\backslash
D_{M_{j}}}\equiv\mathbf{V}_{M_{j}}\backslash\mathbf{V}_{D_{M_{j}}}$, can be
done in an entirely analogous way to Eq. (1), with
$p_{\mathcal{M}}(r_{M_{j}}\leavevmode\nobreak\ |\leavevmode\nobreak\
v_{M^{\prime}_{j}})$ taking the role of the averaging distribution.

The collection of marginals $\mathcal{M}$ forms a _realizable_ set of
distributions if there exists a collection of joint distributions of which
they are marginals. When we enforce that each
$p_{\mathcal{M}}(\cdot\leavevmode\nobreak\ |\leavevmode\nobreak\ \cdot)$ in
$\mathcal{M}$ must match the corresponding $p(\cdot\leavevmode\nobreak\
|\leavevmode\nobreak\ \cdot)$ in $\mathcal{P}$, we guarantee that they are
realizable for the regimes in $\mathcal{P}$. The fact that this collection is
indexed by different regimes forms a generalization of the _marginal polytope_
construction [Wainwright and Jordan, 2008], as usually exploited in the
variational inference literature. We will call our collection the _causal
marginal polytope_. As in belief propagation and their generalizations
[Yedidia et al., 2005], we will relax the problem of defining the set of all
causal models compatible with observation by enforcing only _local coherence_
over pairwise marginal distributions in $\mathcal{M}$: submarginals jointly
implied by pair $(\mathcal{M}_{j_{a}},\mathcal{M}_{j_{b}})$ must agree in each
possible causal regime. In particular, they must agree on regimes that are not
in $\mathcal{P}$, as those will be directly matched to each $\mathcal{M}_{j}$,
a step made more explicit in the next section.

###  3.2 Marginalization Constraints and Implications to Modeling

Our linear program formulation contains the constraint that

  * •

for each $\mathcal{M}_{j}\in\mathcal{M}$, we have that
$\theta_{j.{v_{M^{\prime}_{j}}}}$ must be in the probability simplex;

  * •

for each $\mathcal{M}_{j}\in\mathcal{M}$ where $v_{M^{\prime}_{j}}$ is in its
conditional scope, we have that

$p_{\mathcal{M}_{j}}(v_{\backslash D_{M_{j}}}\leavevmode\nobreak\
|\leavevmode\nobreak\ v_{M^{\prime}_{j}},do(v_{D_{M_{j}}}))=$ $p(v_{\backslash
D_{M_{j}}}\leavevmode\nobreak\ |\leavevmode\nobreak\
v_{M^{\prime}_{j}},do(v_{D_{M_{j}}})),$

if $p(v_{\backslash D_{M_{j}}}\leavevmode\nobreak\ |\leavevmode\nobreak\
v_{M^{\prime}_{j}},do(v_{D_{M_{j}}}))$ is contained in $\mathcal{P}$.

Now, consider a pair of causal marginals $\mathcal{M}_{j_{a}}$ and
$\mathcal{M}_{j_{b}}$ where $v_{c}$ is in the conditional scope of both
models, and where $\mathbf{V}_{M_{j_{a}}}$, $\mathbf{V}_{M_{j_{b}}}$ are the
corresponding causal variables in each marginal. Let
$\mathbf{V}_{M_{j_{ab}}}\equiv\mathbf{V}_{M_{j_{a}}}\cap\mathbf{V}_{M_{j_{b}}}$.
Likewise, let $\mathbf{V}_{D_{M_{j_{ab}}}}\subset\mathbf{V}_{M_{j_{ab}}}$,
with $\mathbf{V}_{\backslash
D_{M_{j_{ab}}}}\equiv\mathbf{V}_{M_{j_{ab}}}\backslash V_{D_{M_{j_{ab}}}}$
defined accordingly. We call _local coherence_ constraints those of the type

$p_{\mathcal{M}_{j_{a}}}(v_{\backslash D_{M_{j_{ab}}}}\leavevmode\nobreak\
|\leavevmode\nobreak\ v_{M^{\prime}_{j}},do(v_{D_{M_{j_{ab}}}}))=$
$p_{\mathcal{M}_{j_{b}}}(v_{\backslash D_{M_{j_{ab}}}}\leavevmode\nobreak\
|\leavevmode\nobreak\ v_{M^{\prime}_{j}},do(v_{D_{M_{j_{ab}}}})).$

These constraints will be non-redundant if $p(v_{\backslash
D_{M_{j_{ab}}}}\leavevmode\nobreak\ |\leavevmode\nobreak\
v_{M^{\prime}_{j}},do(v_{D_{M_{j_{ab}}}}))$ is not in $\mathcal{P}$.

One implication of working with marginal constraints, besides being a
relaxation that trades-off tractability for looser bounds, is that some
conditional independencies cannot be explicitly represented. For instance, if
we postulate that some observable $X_{1}$ removes confounding between some
$X_{2}$ and some $X_{3}$, then a model that excludes $X_{1}$ while including
$X_{2}$ and $X_{3}$ will need to express unmeasured confounding between the
latter two.

The upside, and an opportunity, is that domain experts may feel more
comfortable about expressing knowledge on the _strength_ of direct causal
effects and unmeasured confounding on smaller marginals of a causal system,
than on a full causal structure among all variables. The idea of composite
likelihood, for instance, is to express only particular marginals of a
likelihood function. Although in machine learning that is mostly motivated by
computational tractability, one appeal in statistical modeling is that this
reduces exposure to model misspecification [Varin et al., 2011].

In the context of partial identification methods, Ramsahai [2012] and Silva
and Evans [2016] discuss ways by which “weak directed effects” and “weak
unmeasured confounding” can be encoded. Those constraints were introduced in
the context of non-counterfactual approaches for instrumental variable models
and are hard to generalize for other structures. Here, we show a way of
encoding them in the context of Balke and Pearl’s parameterization.

#### “Weak” direct effects.

Assume that we elicit from an expert or algorithm that the direct effect of
some $X_{j}$ into $X_{i}$ is “weak,” within the system covered by marginal
model $\mathcal{M}$ and conditional value
$V_{\mathcal{M}^{\prime}}=v_{\mathcal{M}^{\prime}}$. This knowledge must be
stated as

$|P_{\mathcal{M}}(V_{i}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{pa_{i}\backslash j}),do(V_{j}=1),v_{\mathcal{M}^{\prime}})-$
$P_{\mathcal{M}}(V_{i}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{pa_{i}\backslash
j}),do(V_{j}=0),v_{\mathcal{M}^{\prime}})|\leq\epsilon_{ij},$

for a particular parent $V_{j}$ of $V_{i}$ in marginal model $\mathcal{M}$,
where $\mathbf{V}_{pa_{i}\backslash j}$ is the set of parents of $V_{i}$ other
than $V_{j}$ and $\epsilon_{ij}$ is a parameter to be elicited from an expert.
See Silva and Evans [2016] for a discussion on its choice, including its
derivation from generalizations of structure learning methods, which is out of
scope of this manuscript. As $P_{\mathcal{M}}(V_{j}=1\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{pa_{i}}))$ is a linear function of the parameters
of $\mathcal{M}$, per Eq. (1), and the above is a linear function of
$P_{\mathcal{M}}(V_{j}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{pa_{i}}))$, the addition of such constraints still imply a linear
programming formulation if the objective function is also linear (which will
be the case for the ATE).

Such constraints are expressed as interventional, as opposed to
counterfactual: they are testable from experiments even if the
parameterization could be interpreted as also implying the existence of cross-
world potential outcomes, a criticism of Balke and Pearl [1994] raised by
[Ramsahai, 2012]. The case $\epsilon_{ij}=0$ corresponds to a conditional
independence constraint in the experimental distribution where we
simultaneously intervene on all parents of $V_{i}$, and a common interpretion
of the lack of an directed edge in a causal DAG.

#### “Weak” ancestral bidirected chains.

Start by assuming that we elicit from an expert or algorithm that the level of
unmeasured confounding between some some ancestor $X_{j}$ of $X_{i}$, and
$X_{i}$ itself, is “weak” within the system covered by marginal model
$\mathcal{M}$ and conditional value
$V_{\mathcal{M}^{\prime}}=v_{\mathcal{M}^{\prime}}$. We first consider
constraints of the type

$|P_{\mathcal{M}}(V_{i}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{pa_{ij}}),do(V_{j}=v_{j}),v_{\mathcal{M}^{\prime}})-$
$P_{\mathcal{M}}(V_{i}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{pa_{ij}}),V_{j}=v_{j},v_{\mathcal{M}^{\prime}})|\leq\epsilon^{\mathsf{C}}_{ij},$

for $v_{j}=0,1$, and where $\mathbf{V}_{pa_{ij}}$ is the union of the parents
of $V_{i}$ and $V_{j}$ in $\mathcal{M}$, other than $\\{V_{i},V_{j}\\}$
themselves. Once again, $\epsilon^{\mathsf{C}}_{ij}$ is part of the model
input specification.

To avoid formulating a non-linear constraint, we could first attempt to
express the second term in the above as

$\displaystyle\frac{P_{\mathcal{M}}(V_{i}=1,V_{j}=v_{j}\leavevmode\nobreak\
|\leavevmode\nobreak\
do(v_{pa_{ij}}),v_{\mathcal{M}^{\prime}})}{P(V_{j}=v_{j}\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{pa_{j}}),v_{\mathcal{M}^{\prime}})},$

That is, the denominator can be treated as a constant, _if_ the regime where
we intervene in $V_{pa_{j}}$ is in the data.

For instance, assume that the graph of the marginal is $V_{1}\rightarrow
V_{2}\rightarrow V_{3}$, with “unconstrained bidirected edges”
$\\{V_{1}\leftrightarrow V_{2},V_{1}\leftrightarrow V_{3}\\}$ and a “weak”
amount of unmeasured confounding left between $V_{2}$ and $V_{3}$. Assume that
we have in $\mathcal{P}$ the regimes corresponding to $do(X_{1}=1)$ and
$do(X_{1}=0)$, and that the conditional scope is empty. We then express a
“weak” $V_{2}\leftrightarrow V_{3}$ as

$|P_{\mathcal{M}}(V_{3}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{1}),do(V_{2}=v_{2}))-$ $P_{\mathcal{M}}(V_{3}=1\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{1}),V_{2}=v_{2})|\leq\epsilon^{\mathsf{C}}_{ij}.$

The second term can be expressed as

$\displaystyle\frac{P(V_{3}=1,V_{2}=v_{2}\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{1}))}{P(V_{2}=v_{2}\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{1}))},$

which is given by $\mathcal{P}$, and hence a constant, making the constraint
linear in the parameters. However, if regime $do(v_{1})$ is _not_ in
$\mathcal{P}$, then the denominator is a function
$P_{\mathcal{M}}(V_{2}=v_{2}\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{1}))$ of the parameters, and the constraint becomes _polynomial_ , as in
Duarte et al. [2021].

However, once we give up on hard unmeasured confounding constraints (i.e.,
lack of bidirected edges), we may instead enquire knowledge about the strength
of _bidirected paths_ that confound $V_{i}$ and $V_{j}$. Consider the
constraint

$|P_{\mathcal{M}}(V_{i}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(v_{pa_{i}}),do(V_{j}=v_{j}),V_{\prec j}=v_{\prec
j},v_{\mathcal{M}^{\prime}})-$ $P_{\mathcal{M}}(V_{i}=1\leavevmode\nobreak\
|\leavevmode\nobreak\ do(v_{pa_{i}}),V_{j}=v_{j},V_{\prec j}=v_{\prec
j},v_{\mathcal{M}^{\prime}})\leq\epsilon^{\mathsf{C}}_{ij},$

for $v_{j}\in\\{0,1\\}$, where $\mathbf{V}_{\prec_{j}}$ are the ancestors
$V_{j}$ in $\mathcal{M}$. For $\epsilon^{\mathsf{C}}_{ij}=0$, this would _not_
necessarily correspond to the lack of bidirected edge $V_{j}\leftrightarrow
V_{i}$. In the previous example, without the edge $V_{1}\rightarrow V_{3}$,
expressing $P(V_{3}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(V_{2}=v_{2}))\approx P(V_{3}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
V_{2}=v_{2})$ captures the idea of weak confounding. However, if edge
$V_{1}\rightarrow V_{3}$ exists, then we are allowing for a possible backdoor
path $V_{3}\leftarrow V_{1}\rightarrow V_{2}$. An alternative is check whether
input knowledge allows for

$P(V_{3}=1\leavevmode\nobreak\ |\leavevmode\nobreak\
do(V_{2}=v_{2}),V_{1}=v_{1})\approx$ $P(V_{3}=1\leavevmode\nobreak\
|\leavevmode\nobreak\ V_{2}=v_{2},V_{1}=v_{1}).$

This would correspond to a constraint on the contribution of the active
backdoor path $V_{3}\leftrightarrow V_{1}\leftrightarrow V_{2}$ given $V_{1}$,
a “weak bidirected chain” that passes through the ancestors of $V_{j}$. The
convenience of this constraint is that $P(V_{3}=1\leavevmode\nobreak\
|\leavevmode\nobreak\ V_{2}=v_{2},V_{1}=v_{1})$ is in $\mathcal{P}$, and hence
can be treated as a constant, making the above constraint linear in the
parameters.

Its interpretation is that the ancestors of $V_{j}$ are taken as a
“approximate covariate adjustment set,” that is, those ancestors would
approximately block confounders between $X_{i}$ and $X_{j}$. Notice that
confounding induced by vertices “downstream” of $X_{j}$ are accounted by the
intervention on the parents of $X_{i}$. That is, if the structure is
$X_{2}\rightarrow X_{3}\rightarrow X_{4}$, with bidirected edges
$X_{2}\leftrightarrow X_{3}$, and $X_{2}\leftrightarrow X_{4}$, then a
statement of “weak confounding” between $X_{2}$ and $X_{4}$ removes the
confounding contribution of the path $X_{2}\leftrightarrow X_{3}\rightarrow
X_{4}$ via the intervention $do(x_{3})$. If desired, we may wish to condition
on a subset of the ancestors of $V_{j}$ as opposed to $V_{\prec j}$.

###  3.3 Learning

In our illustrative examples of the next section, our applications will be
simple enough so that each element of $\mathcal{P}$ can be estimated just
using the frequencies in the data. We do not discuss measures of uncertainty.
One idea is discussed by Duarte et al. [2021], treating elements of
$\mathcal{P}$ themselves as unknown, and taking values within a confidence
interval obtained by a black-box fitting of the marginals. Likewise, the
Bayesian approach Silva and Evans [2016] can be immediately used here: a
sample from the posterior of a black-box model is passed through the
optimization problem (which is, after all, just a functional of $\mathcal{P}$)
resulting in a posterior sample of an objective function of interest.

###  3.4 Summary

To summarize, exploiting a causal marginal polytope has major computational
advantages when we want to avoid constraining the function space of structural
equations while retaining tractability. The downsides are: i) it is a
relaxation, which in general will imply looser bounds - but unless we are
willing to introduce untestable constraints in the function space, it is a
price we need to pay; ii) given that marginals imply a removal of conditional
independencies, we need to rethink which structural knowledge should be
enquired. To maintain linearity of the mathematical program, we suggest weak
directed edges and weak bidirected chains. This is not necessarily a downside:
we may never believe sparsity anyway, even in the full causal system, and it
may actually be easier to extract knowledge from enquiring experts about small
subsystems than a complex causal system of many variables. If such a
collection of marginal models is inconsistent, in the sense that no full
causal model is compatible with them, we can still detect where local
coherence breaks down as the linear program solver should flag infeasibility.

##  4 Experiments

Our goal in this section is to provide numerical examples of an implementation
of the causal marginal polytope in action. Training data is considered to be
large enough so that, as the optimization is exact, whether we trap or not the
true causal effect will be a function of the input knowledge. Hence, we
illustrate how informative (or not) bounds turn out to be, as we vary
assumptions about weaknesses of edges, as well as our choice of marginals.

One of the primary goals of this section is to highlight how knowledge
elicitation in this marginal setting may lead to logical but potentially
counter-intuitive results. For instance, it is to be expected that bounds
obtained by using 3-dimensional marginals of a 4-dimensional causal model
would be at least as wide as, and in fact contain, the bounds found by running
a linear program on the 4-variable model. However, this is predicated on the
3-dimensional models being submodels of the corresponding marginal in the
4-variable global model. This would be the case if knowledge were to be
specified only in terms of the global Markov properties. However, besides the
computational gains, much of the modeling motivation for this work is that
sparsity constraints are not necessarily desirable, and that knowledge
elicitation from margins is an asset to be added to the modeling toolkit of a
practitioner. Model spaces are a function of the restriction parameters
$\\{\epsilon_{ij},\epsilon_{ij}^{\mathsf{C}}\\}$. We should not expect that
human expert judgement calls will entail that restriction parameters chosen
from a marginal perspective will imply a larger feasible region to the linear
program compared to restriction parameters in the global model elicited from
the same individual (if this was even possible to compute). Being allowed to
think in terms of margins is a contribution that we believe will be of
practical value for experts.

###  4.1 Setup

Data As our starting point for
$\mathcal{M}_{global}=\\{X_{1},X_{2},X_{3},X_{4}\\}$ we assume the DAG in
Figure 1, i.e. a fully connected graph with full confounding. We construct a
simulation from which we can sample a single example SCM with random
structural equations. We provide the observational data regime as well as a
subset of all possible single and double intervention data regimes to our
bounding model:
$do(\emptyset),do(x_{2}),do(x_{3}),do(x_{2},x_{3}),do(x_{1}=0,x_{3}=0),do(x_{1}=0,x_{3}=1)$
As such, some interventions will be identifiable, and some will be partially
informed from overlapping interventional information or the observational
regime, see Figure 2.

Figure 1: $\mathcal{M}_{global}$: The red bold directed and bidirected edges
will be constrained by expert knowledge about the causal system. Figure 2:
Unconstrained bounds on all single and double interventions. True effects are
marked with a star. Intervention 3 to 12 are identified from just data alone.

Model We choose margins of length 3 from $\mathcal{M}_{global}$, i.e. $M=\\{$
$M_{1}=\\{X_{1},X_{2},X_{3}\\},$ $M_{2}=\\{X_{1},X_{2},X_{4}\\},$
$M_{3}=\\{X_{1},X_{3},X_{4}\\},$ $M_{4}=\\{X_{2},X_{3},X_{4}\\}\\}$ We enforce
the overlap constraint on the overlap sets
$\\{\\{X_{1},X_{2}\\},\\{X_{2},X_{4}\\},\\{X_{1},X_{3}\\},\\{X_{3},X_{4}\\}\\}$.
Given assumptions from Figure 1, we then enforce a weak edge $X_{1}\rightarrow
X_{4}$ and edge $X_{1}\leftrightarrow X_{4}$ in margins $M_{2}$ and $M_{3}$,
as $M_{1}$ and $M_{4}$ do not contain both $X_{1}$ and $X_{4}$.

###  4.2 Results

We will first consider the impact of the constraints individually, then
together. True intervention effects are shown as dotted lines. Interventions
that are not affected by the constraints are not shown. Figure 2 shows the
bounds and true intervention effects for all possible single and double
interventions of the unconstrained causal marginal polytope.

Overlap Figure 3 shows the impact of the overlap constraints on bounds, i.e.
in this case tightening the lower bound on $P(X_{4}=1|do(X_{1}=0,X_{2}=1)$
from 0.3 to 0.4, and $P(X_{4}=1|do(X_{1}=0)$ from 0.3 to 0.6.

Weak directed edge Figure 4 shows the impact of a decreasing $\epsilon$ for
constraining $X_{1}\rightarrow X_{4}$. Notably, $\epsilon\leq 0.03$ does not
almost achieve identification of $P(X_{4}=1|do(X_{1}=0,X_{3}=0)$ and
$P(X_{4}=1|do(X_{1}=0,X_{3}=1)$, as the crosses indicate invalid upper bounds
that were tightened below the true causal effect (dotted line). This
demonstrates an $\epsilon$ choice assuming a too weak $X_{1}\rightarrow
X_{4}$.

Weak bidirected edge Figure 5 shows the impact of assuming an $\epsilon$-weak
$X_{1}\leftrightarrow X_{4}$. The linear program is infeasible from
$\epsilon<0.12$, i.e. using the given data the model falsifies the assumption
of such a strong choice of $\epsilon$.

Overlap with weak directed and bidirect edges Figure 6 finally shows all three
constraints together. The single intervention $P(X_{4}=1|do(X_{1}=1)$ has been
slightly tightened further compared to Figure 5 while $P(X_{4}=1|do(X_{1}=0)$
is not shown anymore, i.e. the combination of the constraints tightened it to
a space in the polytope where choices of $\epsilon\geq 0.12$ do not affect it.
For simplicity of exposition, we set $\epsilon$ to the same number for each
weak edge constraint, but we can choose $\epsilon$ individually for each weak
edge (e.g. $\epsilon^{(X_{1}\rightarrow X_{4})}$ and
$\epsilon^{(X_{1}\leftrightarrow X_{4})}$). Figure 7 shows an example where we
set $\epsilon^{(X_{1}\rightarrow X_{4})}=0.06$ and vary
$\epsilon^{(X_{1}\leftrightarrow X_{4})}$, leading to more tightening, but
also two invalid bounds at $\epsilon=0.12$.

Figure 3: Tightening of bounds by constraining the overlap of the margins.
Figure 4: A weak directed edge tightens bounds. Invalid bounds are marked with
a cross. Figure 5: A weak bidirected edge tightens bounds. Lower epsilon
values are falsified by infeasability. Figure 6: Weak directed and bidirected
edges together with overlap constraint. Figure 7: Bounds with all constraints
in the model, with $\epsilon^{(X_{1}\rightarrow X_{4})}=0.06$ and
$\epsilon^{(X_{1}\leftrightarrow X_{4})}$ decreasing x-axis. Figure 8: All
constraints on the full margin of $\mathcal{M}_{global}$ in Figure 1.

Comparing with full parameterisation Figure 8 shows a causal marginal polytope
with a single margin of size 4, i.e. building a model directly from
$\mathcal{M}_{global}$. With $2^{15}$ parameters it is generally the only full
parameterisation of bounding model that is still computationally feasible. At
a first glance, the bounds seem more loose compared to the causal marginal
polytope. Importantly, such a comparison needs to be carefully done as choices
of $\epsilon$ have unique interpretations for each margin, including the full
margin $\mathcal{M}_{global}$, as mentioned in the introduction to this
Section.

###  4.3 Experiments with $n>4$

For a $n=6$ SCM , Figure 9 shows bounds on interventions assuming
$\epsilon$-weak edges $X_{1}\rightarrow X_{4}$, $X_{1}\rightarrow X_{6}$,
$X_{3}\rightarrow X_{5}$, $X_{1}\leftrightarrow X_{4}$ and
$X_{2}\leftrightarrow X_{6}$. We provide the same data regimes as in the
previous example. A comparison with a full margin model is computationally
infeasible as it would require a model with $2^{31}$ parameters, compared to
the computationally feasible $2^{1+2+4}\times 20=2560$ parameters of our
causal marginal polytope with 20 margins each $|M_{i}|=3$. In practice, an SCM
with 10 variables will require only 15360 parameters and, given careful choice
of the $\epsilon$-weak edges, can yield informative bounds. Importantly, the
examples above show causal systems where all three constraints are relevant to
the margins, which is not always the case. Furthermore, the available data
impacts the tightness of the bounds, with double intervention regimes being
more useful for informing bounds on unobserved single intervention effects.
Finally, to achieve the tightest bounds on the causal marginal polytope, a
thorough evaluation should always individually choose $\epsilon$ for each
constraining of weak directed and bidirected edges.

Figure 9: Bounds on 6 variable graph.

##  5 Conclusion

We introduced the causal marginal polytope to address two needs: computational
scalability and an alternative way to elicit causal knowledge that will not
require thinking in terms of a complete system. There are several uses which
were not considered in this manuscript. For instance, it can guide
experimental design where a small combination of design points can lead to
information about combinations of promising smaller-dimensional (and hopefully
cheaper) interventions. How to extend this idea to continuous variables, as a
way of reducing variance in the complex Monte-Carlo based methods of e.g.
Kilbertus et al. [2020], Hu et al. [2021], Xia et al. [2021], is another
direction of research.

## References

  * Balke [1995] A. Balke.  _Probabilistic Counterfactuals: Semantics, Computation, and Applications_.  1995\.
  * Balke and Pearl [1994] A. Balke and J. Pearl.  Counterfactual probabilities: Computational methods, bounds and applications.  _Proceedings of the Tenth Conference on Uncertainty in Artificial Intelligence (UAI1994)_ , pages 46–54, 1994.
  * Chickering and Pearl [1996] D. Chickering and J. Pearl.  A clinician’s apprentice for analyzing non-compliance.  _Proceedings of the 12th National Conference on Artificial Intelligence_ , pages 1268–1276, 1996.
  * Dawid [2003] A.P. Dawid.  Causal inference using influence diagrams: the problem of partial compliance.  In P.J. Green, N.L. Hjort, and S. Richardson, editors, _Highly Structured Stochastic Systems_ , pages 45–65. Oxford University Press, 2003.
  * Drton and Richardson [2008] M. Drton and T. Richardson.  Binary models for marginal independence.  _Journal of the Royal Statistical Society, Series B_ , 70:287–309, 2008.
  * Duarte et al. [2021] G. Duarte, N. Finkelstein, D. Knox, J. Mummolo, and I. Shpitser.  An automated approach to causal inference in discrete settings.  _arXiv:2109.13471_ , 2021.
  * Evans [2012] R. Evans.  Graphical methods for inequality constraints in marginalized DAGs.  _Proceedings of the 22nd Workshop on Machine Learning and Signal Processing_ , 2012.
  * Evans [2018] R. Evans.  Margins of discrete Bayesian networks.  _Annals of Statistics_ , 46:2623–2656, 2018.
  * Finkelstein and Shptiser [2020] N. Finkelstein and I. Shptiser.  Deriving bounds and inequality constraints using logical relations among counterfactuals.  _Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), PMLR_ , 124, 2020.
  * Gunsilius [2020] F. Gunsilius.  A path-sampling method to partially identify causal effects in instrumental variable models.  _arXiv:1910.09502_ , 2020.
  * Gunsilius [2021] F. Gunsilius.  Nontestability of instrument validity under continuous treatments.  _Biometrika_ , 108(4):989–995, 2021.
  * Hu et al. [2021] Y. Hu, Y. Wu, L. Zhang, and X. Wu.  A generative adversarial framework for bounding confounded causal effects.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35:12104–12112, 2021.
  * Kilbertus et al. [2020] N. Kilbertus, M.J . Kusner, and R. Silva.  A class of algorithms for general instrumental variable models.  _Advances in Neural Information Processing Systems_ , 33:20108–20119, 2020.
  * Navascues and Wolfe [2020] M. Navascues and E. Wolfe.  The inflation technique completely solves the causal compatibility problem.  _Journal of Causal Inference_ , 8:70–91, 2020.
  * Pearl [1995] J. Pearl.  On the testability of causal models with latent and instrumental variables.  _Proceedings of the 11th Conference on Uncertainty in Artificial Intelligence (UAI ’95)_ , pages 435–443, 1995.
  * Pearl [2009] J. Pearl.  _Causality: Models, Reasoning and Inference, 2nd edition_.  Cambridge University Press, 2009.
  * Ramsahai [2012] R. Ramsahai.  Causal bounds and observable constraints for non-deterministic models.  _Journal of Machine Learning Research_ , pages 829–848, 2012.
  * Richardson [2003] T. Richardson.  Markov properties for acyclic directed mixed graphs.  _Scandinavian Journal of Statistics_ , 30:145–157, 2003\.
  * Silva and Evans [2016] R. Silva and R. Evans.  Causal inference through a witness protection program.  _Journal of Machine Learning Research_ , 17(56):1–53, 2016.
  * Spirtes et al. [2000] P. Spirtes, C. Glymour, and R. Scheines.  _Causation, Prediction and Search_.  Cambridge University Press, 2000.
  * Tian and Pearl [2002] J. Tian and J. Pearl.  On the testable implications of causal models with hidden variables.  _Proceedings of the 18th Conference on Uncertainty in Artificial Intelligence_ , 2002.
  * Varin et al. [2011] C. Varin, N. Reid, and D. Firth.  An overview of composite likelihood methods.  _Statistica Sinica_ , 21:5–42, 2011.
  * Wainwright and Jordan [2008] M. Wainwright and M. Jordan.  Graphical models, exponential families, and variational inference.  _Foundations and Trends in Machine Learning_ , 1:1–305, 2008.
  * Xia et al. [2021] K. Xia, K.-Z. Lee, Y. Bengio, and E. Bareinboim.  The causal-neural connection: Expressiveness, learnability, and inference.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Yedidia et al. [2005] J. Yedidia, W. Freeman, and Y. Weiss.  Constructing free-energy approximations and generalized belief propagation algorithms.  _IEEE Transactions on Information Theory_ , 51:2282–2312, 2005.


# Revisiting the Siegert relation for the partially coherent regime of
nanolasers

M. Drechsler  Institute for Theoretical Physics, University of Bremen, Bremen,
Germany  Bremen Center for Computational Materials Science, University of
Bremen, Bremen, Germany  F. Lohof  Institute for Theoretical Physics,
University of Bremen, Bremen, Germany  Bremen Center for Computational
Materials Science, University of Bremen, Bremen, Germany  C. Gies
<EMAIL_ADDRESS>http://www.itp.uni-bremen.de/ag-gies/  Institute for
Theoretical Physics, University of Bremen, Bremen, Germany  Bremen Center for
Computational Materials Science, University of Bremen, Bremen, Germany

###### Abstract

The Siegert relation connects the first- and second-order coherence properties
of light. While strictly valid only in the thermal regime and in the absence
of correlations, this relation is routinely extended to the partially coherent
regime in the study of high-$\beta$ nanolasers, where it aids in the
identification of the lasing threshold. We re-evaluate the use of a
generalized Siegert relation in different device regimes. A full two-time
quantum-optical theory is derived as a reference for obtaining first- and
second-order correlation functions $g^{(1)}(\tau)$ and $g^{(2)}(\tau)$ in the
steady state. We find that even in the partially coherent regime the
generalized Siegert relation is well suited as an approximation to
$g^{(2)}(\tau)$ as long as emitter correlation are negligible, but does not
apply well in the quantum regime of few-emitter nanolasers, or to devices
featuring sub- and superradiant emission.

††preprint: AIP/123-QED

Photon-correlation functions are a powerful tool in the analysis of light, as
they carry signatures of the light’s photon statistics, its spatial and
temporal coherence, and they contain information on the emission process
itself Loudon (2000). In particular, the first-order coherence function
$g^{(1)}(\tau)$ is related to temporal fluctuations in the electric field and
the spectral properties of the light, while the second-order coherence
function $g^{(2)}(\tau)$ contains information on intensity fluctuations
Glauber (1963). These coherence functions are generally an accepted measure of
the coherence properties in the emission of nanolaser devices. Nevertheless,
it should be said that small lasers may exhibit peculiar properties that arise
from many-body effects, correlations, and fluctuations that go beyond what is
commonly captured in the first and second moments of the photon statistics
Wang, Puccioni, and Lippi (2015); Roy-Choudhury, Haas, and Levi (2009);
Schlottmann  _et al._ (2018). Analysis of second-order correlations is used in
cavity quantum electrodynamics (cQED) Hamsen  _et al._ (2017); Snijders  _et
al._ (2018) and for the study of emission characteristics and noise spectra of
LED and laser light sources Paik  _et al._ (2019); Megyeri  _et al._ (2018);
Shan  _et al._ (2021); Blazek  _et al._ (2011); Wang  _et al._ (2021). In the
study of high-$\beta$ nanolasers, measurements of first- and second-order
coherence are commonly used to reveal the transition from thermal to coherent
emission as they carry fingerprints of the corresponding photon statistics
Jagsch  _et al._ (2018); Ma and Oulton (2019); Kreinberg  _et al._ (2020);
Ulrich  _et al._ (2007); Takemura  _et al._ (2019); Mork and Lippi (2018).

The Siegert relation Siegert (1943) relates field and intensity correlations
for light that exhibit a Gaussian field-amplitude distribution via

$\displaystyle g^{(2)}(\tau)=1+|g^{(1)}(\tau)|^{2}.$ (1)

Its classical derivation assumes a large ensemble of emitters and uncorrelated
emission events Voigt and Hess (1994); Ferreira  _et al._ (2020). However, a
quantum mechanical treatment relies on the Wick’s theorem for systems in
thermal equilibrium and weak interaction between light and matter Fricke
(1996) (see supplementary material). The Siegert relation has been employed to
analyze experiments in astronomy Foellmi (2009), molecular physics Borycki,
Kholiqov, and Srinivasan (2016), biology Stetefeld, McKenna, and Patel (2016),
and quantum optics Aßmann  _et al._ (2010). In particular, it has become an
established tool in the characterization of the laser threshold in nanolasers
Ulrich  _et al._ (2007); Chow, Jahnke, and Gies (2014).

In studies of high-$\beta$ semiconductor nanolasers, the transition from
thermal to coherent emission is commonly detected by measurements of the
second-order correlation function $g^{(2)}(\tau)$ with a Hanbury-Brown and
Twiss (HBT) setup that uses avalanche photo diodes or superconducting single-
photon detectors (SSPD) for recording correlated photon-emission events
Kreinberg  _et al._ (2020, 2017); Takemura  _et al._ (2019); Elvira  _et al._
(2011); Yokoo  _et al._ (2017). The temporal detector resolution in such
experiments is insufficient to resolve photon bunching peaks in
$g^{(2)}(\tau)$ due to small coherence times of thermal and only partially
coherent light. To address this issue, a generalized Siegert relation (GSR) is
assumed to hold in the partially coherent regime, which is characterized by
$g^{(2)}(0)<2$, yet not coherent emission with $g^{(2)}(0)=1$ Blazek  _et al._
(2011)

$\displaystyle g^{(2)}(\tau)=1+a|g^{(1)}(\tau)|^{2}.$ (2)

It asserts that intensity correlations captured by $g^{(2)}(\tau)$ are
directly related to first-order correlations and decay on a timescale given by
the coherence time

$\displaystyle\tau_{\mathrm{c}}=\int_{-\infty}^{\infty}|g^{(1)}(\tau)|^{2}\mathrm{d}\tau,$
(3)

which defines the timescale on which $g^{(1)}(\tau)$ vanishes, e.g. in weakly
interacting systems we have $|g^{(1)}(\tau)|=\exp(-|\tau|/\tau_{\mathrm{c}})$
Loudon (2000). The factor $a=g^{(2)}(0)-1$ accounts for the reduced value of
$g^{(2)}(0)$ Blazek  _et al._ (2011). Using this relation allows to
reconstruct the zero-delay two-photon correlation $g^{(2)}(0)$ from detector-
limited HBT measurements. In this letter we investigate under which conditions
the GSR is applicable and where it might fail. For this purpose we employ a
fully quantum-optical laser model and calculate the two-time coherence
functions microscopically, allowing us to explicitly verify the validity of
the GSR for different classes of nanolasers. Finally, we discuss the use of
the GSR in the analysis of semiconductor nanolaser emission and its
implications for the interpretation of experimental results.

We begin by deriving equations of motion for the two-time correlation
functions $g^{(1)}(\tau)$ and $g^{(2)}(\tau)$ directly from a microscopic
Hamiltonian. For small systems the full dynamics can be calculated with the
Lindblad-von Neumann equation. For many emitters this is not feasible, and we
use an approach based on the Heisenberg equation of motion for operators
together with the cluster expansion approximation Fricke (1996); Leymann,
Foerster, and Wiersig (2014); Gies  _et al._ (2007); Kira and Koch (2008). We
consider a system of $N$ individual two-level emitters that couple uniformly
to a single mode of the electromagnetic field. The light-matter interaction is
described in dipole approximation and the corresponding Hamilton operator is
given by

$\displaystyle H=\frac{\varepsilon}{2}\sum_{i}^{N}\sigma^{(i)}+\hbar\omega
b^{\dagger}b+\sum_{j}^{N}\left(gb^{\dagger}\sigma^{(j)}_{-}+g^{*}b\sigma^{(j)}_{+}\right),$
(4)

where $b$ and $b^{\dagger}$ are photon annihilation and creation operators,
and $\sigma^{(i)}$ denote spin-$\frac{1}{2}$ particle operators acting on the
Hilbert space of the two-level system $i$. The $\varepsilon$ denotes the
transition energy of a single two-level system, $\omega$ is the frequency of
the cavity mode, and $g$ is the light-matter coupling strength. Resonator
losses, relaxation, and the pump processes are considered using Lindblad terms
with constant rates $\kappa$, $\gamma$ and $P$, respectively. Similar models
have been used to describe nanolasers with many Kreinberg  _et al._ (2020),
few Moody  _et al._ (2018) and even single emitters in a microcavity Strauf
and Jahnke (2011); Reithmaier  _et al._ (2004); Gies  _et al._ (2017); Laussy,
del Valle, and Tejedor (2009).

Figure 1: Laser with many emitters subjected to fast dephasing (regime A):
$N=500$, $g=0.03\,\frac{1}{\mathrm{ps}}$,
$\kappa=0.04\,\frac{1}{\mathrm{ps}}$, $\Gamma=3\,\frac{1}{\mathrm{ps}}$ and
$\beta=0.1$. (a) Intracavity photon number and second-order correlation
$g^{(2)}(0)$ as function of pump power. The orange dashed line shows the
reconstructed $g^{(2)}(0)$ calculated according to Eq. (8). Colored markers
indicate pump rates at which $g^{(2)}(\tau)$ evaluated in panels below.
(b)-(d) Two-time values $g^{(2)}(\tau)$ at fixed pump rates. Dashed lines show
predictions by the GSR according to Eq. (2). The inset in (b) highlights small
differences between $g^{(2)}(\tau)$ and Eq. (2) due to residual emitter
correlations.

The equation-of-motion approach couples quantities that account for
correlations between particles up to arbitrary order in the interaction
strength. The arising hierarchy of coupled equations is truncated in the
cluster-expansion approximation Kira  _et al._ (1999) by neglecting
correlations of five or more particles, which gives rise to a closed set of
equations (see supplementary material). We use the quantum regression theorem
Gardiner and Zoller (2004) to derive equations of motion for two-time photon-
correlation functions using the cluster-expansion. From this we obtain a
relation between $g^{(1)}(\tau)$ and $g^{(2)}(\tau)$ that contains the Siegert
relation (1) and a correction term

$\displaystyle g^{(2)}(\tau)=1+|g^{(1)}(\tau)|^{2}+\delta g^{(2)}(\tau).$ (5)

Eq. (5) is a key finding of our work and twofold important. First, we can show
that it recovers the Siegert relation (1) in thermal equilibrium using Wick’s
theorem for the Matsubara Green function Bruus and Flensberg (2004). For a
non-interacting system in thermal equilibrium, all pure correlations involving
three or more one-particle operators vanish and, in this case, $\delta
g^{(2)}(\tau)=0$. Second, our ansatz enables us to numerically evaluate
photon-correlation functions to show that under certain conditions corrections
$\delta g^{(2)}(\tau)$ can become sizeable.

In the following, we consider three specific parameter regimes. The first
regime represents a class of nanolasers with many emitters subjected to fast
dephasing (regime A), such as quantum-well nanolasers as considered in
[Kreinberg  _et al._ , 2020]. Second, we investigate a regime with fewer
emitters, in which correlations between individual emitters can form due to
their interaction with a common cavity mode (regime B) Jahnke  _et al._
(2016); Leymann  _et al._ (2015). Finally, we consider the limiting regime of
cQED (regime C), in which very few emitters couple strongly to a single laser
mode as in [Reithmaier  _et al._ , 2004]. We use our quantum optical model to
analyze the behavior of $g^{(1)}(\tau)$ and $\delta g^{(2)}(\tau)$ in all
three regimes and subsequently test the validity of the GSR and its
applicability for the characterization of nanolasers.

For regime A, Fig. 1(a) displays input-output characteristics and the steady
state value of $g^{(2)}(0)$ as a function of pump power. It shows the
threshold jump in the photon number (gray) and the transition of the two-
photon correlations (black) from the thermal to the coherent regime Chow,
Jahnke, and Gies (2014); Jagsch  _et al._ (2018). Fig 1(b)–(d) show the
correlation function $g^{(2)}(\tau)$ (solid lines) calculated by solving the
equations-of-motion for the correlations in Eq. (5). For comparison, we plot
the result obtained from the GSR in Eq. (2) (dashed lines), which allows
comparing the time scales on which first- and second-order correlations decay.
We observe that in regime A the GSR gives a good approximation at all pump
rates, and in particular in the partially coherent regime for
$g^{(2)}(0)>1.3$. As we observe in the inset of Fig. 1(b), there are small
deviations from the GSR even deep in the thermal regime ($g^{(2)}(0)=2$). This
deviation is related to a residual influence of emitter correlations in the
system that are treated consistently in our quantum-optical model but are
largely suppressed by fast dephasing.

The effect of emitter correlations becomes more evident when the number of
emitters is reduced and emitter correlations are not strongly suppressed by
dephasing, which is the case in regime B. The presence of emitter correlations
gives rise to sub- and superradiant coupling, as observed e.g. in [Kreinberg
_et al._ , 2017] and [Leymann  _et al._ , 2015].

Figure 2: Laser representing regime B: $N=100$,
$g=0.008\,\frac{1}{\mathrm{ps}}$, $\kappa=0.1\,\frac{1}{\mathrm{ps}}$ and
$\beta=0.2$. (a) Intracavity photon number as function of pump power and
second-order correlation $g^{(2)}(0)$ showing superthermal emission below
threshold due to correlations. The black dashed line marks the value
$g^{(2)}(0)=2$. The orange dashed line shows the reconstructed $g^{(2)}(0)$
calculated according to Eq. (8). Colored markers indicate pump rates at which
$g^{(2)}(\tau)$ evaluated in panels below. (b)-(d) Two-time values
$g^{(2)}(\tau)$ at fixed pump rates. Dashed lines show predictions by the GSR
according to Eq. (2), which fails to reproduce the antibunching
$g^{(2)}(\tau)<1$ at finite delay $\tau$. The inset in (b) shows
$g^{(2)}(\tau)$ calculated without emitter correlations which restores
agreement with the GSR.

In Fig. 2 we show results obtained for a prototypical device operating in
regime B. It indeed exhibits superthermal bunching features at low excitation
(Fig. 2(a)), which have been identified as a signature in the subradiant
regime Leymann  _et al._ (2015). Here, the assumption of individual,
uncorrelated emitters is not justified, and consequently we observe deviations
from the GSR in panel (b)–(d) at all pump rates. In particular,
$g^{(2)}(\tau)$ shows anti-correlated behavior $g^{(2)}(\tau)<1$ for finite
delays $\tau$ that the GSR fails to capture since, by definition, we have
$|g^{(1)}(\tau)|^{2}>0$ for all $\tau$. Nevertheless, the GSR correctly
quantifies the timescale, on which the second-order correlations decay below
the threshold. In Fig. 2(c) we have $g^{(2)}(0)=2$, which is commonly
associated with thermal emission where the Siegert relation is supposed to
hold. Still, we see a significant deviation from the results of the full
calculation, which is an example of the fact that the zero-delay second-order
coherence is an indicator, but not a proof for a particular state of the
electromagnetic field.

To prove that the failure of the GSR is indeed caused by emitter correlations,
we take advantage of the possibility to suppress their effect in our
theoretical model. While emitter correlations are inseparable from the rest of
the dynamics in a Lindblad-von-Neumann equation, our equation-of-motion
approach has the feature that it allows to specifically neglect terms
accounting for emitter correlations. In the inset of Fig. 2(b) we show
$g^{(2)}(\tau)$ for this case. We observe that the superthermal bunching for
$\tau=0$ vanishes, and we recover the physics represented by the Siegert
relation. This highlights that the GSR should only be used for the
characterization of nanolaser devices when it is clear that emitter-
correlation effects related to sub- and superradiance are largely suppressed.

Finally, Fig. 3 shows results for regime C, which is a system of three
emitters coupled to a single cavity mode. In this cQED limit, lasing is
typically possible in the strong-coupling regime Gies  _et al._ (2017), and we
treat this system exactly by solving the Lindblad-von-Neumann equation, which
contains correlations up to arbitrary order in the interaction strength. The
parameters are chosen such that the light-matter interaction exceeds all
dissipation rates in the system, which is necessary to achieve lasing Gies
_et al._ (2017). Here we observe even stronger superthermal bunching below
threshold than in the previous case. The presence of strong light-matter
coupling gives rise to Rabi oscillations, which manifest themselves also in
the coherence properties (see Fig. 3(b)–(d)). Clearly, the validity of the GSR
is no longer given in this regime, and it is no longer justified to use it for
the estimation of the time scales, on which the second-order correlations
decay.

In the remaining part of this letter, we discuss the applicability of the GSR
in experimental characterization of semiconductor nanolasers where the second-
order photon-correlation $g^{(2)}(0)$ acts as a fingerprint of the light’s
photon statistics. It is determined in an HBT setup Jagsch  _et al._ (2018);
Hayenga  _et al._ (2016); Kreinberg  _et al._ (2020); Ulrich  _et al._ (2007),
in which the used single-photon detectors exhibit a timing jitter, i.e. an
uncertainty in photon arrival times, limiting the temporal resolution of the
measurement Korzh  _et al._ (2020). Together with the small coherence time of
the emission for thermal and partially coherent light, this prevents the
resolution of photon bunching at zero time delay in this regime. Fig. 4
illustrates the averaging procedure caused by the detector’s temporal
resolution using the example of a thermal $g^{(2)}(\tau)$ (blue line). The
timing jitter redistributes the recorded arrival times of photons, which is
modelled by a convolution of the correlation function $g^{(2)}(\tau)$

$\displaystyle\tilde{g}^{(2)}(\tau)=\int_{-\infty}^{\infty}g^{(2)}(t)F(\tau-t)\mathrm{d}t,$
(6)

with a detector response function (DRF) $F(t)$ shown in the inset of Fig. 4,
where the temporal resolution $\Delta t$ of the detector defines the FWHM of
the DRF. The result of the convolution corresponds to the measured
$\tilde{g}^{(2)}(\tau)$ (mauve line in Fig. 4). Clearly, the convoluted
$\tilde{g}^{(2)}(0)$ significantly differs from $g^{(2)}(\tau)$, showing
strongly reduced signatures of photon bunching $\tilde{g}^{(2)}(0)<2$ at
thermal emission. Crucially, the areas under the curves are invariant under
the convolution due to the normalization of the DRF.

Figure 3: Laser with few emitters operating in regime C: $N=3$,
$g=0.3\,\frac{1}{\mathrm{ps}}$, $\kappa=0.07\,\frac{1}{\mathrm{ps}}$,
$\beta=1$. (a) Intracavity photon number and second-order correlation
$g^{(2)}(0)$ as function of pump power. Colored markers indicate pump rates at
which $g^{(2)}(\tau)$ evaluated in panels below. (b)-(d) Two-time values
$g^{(2)}(\tau)$ at fixed pump rates show pronounced oscillations due to strong
coupling conditions. Dashed lines show predictions by the GSR according to Eq.
(2), which fail to reproduce the time scales of $g^{(2)}(\tau)$ in this
regime. Figure 4: Illustration of the convolution of $g^{(2)}(\tau)$. Here, a
thermal $g^{(2)}(\tau)$ ($g^{(2)}(0)=2$) was convoluted with a DRF modelled by
a Gaussian function. The FWHM of the DRF function is given by $\Delta
t=20\,\mathrm{ps}$ and quantifies the timing jitter of the detector.

By integrating $g^{(2)}(\tau)-1$ we find

$\displaystyle\tau_{\mathrm{area}}=\int_{-\infty}^{\infty}(g^{(2)}(\tau)-1)\mathrm{d}\tau=a\int_{-\infty}^{\infty}|g^{(1)}(\tau)|^{2}\mathrm{d}\tau=a\tau_{\mathrm{c}},$
(7)

where we used Eq. (2) and Eq. (3) for the coherence time $\tau_{\mathrm{c}}$.
Thus, assuming that the GSR is valid for all pump rates, it can be used to
reconstruct the zero-delay second-order coherence $g^{(2)}(0)$. We can
rearrange equation (7) and insert it into Eq. (2) to obtain the reconstructed
value of $g^{(2)}(0)$:

$\displaystyle
g^{(2)}_{\mathrm{rec.}}(0)=1+\frac{\tau_{\mathrm{area}}}{\tau_{\mathrm{c}}},$
(8)

i.e., we can determine $g^{(2)}_{\mathrm{rec.}}(0)$ from quantities that are
not subject to the temporal resolution of the detector. As the convolution
with the DRF does not change $\tau_{\mathrm{area}}$, we can integrate
$\tilde{g}^{(2)}(\tau)$ directly to obtain $\tau_{\mathrm{area}}$.

We illustrate the reconstruction process by adding the results for
$g^{(2)}_{\mathrm{rec.}}(0)$ to Figs. 1(a) and 2(a) as orange dashed line. As
we have discussed above, the GSR gives a good approximation to the real
$g^{(2)}(\tau)$ for the thermal and partially coherent regimes only if emitter
correlations are negligible. For regime A, where these conditions are
fulfilled the described procedure leads to good agreement with the actual
$g^{(2)}(0)$. Below $g^{(2)}(0)=1.3$ we see a deviation, indicating that the
GSR ceases to be valid beyond this point due to oscillations appearing in
$g^{(2)}(\tau)$. However, in a real experiment the photon bunching of the
thermal regime and most of the transition regime $g^{(2)}(0)<2$ can be
reliably reconstructed with this method. For regime B, the reconstructed
$g^{(2)}(0)$ severely underestimates the true value which is linked to the
anti-correlated behavior of $g^{(2)}(\tau)$ at finite $\tau$.

In conclusion, we have reviewed the validity of the GSR and its application in
the characterization of high-$\beta$ nanolasers with semiconductor based gain
media. We have developed a quantum-optical approach to access two-time second-
order correlations functions taking into account emitter correlations
consistently, giving rise to a correction to the Siegert relation that makes
it valid in all emission regimes. In nanolasers with extended inhomogeneous
gain media, such as semiconductor quantum wells or ensembles of quantum dots,
we expect that fast dephasing suppresses emitter correlations, such that the
GSR is largely valid and can be used for reconstructing
$g^{(2)}_{\mathrm{rec.}}(0)$. In the presence of emitter correlations, strong
deviations from the exact calculations occur, leading us to conclude that it
is imperative to verify the absence of such correlations prior to invoking the
GSR. For future studies our equation of motion technique can be extended to
access $g^{(2)}(\tau)$ for a range of different nanolaser devices, such that
relying on the GSR becomes unnecessary.
See the supplementary material for the details regarding the quantum
mechanical treatment of the Siegert relation and the derivation of the
quantum-optical laser model.
F.L. acknowledges funding by the central research development fund (CRDF) of
the University of Bremen.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

## References

  * Loudon (2000) R. Loudon,  _The quantum theory of light_ , 3rd ed. (Oxford University Press, Oxford ; New York, 2000).
  * Glauber (1963) R. J. Glauber, “The quantum theory of optical coherence,” Physical Review 130, 2529–2539 (1963).
  * Wang, Puccioni, and Lippi (2015) T. Wang, G. P. Puccioni, and G. L. Lippi, “Dynamical Buildup of Lasing in Mesoscale Devices,” Scientific Reports 5, 15858 (2015).
  * Roy-Choudhury, Haas, and Levi (2009) K. Roy-Choudhury, S. Haas, and A. F. J. Levi, “Quantum fluctuations in small lasers,” Phys. Rev. Lett. 102, 053902 (2009).
  * Schlottmann  _et al._ (2018) E. Schlottmann, M. von Helversen, H. A. M. Leymann, T. Lettau, F. Krüger, M. Schmidt, C. Schneider, M. Kamp, S. Höfling, J. Beyer, J. Wiersig, and S. Reitzenstein, “Exploring the photon-number distribution of bimodal microlasers with a transition edge sensor,” Phys. Rev. Applied 9, 064030 (2018).
  * Hamsen  _et al._ (2017) C. Hamsen, K. N. Tolazzi, T. Wilk, and G. Rempe, “Two-photon blockade in an atom-driven cavity QED system,” Physical Review Letters 118, 133604 (2017).
  * Snijders  _et al._ (2018) H. Snijders, J. Frey, J. Norman, H. Flayac, V. Savona, A. Gossard, J. Bowers, M. van Exter, D. Bouwmeester, and W. Löffler, “Observation of the unconventional photon blockade,” Physical Review Letters 121, 043601 (2018).
  * Paik  _et al._ (2019) E. Y. Paik, L. Zhang, G. W. Burg, R. Gogna, E. Tutuc, and H. Deng, “Interlayer exciton laser of extended spatial coherence in atomically thin heterostructures,” Nature 576, 80–84 (2019).
  * Megyeri  _et al._ (2018) B. Megyeri, G. Harvie, A. Lampis, and J. Goldwin, “Directional bistability and nonreciprocal lasing with cold atoms in a ring cavity,” Physical Review Letters 121, 163603 (2018).
  * Shan  _et al._ (2021) H. Shan, L. Lackner, B. Han, E. Sedov, C. Rupprecht, H. Knopf, F. Eilenberger, K. Yumigeta, K. Watanabe, T. Taniguchi, S. Klembt, S. Höfling, A. V. Kavokin, S. Tongay, C. Schneider, and C. Antón-Solanas, “Coherent light emission of exciton-polaritons in an atomically thin crystal at room temperature,” arXiv:2103.10459 (2021).
  * Blazek  _et al._ (2011) M. Blazek, S. Hartmann, A. Molitor, and W. Elsaesser, “Unifying intensity noise and second-order coherence properties of amplified spontaneous emission sources,” Optics Letters 36, 3455 (2011).
  * Wang  _et al._ (2021) T. Wang, J. Zou, G. P. Puccioni, W. Zhao, X. Lin, H. Chen, G. Wang, and G. L. Lippi, “Methodological investigation into the noise influence on nanolasers’ large signal modulation,” Optics Express 29, 5081 (2021).
  * Jagsch  _et al._ (2018) S. T. Jagsch, N. V. Triviño, F. Lohof, G. Callsen, S. Kalinowski, I. M. Rousseau, R. Barzel, J.-F. Carlin, F. Jahnke, R. Butté, C. Gies, A. Hoffmann, N. Grandjean, and S. Reitzenstein, “A quantum optical study of thresholdless lasing features in high-$\beta$ nitride nanobeam cavities,” Nature Communications 9, 564 (2018).
  * Ma and Oulton (2019) R.-M. Ma and R. F. Oulton, “Applications of nanolasers,” Nature Nanotechnology 14, 12–22 (2019).
  * Kreinberg  _et al._ (2020) S. Kreinberg, K. Laiho, F. Lohof, W. E. Hayenga, P. Holewa, C. Gies, M. Khajavikhan, and S. Reitzenstein, “Thresholdless tansition to coherent emission at telecom wavelengths from coaxial nanolasers with excitation power dependent $\beta$-factors,” Laser & Photonics Reviews 14, 2000065 (2020).
  * Ulrich  _et al._ (2007) S. M. Ulrich, C. Gies, S. Ates, J. Wiersig, S. Reitzenstein, C. Hofmann, A. Löffler, A. Forchel, F. Jahnke, and P. Michler, “Photon statistics of semiconductors microcavity lasers,” Physical Review Letters 98, 043906 (2007).
  * Takemura  _et al._ (2019) N. Takemura, M. Takiguchi, E. Kuramochi, A. Shinya, T. Sato, K. Takeda, S. Matsuo, and M. Notomi, “Lasing thresholds and photon statistics in high-$\beta$ buried multiple quantum well photonic crystal nanocavity lasers,” Physical Review A 99, 053820 (2019).
  * Mork and Lippi (2018) J. Mork and G. L. Lippi, “Rate equation description of quantum noise in nanolasers with few emitters,” Applied Physics Letters 112, 141103 (2018).
  * Siegert (1943) A. J. F. Siegert,  _On the fluctuations in signals returned by many independently moving scatterers_ (Radiation Laboratory, Massachusetts Institute of Technology, Cambridge, Mass., 1943).
  * Voigt and Hess (1994) H. Voigt and S. Hess, “Comparison of the intensity correlation function and the intermediate scattering function of fluids: a molecular dynamics study of the Siegert relation,” Physica A: Statistical Mechanics and its Applications 202, 145–164 (1994).
  * Ferreira  _et al._ (2020) D. Ferreira, R. Bachelard, W. Guerin, R. Kaiser, and M. Fouché, “Connecting field and intensity correlations: The Siegert relation and how to test it,” American Journal of Physics 88, 831–837 (2020).
  * Fricke (1996) J. Fricke, “Transport equations including many-particle correlations for an arbitrary quantum system: A general formalism,” Annals of Physics 252, 479–498 (1996).
  * Foellmi (2009) C. Foellmi, “Intensity interferometry and the second-order correlation function $g^{(2)}$ in astrophysics,” Astronomy & Astrophysics 507, 1719–1727 (2009).
  * Borycki, Kholiqov, and Srinivasan (2016) D. Borycki, O. Kholiqov, and V. J. Srinivasan, “Interferometric near-infrared spectroscopy directly quantifies optical field dynamics in turbid media,” Optica 3, 1471–1476 (2016).
  * Stetefeld, McKenna, and Patel (2016) J. Stetefeld, S. A. McKenna, and T. R. Patel, “Dynamic light scattering: a practical guide and applications in biomedical sciences,” Biophysical Reviews 8, 409–427 (2016).
  * Aßmann  _et al._ (2010) M. Aßmann, F. Veit, J.-S. Tempel, T. Berstermann, H. Stolz, M. van der Poel, J. M. Hvam, and M. Bayer, “Measuring the dynamics of second-order photon correlation functions inside a pulse with picosecond time resolution,” Optics Express 18, 20229 (2010).
  * Chow, Jahnke, and Gies (2014) W. W. Chow, F. Jahnke, and C. Gies, “Emission properties of nanolasers during the transition to lasing,” Light: Science & Applications 3, e201 (2014).
  * Kreinberg  _et al._ (2017) S. Kreinberg, W. W. Chow, J. Wolters, C. Schneider, C. Gies, F. Jahnke, S. Höfling, M. Kamp, and S. Reitzenstein, “Emission from quantum-dot high-$\beta$ microcavities: transition from spontaneous emission to lasing and the effects of superradiant emitter coupling,” Light: Science & Applications 6, e17030 (2017).
  * Elvira  _et al._ (2011) D. Elvira, X. Hachair, V. B. Verma, R. Braive, G. Beaudoin, I. Robert-Philip, I. Sagnes, B. Baek, S. W. Nam, E. A. Dauler, I. Abram, M. J. Stevens, and A. Beveratos, “Higher-order photon correlations in pulsed photonic crystal nanolasers,” Physical Review A 84, 061802 (2011).
  * Yokoo  _et al._ (2017) A. Yokoo, M. Takiguchi, M. D. Birowosuto, K. Tateno, G. Zhang, E. Kuramochi, A. Shinya, H. Taniyama, and M. Notomi, “Subwavelength nanowire lasers on a silicon photonic crystal operating at telecom wavelengths,” ACS Photonics 4, 355–362 (2017).
  * Leymann, Foerster, and Wiersig (2014) H. A. M. Leymann, A. Foerster, and J. Wiersig, “Expectation value based equation-of-motion approach for open quantum systems: A general formalism,” Physical Review B 89, 085308 (2014).
  * Gies  _et al._ (2007) C. Gies, J. Wiersig, M. Lorke, and F. Jahnke, “Semiconductor model for quantum-dot-based microcavity lasers,” Physical Review A 75, 013803 (2007).
  * Kira and Koch (2008) M. Kira and S. W. Koch, “Cluster-expansion representation in quantum optics,” Physical Review A 78, 022102 (2008).
  * Moody  _et al._ (2018) G. Moody, M. Segnon, I. Sagnes, R. Braive, A. Beveratos, I. Robert-Philip, N. Belabas, F. Jahnke, K. L. Silverman, R. P. Mirin, M. J. Stevens, and C. Gies, “Delayed formation of coherence in the emission dynamics of high-Q nanolasers,” Optica 5, 395 (2018).
  * Strauf and Jahnke (2011) S. Strauf and F. Jahnke, “Single quantum dot nanolaser,” Laser & Photonics Reviews 5, 607–633 (2011).
  * Reithmaier  _et al._ (2004) J. P. Reithmaier, G. Sęk, A. Löffler, C. Hofmann, S. Kuhn, S. Reitzenstein, L. V. Keldysh, V. D. Kulakovskii, T. L. Reinecke, and A. Forchel, “Strong coupling in a single quantum dot-semiconductor microcavity system,” Nature 432, 197–200 (2004).
  * Gies  _et al._ (2017) C. Gies, F. Gericke, P. Gartner, S. Holzinger, C. Hopfmann, T. Heindel, J. Wolters, C. Schneider, M. Florian, F. Jahnke, S. Höfling, M. Kamp, and S. Reitzenstein, “Strong light-matter coupling in the presence of lasing,” Phys. Rev. A 96, 023806 (2017).
  * Laussy, del Valle, and Tejedor (2009) F. P. Laussy, E. del Valle, and C. Tejedor, “Luminescence spectra of quantum dots in microcavities. I. Bosons,” Physical Review B 79, 235325 (2009).
  * Kira  _et al._ (1999) M. Kira, F. Jahnke, W. Hoyer, and S. W. Koch, “Quantum theory of spontaneous emission and coherent effects in semiconductor microstructures,” Progress in Quantum Electronics 23, 189–279 (1999).
  * Gardiner and Zoller (2004) C. Gardiner and P. Zoller,  _Quantum Noise_ , 3rd ed. (Springer Berlin, Heidelberg, 2004).
  * Bruus and Flensberg (2004) H. Bruus and K. Flensberg,  _Many-Body Quantum Theory in Condensed Matter Physics: An Introduction_ (OUP Oxford, 2004).
  * Jahnke  _et al._ (2016) F. Jahnke, C. Gies, M. Aßmann, M. Bayer, H. A. M. Leymann, A. Foerster, J. Wiersig, C. Schneider, M. Kamp, and S. Höfling, “Giant photon bunching, superradiant pulse emission and excitation trapping in quantum-dot nanolasers,” Nature Communications 7, 11540 (2016).
  * Leymann  _et al._ (2015) H. A. M. Leymann, A. Foerster, F. Jahnke, J. Wiersig, and C. Gies, “Sub- and superradiance in nanolasers,” Physical Review Applied 4, 044018 (2015).
  * Hayenga  _et al._ (2016) W. E. Hayenga, H. Garcia-Gracia, H. Hodaei, C. Reimer, R. Morandotti, P. LiKamWa, and M. Khajavikhan, “Second-order coherence properties of metallic nanolasers,” Optica 3, 1187–1193 (2016).
  * Korzh  _et al._ (2020) B. Korzh, Q.-Y. Zhao, J. P. Allmaras, S. Frasca, T. M. Autry, E. A. Bersin, A. D. Beyer, R. M. Briggs, B. Bumble, M. Colangelo, G. M. Crouch, A. E. Dane, T. Gerrits, A. E. Lita, F. Marsili, G. Moody, C. Peña, E. Ramirez, J. D. Rezac, N. Sinclair, M. J. Stevens, A. E. Velasco, V. B. Verma, E. E. Wollman, S. Xie, D. Zhu, P. D. Hale, M. Spiropulu, K. L. Silverman, R. P. Mirin, S. W. Nam, A. G. Kozorezov, M. D. Shaw, and K. K. Berggren, “Demonstration of sub-3 ps temporal resolution with a superconducting nanowire single-photon detector,” Nature Photonics 14, 250–255 (2020).