# Low probability states, data statistics, and entropy estimation

Damián G. Hernández 1,4,∗ Ahmed Roman 1,∗ Ilya Nemenman1,2,3 1Physics
Department, Emory University, Atlanta, Georgia, USA
2 Biology Department, Emory University, Atlanta, Georgia, USA
3 Initiative for Theory and Modeling of Living Systems, Emory University,
Atlanta, Georgia, USA
4 Department of Medical Physics, Centro Atómico Bariloche and Instituto
Balseiro, 8400 San Carlos de Bariloche, Argentina
∗ D.G.H. and A.R. contributed equally to this work.

(August 28, 2024)

###### Abstract

A fundamental problem in analysis of complex systems is getting a reliable
estimate of entropy of their probability distributions over the state space.
This is difficult because unsampled states can contribute substantially to the
entropy, while they do not contribute to the Maximum Likelihood estimator of
entropy, which replaces probabilities by the observed frequencies. Bayesian
estimators overcome this obstacle by introducing a model of the low-
probability tail of the probability distribution. Which statistical features
of the observed data determine the model of the tail, and hence the output of
such estimators, remains unclear. Here we show that well-known entropy
estimators for probability distributions on discrete state spaces model the
structure of the low probability tail based largely on few statistics of the
data: the sample size, the Maximum Likelihood estimate, the number of
coincidences among the samples, the dispersion of the coincidences. We derive
approximate analytical entropy estimators for undersampled distributions based
on these statistics, and we use the results to propose an intuitive
understanding of how the Bayesian entropy estimators work.

entropy estimation; Pitman-Yor Mixture estimator; Nemenman-Shafee-Bialek
estimator; Bayesian inference; coincidences; tail probability.

††preprint: AAPM/123-QED

##  I Introduction

Estimating entropy – that is, the measure of uncertainty [1, 2] – of a random
variable from its samples is often a key question in analysis of complex
systems. This estimation from a finite (and often small) set of samples is a
hard problem, especially for high dimensional systems, where the number of
states that a variable can take quickly overwhelms the number of samples $N$.
Then many of the states, hereafter called low probability states, have
probability $<1/N$. Collectively, we refer to all of these states as the tail
of the probability distribution. While there may be a lot of samples in the
tail, each low probability state will not be sampled typically, or will be
sampled at most once. Because of the tail, the entropy estimator that replaces
probabilities of states by their empirical frequencies (the so called naive or
Maximum Likelihood estimator [3]) has a large sample size dependent bias [4].
Corrections have been derived to overcome this bias [5, 6, 7], but these tend
to be valid only in the well-sampled regime. Outside of this regime, Bayesian
[8, 9, 10] and some non-parametric [11, 12, 13] estimators may still result in
low bias estimates by imposing a priori assumptions on the probabilities of
the low-probability states.

Although these Bayesian and non-parametric estimators perform well on some
data sets, it is known that no estimator can be universally unbiased in this
regime [4, 14]. Thus it is crucial to understand how these estimators extract
information about entropy from data, and hence when they will fail.
Unfortunately, such theoretical understanding is missing for many estimators.
Ma was the first to point out that estimation of entropy is possible for
poorly-sampled uniform distributions by analysing a particular statistics of
the data: coincidences [15]. Nemenman extended the theoretical idea that
coincidences determine entropy to non-uniform distributions obeying some
Bayesian priors [16]. However, a similar theoretical understanding is still
missing in a broader context, and it remains unclear which statistics of data,
in addition to the number of coincidences, may contribute to entropy
estimation and why.

In this paper, we analytically investigate two Bayesian estimators: that of
Nemenman, Shafee and Bialek [9, 17] and of Archer and Pillow [10]. We focus on
the regime, which is arguably the most important for real life applications,
where the number of states with at least one sample, $K_{1}$, is similar to
the total number of samples, $K_{1}\sim N\gg 1$, and yet $K_{1}<N$, so that
there are coincidences in the data. Outside of this regime, the probability
distribution is either well-sampled (so that many different methods for
entropy estimation would work), or there are no coincidences at all (so that
entropy estimation is impossible). In our regime of interest, we show that the
result of the estimation by the studied estimators depends on the Maximum
Likelihood entropy estimate $S_{0}$, the number of coincidences, and also on
two measures of dispersion of coincidences. The first of these, $K_{2}$, is
the number of states with at least two samples. The second, which we call
$Q_{1}$, characterizes the spread of coincidences over states with three or
more samples.

We show that values of these statistics are related to the structure of the
tails of the probability distribution that is assumed by the estimators.
Specifically, a short, exponential, tail is more likely to be inferred by the
estimators when there many coincidences or they are dispersed. If the number
of coincidences is intermediate, and the coincidences are concentrated, then
the estimators infer a long tail. In between these two regions, a mixed tail
dominates. We show that the studied estimators correct Maximum Likelihood, and
that the correction is larger when there are fewer coincidences and they are
concentrated, which in turn happens with a large exponential tail or a slowly-
decaying long tail. This understanding relates the observable data statistics
to assumptions that Bayesian estimators make about the underlying probability
distributions (see Fig. 1), and hence provides an intuitive explanation for
how these estimators work and, crucially, when they fail.

##  II Overview of Bayesian entropy estimation

Given a probability distribution $\\{q_{x}\\}=\bm{q}$ for a discrete one-
dimensional random variable $X$, its entropy is defined as [1]

$S(\bm{q})=-\sum_{x}q_{x}\log q_{x}.$ (1)

Note that we use the natural logarithm throughout this paper, and hence
entropy is measured in nats. One is often faced with a problem when $S$ must
be estimated for unknown $q_{x}$ from a set of $N$ samples
$\\{x_{1},\dots,x_{N}\\}$ from the probability distribution. The Maximum
Likelihood estimator of entropy, $S_{0}$, is then defined by replacing the
probabilities with frequencies $q_{x}\to\hat{q}_{x}=n_{x}/N$,

$S_{0}=S(\hat{\bm{q}})=-\sum_{x}\frac{n_{x}}{N}\log\frac{n_{x}}{N}.$ (2)

States with zero frequencies in the sample do not contribute to $S_{0}$
resulting typically in underestimation of the entropy [4]. In general, because
of this low probability tail, estimation of entropy from data is very hard
when the number of samples is smaller than the number of effective states of
the variable, $N\ll\exp(S)$.

Figure 1: Relation between assumptions about the tail structure and the
statistics that determine entropy estimation. The set of unsampled states,
$q_{i}\leq 1/N$, which we refer to as the tail, may contribute substantially
to the entropy. However, the Maximum Likelihood estimation overlooks this
contribution. If the rank ordered plot of the tail is exponential with the
scale $\alpha$ (top panel), then the tail has effectively $\alpha$ states,
which contribute $\delta S\sim\log\alpha$ to the entropy. While the tail
cannot be observed directly, it pulls samples from the head of the
distribution, so that the number of coincidences, $\Delta$, in the head
decreases as $\alpha$ grows. Thus one can estimate $\alpha$ and hence the
entropy itself from $\Delta$. Alternatively, if the rank-ordered plot of the
tail has a power law structure with the exponent $-1/d$, then the tail does
not have a finite effective size (bottom panels). Then its contribution to
entropy depends on $d$ as $\delta S\sim(1-d)^{-1}$. In this case, one can
estimate $d$, and hence the entropy, from the dispersion of the coincidences,
which depends, in part, on how many samples happen once or more, $K_{1}$, or
twice or more, $K_{2}$, in the dataset.

Bayesian estimators address the problem by imposing various a priori
assumptions $p(\bm{q})$. One then uses Bayes theorem to infer the a posteriori
distribution of $\bm{q}$, and finally integrates over $\bm{q}$ to get the a
posteriori distribution or moments of entropy. Specifically, the mean
posterior entropy $\hat{S}=\langle S|\bm{n}\rangle$ given the counts
$\bm{n}=\\{n_{x}\\}$ of how many times state $x$ was sampled is given by

$\displaystyle\hat{S}$ $\displaystyle=\langle S|\bm{n}\rangle=\int
S(\bm{q})p(S|\bm{q})p(\bm{q}|\bm{n})d\bm{q}$ $\displaystyle=\int
S(\bm{q})\delta\left(S+\sum_{x}q_{x}\log q_{x}\right)p(\bm{q}|\bm{n})d\bm{q},$
(3)

where $p(\bm{q}|\bm{n})$ is the posterior over $\bm{q}$ under some prior
$p(\bm{q})$,

$p(\bm{q}|\bm{n})=\frac{p(\bm{n}|\bm{q})p(\bm{q})}{p(\bm{n})}=\frac{\prod_{x}q_{x}^{n_{x}}p(\bm{q})}{p(\bm{n})}.$
(4)

For distributions with known finite size $\mathcal{A}$ of the space of the
possible outcomes (aka the alphabet size), the Dirichlet distribution is often
chosen as a prior due to its conjugacy with the categorical distribution:

$p(\bm{q})=\text{Dirichlet}(\bm{q}|\lambda)\propto\prod_{i=1}^{\mathcal{A}}q_{i}^{\lambda},$
(5)

where $\lambda$ is known as the concentration parameter.

Note that any chosen prior $p(\bm{q})$ implicitly imposes assumptions on the
structure of the low probability tail (and hence its contribution to the
entropy) based on the observed statistics of the well-sampled part of the
probability distribution. However, these implicit assumptions usually are not
made explicit, and they remain mysterious even for most commonly used Bayesian
estimators. Lifting this veil is the goal of this work.

###  II.1 The Nemenman-Shafee-Bialek (NSB) Estimator

Nemenman et al. [9] showed that, for variables with the finite alphabet size
$\mathcal{A}$, Dirichlet priors on $\bm{q}$ with a fixed value for the
concentration parameter $\lambda$ correspond to highly concentrated a priori
distribution on entropy, which persists for large sample sizes. This bias
induces incorrect entropy estimates, which nonetheless have low variance and
hence are certain about their outputs. To address this issue, Ref. [9]
suggested a Dirichlet-mixture prior

$p_{\rm NSB}(\bm{q})=\int\text{Dirichlet}(q|\lambda)p_{\rm
prior}(\lambda)d\lambda,$ (6)

where $p(\lambda)$ are the mixture weights determined by

$p_{\rm prior}(\lambda)\propto\partial_{\lambda}\langle
S|\lambda\rangle=\mathcal{A}\psi_{1}(\mathcal{A}\lambda+1)-\psi_{1}(\lambda+1),$
(7)

and where $\langle S|\lambda\rangle$ is the a priori expected entropy under
the $\text{Dirichlet}(\bm{q}|\lambda)$ prior, and $\psi_{1}(\cdot)$ is the
tri-gamma function [18]. This choice of weights implies a nearly uniform a
priori distribution for the entropy $S$ on the interval $[0,\log\mathcal{A}].$
The resulting entropy estimate is then

$\displaystyle\hat{S}_{\text{NSB}}$ $\displaystyle=\langle
S|\bm{n}\rangle=\int\int
S(\bm{q})p(\bm{q}|\bm{n},\lambda)p(\lambda|\bm{n})d\bm{q}d\lambda$
$\displaystyle=\int\langle
S|\bm{n},\lambda\rangle\frac{p(\bm{n}|\lambda)p_{\rm
prior}(\lambda)}{p(\bm{n})}d\lambda.$ (8)

Here $\langle S|\bm{n},\lambda\rangle$ is the posterior mean entropy under the
prior $\text{Dirichlet}(\bm{q}|\lambda)$, and $p(\bm{n}|\lambda)$ is the
evidence (which has a Polya distribution) [19],

$\displaystyle p(\bm{n}|\lambda)$ $\displaystyle=\int
p(\bm{n}|\bm{q})p(\bm{q}|\lambda)d\bm{q}$
$\displaystyle=\frac{N!\Gamma(\mathcal{A}\lambda)}{\Gamma(\lambda)^{\mathcal{A}}\Gamma(N+\mathcal{A}\lambda)}\prod_{i=1}^{\mathcal{A}}\frac{\Gamma(n_{i}+\lambda)}{n_{i}!}$
(9)

where $\Gamma(\cdot)$ is the gamma function [18]. Using the analytical
expressions for the first two moments of posterior mean entropy $\langle
S|\bm{n},\lambda\rangle$ (available from Refs. [8, 9]), one then uses one-
dimensional numerical integration over $\lambda$ to obtain
$\hat{S}_{\text{NSB}}$.

###  II.2 The Dirichlet and the Pitman-Yor Processes

When the size of the state space is unknown or infinite, the standard NSB
construction does not work. Then one commonly uses one of the following two
stochastic processes to construct a prior $p(\bm{q})$ over a countably
infinite state space: the Pitman-Yor Process (PYP) [20] and its special case,
the Dirichlet Process (DP) [21]. To specify these processes, one requires two
inputs: a parameter vector and a base distribution. Parameters of the Pitman-
Yor process are known as the discount parameter $d$, $0\leq d<1$, and the
concentration parameter $\alpha$. The parameters control the shape of typical
distributions generated by the process. Specifically, $d$ controls the
structure of the low probability tail of $\bm{q}$, so that the tail typically
decays as $q_{x}\propto x^{-1/d}$. The concentration parameter $\alpha$
control the probability mass near the head of the distribution. In the limit
$d\rightarrow 0$, $\text{PYP}(d,\alpha)$ becomes the Dirichlet Process,
$\text{DP}(\alpha)$. In other words, the Dirichlet Process generates
distributions with short tails.

When the base distribution is the Beta distribution, one draws samples
$q_{x}\sim\text{PYP}(d,\alpha)$ via the so called stick-breaking process [22],
which uses an infinite sequence of independent Beta-distributed random
variables $\beta_{x}\sim\text{Beta}(1-d,\alpha+xd)$, so that

$\tilde{q}_{x}=\beta_{x}\prod_{y=1}^{x-1}(1-\beta_{y}).$ (10)

Thus obtained $\tilde{\bm{q}}$ are not strictly decreasing with $x$, and so
one obtains a strictly non-increasing distribution $\bm{q}$ from them by rank
ordering.

###  II.3 Expectations over DP and PYP Posteriors

Previous studies [23] showed that PYP priors (for multinomial observations)
yield a posterior $p(\bm{q}|\bm{n},\alpha,d)$, which consists of two parts:
probability of $K_{1}$ states that exist in the sample with the counts of, at
least, one, and probability of states that are not sampled. We will denote the
set of states with nonzero counts as $\mathbb{K}$, and its cardinality is
$K_{1}=||\mathbb{K}||$. Then the first term of the posterior is given by the
Dirichlet distribution,
$p(\bm{q}\in\mathbb{K}|\bm{\mu})\propto\prod_{x}q_{x}^{\mu_{x}}$, where
$\bm{\mu}$ is a concentration vector
$\bm{\mu}=(n_{1}-d,\cdots,n_{K_{1}}-d,\alpha+K_{1}d)$. This leaves the
probability of $q_{*}=1-\sum_{x\in\mathbb{K}}q_{x}$ for the unobserved states.
In other words, the states with nonzero counts contribute the following to the
posterior:

$\displaystyle p(\bm{q}\in\mathbb{K}|\bm{n})$
$\displaystyle=p(q_{1},\cdots,q_{K_{1}},q_{*}|\bm{n})$
$\displaystyle=\text{Dirichlet}(n_{1}-d,\cdots,n_{K_{1}}-d,\alpha+K_{1}d)$
$\displaystyle\propto
q_{*}^{\alpha+K_{1}d}\prod_{i=1}^{K_{1}}q_{i}^{n_{i}-d}.$ (11)

For the states that have no samples, the posterior is equal to the prior. Thus
their contribution to the posterior is the Pitman-Yor Process, normalized by
their total probability being $q^{*}$:

$p(\bm{q}\not\in\mathbb{K})=p(q_{K_{1}+1},q_{K_{1}+2},\cdots)=q^{*}\text{PYP}(d,\alpha+K_{1}d).$
(12)

Overall, this yields a closed form solution for the posterior mean and
variance of the entropy $S$. Specifically, the resulting posterior mean
$\langle S|\bm{n},\alpha,d\rangle$ is

$\begin{split}\langle
S|\bm{n},\alpha,d\rangle=\psi(\alpha+N+1)-\frac{\alpha+K_{1}d}{\alpha+N}\psi(1-d)\\\
-\frac{1}{\alpha+N}\left(\sum_{x=1}^{K_{1}}(n_{x}-d)\psi(n_{x}-d+1)\right),\end{split}$
(13)

where $\psi(x)=\partial_{x}\log\Gamma(x)$ is the di-gamma function [18].
Unfortunately, this is usually not a good estimate of entropy since, for fixed
$\alpha$ and $d$, the prior $\text{PYP}(d,\alpha)$ on $\bm{q}$ corresponds to
a highly concentrated a priori distribution on entropy, just like was noted
before in the context of the NSB estimator. To counter this, Archer and Pillow
[10] followed the NSB prescription and introduced a prior (mixture) over the
parameters of $PYP(d,\alpha)$, $p_{\rm prior}(\alpha,d)$, which uniformized
the induced prior over entropy (with the caveat that, for a distribution on a
countable alphabet, the entropy may be infinite, and hence strict uniform
distribution over entropy is impossible). Specifically, they used

$\displaystyle p_{\rm
prior}(\alpha,d)=p(\gamma)=e^{-10/(1-\gamma)},\quad\quad{\rm where}$ (14)
$\displaystyle\gamma=(\psi(1)-\psi(1-d))/(\psi(\alpha+1)-\psi(1-d)),$ (15)

and then they confirmed numerically that this choice of the prior leads to
good estimates of entropy for various test data sets. In other words, they
proposed a new estimate of entropy, the Pitman-Yor Mixture (PYM):

$\displaystyle\hat{S}_{PYM}$ $\displaystyle=\langle
S|\bm{n}\rangle=\int\langle S|\bm{n},\alpha,d\rangle p_{\rm
posterior}(\alpha,d|\bm{n})d(\alpha,d)$ $\displaystyle=\int\langle
S|\bm{n},\alpha,d\rangle\frac{p(\bm{n}|\alpha,d)p_{\rm
prior}(\alpha,d)}{p(\bm{n})}d(\alpha,d),$ (16)

where $\langle S|\bm{n},\alpha,d\rangle$ is given in Eq. (13). The evidence
$p(\bm{n}|\alpha,d)$ is then given by (see Ref. [10] for a detailed
derivation)

$p(\bm{n}|\alpha,d)=\frac{\Gamma(1+\alpha)\prod_{l=1}^{K_{1}}(\alpha+ld)\prod_{x=1}^{K_{1}}\Gamma(n_{x}-d)}{\Gamma(1-d)^{K_{1}}\Gamma(\alpha+N)}.$
(17)

Note that taking $d\to 0$ in Eqs. (16 and 17) and making the identification
$\alpha=\mathcal{A}\lambda$ in the limits $\lambda\rightarrow 0$ and
$\mathcal{A}\rightarrow\infty$ such that $\alpha$ is finite, result in a
countably-infinite analogue of the NSB estimator.

##  III Determining data statistics that define entropy estimates

In the section, we approximate the likelihood function of the Pitman-Yor
process, Eq. (17), analytically in terms of coincidence-based data statistics.
We then numerically show that the resulting analytical entropy estimates are
close to the exact Pitman-Yor Mixture estimator. We focus on the regime where
the Maximum Likelihood entropy estimator fails dramatically. For this, we
study random variables with many accessible states in the regime where the
number of unique samples, $K_{1}$, is of the order of the total sample size
$N$. This regime corresponds to $K_{1}\lesssim N\leq\exp(S)$, where $N$ is the
number of samples and $S$ is the true entropy.

We start by considering the log-likelihood function, which is the logarithm of
the evidence $p(\bm{n}|\alpha,d)$ in Eq. (17):

$\mathcal{L}({\bf
n}|\alpha,d)=\log\Gamma(1+\alpha)-\log\Gamma(N+\alpha)+\log\Gamma\left(\frac{\alpha}{d}+K_{1}\right)\\\
-\log\Gamma\left(\frac{\alpha}{d}+1\right)+\sum_{i=1}^{K_{1}}\log\Gamma(n_{i}-d)-K_{1}\log\Gamma(1-d).$
(18)

We now define $K_{m}$ as the number of states with at least $m$ counts in the
total sample of size $N$, $K_{m}=\sum_{n_{i}\geq m}1$. We denote by $m_{f}$
the largest occupancy of any state in the sample. Further, we define
$\mathcal{K}$ as the vector, whose $m$th element is $K_{m}$. We note that, for
any function $f(n)$,

$\sum_{i}f(n_{i})=\sum_{m}(K_{m}-K_{m+1})f(m).$ (19)

Thus, in particular, the log-likelihood $\mathcal{L}({\bf n}|\alpha,d)$ can be
viewed as $\mathcal{L}({\mathcal{K}}|\alpha,d)$. With this, we can expand Eq.
(18) around $d=0$ to get (see Appendix VI.1 for details):

$\mathcal{L}(\bm{n}|\alpha,d)\approx\mathcal{L}_{a}(\mathcal{K}|\alpha,d)\equiv\log\Gamma(1+\alpha)\\\
-\log\Gamma(N+\alpha)+\log\Gamma\left(\frac{\alpha}{d}+K_{1}\right)-\log\Gamma\left(\frac{\alpha}{d}+1\right)\\\
+(K_{1}-1)\log d+K_{2}\log(1-d)-Q_{1}d+\mathcal{O}(d^{2}),$ (20)

where

$Q_{1}=\sum_{m=3}^{m_{f}}\frac{K_{m}}{m-1},$ (21)

and the subscript $a$ denotes the $d\to 0$ asymptotic nature of the
expression.

Figure 2: Comparison between PYM and related estimators and their
approximations for distributions with different tails. The upper panels (a-c)
show the distributions, whose entropy is being estimated. The lower panels
(d-f) show the corresponding entropy estimates as a function of the number of
samples, averaged over ten sets of samples. The full estimators, PYM and NSB
(with a large alphabet size $\mathcal{A}=20K_{1}$), almost overlap with our
approximations, aPYM and aNSB. In all panels, we show results for Maximum
Likelihood (black), NSB (blue), aNSB (dashed blue), PYM (orange), aPYM (dashed
orange), and $\hat{S}_{\text{long}}$ (green) estimators. The dashed gray line
represents the true value of entropy for each of the studied distributions.

By rewriting the Maximum Likelihood estimate $S_{0}$ of Eq. (2) in terms of
coincidences (see Appendix VI.2), using the identity Eq. (19), and
approximating certain terms that are finite in the limit $d\rightarrow 1$ via
a Taylor expansion around $d\ll 1$, the mean posterior entropy, Eq. (13),
results in (see Appendix VI.3):

$\langle S|{\bf n},\alpha,d\rangle\approx\langle
S|{\mathcal{K}},\alpha,d\rangle_{a}\equiv\psi(N+\alpha+1)\\\
-\left(\frac{\alpha+K_{1}}{\alpha+N}\right)\psi(1-d)+\frac{1}{\alpha+N}\Bigg{[}.N(S_{0}-\log
N)-K_{1}\\\ +K_{2}(\log 4-1-\psi(2-d))+Q_{1}d\\\
+\mathcal{O}\left(d^{2},\sum_{m=3}\frac{K_{m}}{(m-1)^{2}}\right)\Bigg{]},$
(22)

where $\mathcal{O}(d^{2},\sum_{m=3}K_{m}/m^{2})$ means that we kept terms that
are at most linear in $d$ and at most proportional to
$\sum_{m=3}\frac{K_{m}}{(m-1)}$. Interestingly, within this approximation, the
log-likelihood and the posterior mean entropy depend on the sample size $N$,
the Maximum Likelihood entropy estimate $S_{0}$, and the three characteristics
of the coincidence vector: $K_{1},K_{2}$ and $Q_{1}$.

The final step in approximating the estimator $\hat{S}_{PYM}$, Eq. (16), is to
integrate the expected entropy for fixed hyper-parameters $\langle
S|{\mathcal{K}},\alpha,d\rangle_{a}$ over the posterior $p_{\rm
posterior}(\alpha,d|{\bf n})\propto p({\bf n}|\alpha,d)p_{\rm
prior}(\alpha,d)$ to form the Pitman-Yor mixture. Then the variance of the
resulting estimator is dominated by the contribution from the uncertainty in
the posterior distribution of the parameters $\alpha,d$, which is about $80\%$
of the total variance in our simulations.

This procedure of replacing $\langle S|\bm{n},\alpha,d\rangle$ with the
asymptotic expression $\langle S|{\mathcal{K}},\alpha,d\rangle_{a}$ in Eq.
(16) leads to a new estimator of entropy, which we call approximate PYM
estimator, or aPYM. This estimator is fully determined by just few data
statistics, $N$, $S_{0}$, $K_{1}$, $K_{2}$, and $Q_{1}$. There are also two
limiting cases of this estimator. First, by taking $d\to 0$ in Eqs. (20, 22),
we define the approximate version of the NSB limit of the PYM estimator on a
countably infinite number of possible outcomes, which we denote as aNSB. At
the other extreme, taking $\alpha\to 0$ in Eqs. (20, 22), corresponds to a
prior that favors distributions with long tails. We denote the corresponding
estimator as $\hat{S}_{\text{long}}$.

Figure 3: a: Phase diagram of the dominant tail hypothesis selected by the PYM
estimator as a function of various statistics of the data sample. The explored
statistics are the fraction of coincidences in the sample, $\Delta/N$, and
dispersion of the coincidences, $K_{2}/K_{1}$. This diagram is evaluated at
the third crucial data statistics set at
$Q_{1}=0.3\,Q_{\text{max}}=0.3(\Delta-K_{2})/2$. b: Schematic diagram that
illustrates how sample sets with different $\Delta$, $K_{1}$, and $K_{2}$ may
look like. An empty or gray circle above a state $x_{i}$ represent a single
sample for that state. Gray circles denote coincidences.

The above observation that, in the undersampled regime where
$\exp(S/2)<N<\exp(S)$, the PYM entropy estimator and its relatives are
determined approximately by just few statistics of the data,
$\\{N,S_{0},K_{1},K_{2},Q_{1}\\}$, is the main result of our paper. To
corroborate this, we explore the quality of the approximation numerically for
different distributions $\bm{q}$. Figure 2 presents results for three
distributions with different structures of tails, generated from the Pitman-
Yor Process: a distribution with an exponential tail (Fig. 2a, ${\rm
PYP}(d=0,\alpha=400)={\rm DP}(400)$), one with a mixed tail (Fig. 2b: ${\rm
PYP}(d=0.4,\alpha=100)$), and one with a long tail (Fig. 2c: ${\rm
PYP}(d=0.6,\alpha=0)$). In the lower panels we show the results of estimating
entropy for different dataset sizes using the ML estimator, the PYM estimator,
the NSB estimator with a large alphabet size $\mathcal{A}=20K_{1}$, and the
three approximations: aPYM, aNSB, and $\hat{S}_{\text{long}}$. All results are
averaged over ten sets of random samples. In all cases, the differences
between NSB and aNSB on the one hand, and PYM and aPYM on the other are
negligible, supporting the accuracy of the approximation. All four of these
estimators produce high quality estimates for all sample sizes. Further, we
also checked that the approximation of the posterior error of the estimators
is close to that of the full versions (not shown). In contrast,
$\hat{S}_{\text{long}}$ only performs well when the distribution has a long
tail, and the Maximum Likelihood never works well.

##  IV Tail-hypothesis and entropy estimation phase diagrams

The above discussion shows that the PYM estimator and its relatives work by
first estimating the most likely $\alpha$ and $d$ from the sampled data, and
then using these estimated parameters to approximate the structure of the low
probability tail (from short, to long) and hence of its contributions to the
entropy. We further showed that, in the regime of interest, the log-likelihood
of $\alpha$ and $d$ is dominated by just few statistics: $N$, $S_{0}$,
$K_{1}$, $K_{2}$, and $Q_{1}$. It is thus illustrative to understand, which
combinations of these statistics select which hypothesis on the structure of
the tail. Building the corresponding phase diagram of the selected tail
structure as a function of the data statistics is the goal of this Section.

We will consider three classes of tails: exponential ($d=0$ selected, denoted
as hypothesis $H=1$), long tail ($\alpha=0$ selected, denoted as hypothesis
$H=2$), and a mixed tails (arbitrary $\alpha$ and $d$, denoted as $H=3$). Our
goal is then to evaluate which of the three tail hypotheses has a higher
probability given the data. Long and short tail hypotheses have one parameter
each, while the mixed tail hypothesis has two parameters and contains the
other two hypotheses as special cases. Thus when evaluating the log-
likelihoods of each of the hypotheses, we must penalize them for having a
different number of parameters, which we do using Bayesian Information
Criterion [24]. To do this, we evaluate the likelihoods

$\mathcal{L}_{H}=\log p({\mathcal{K}}|\hat{\alpha},\hat{d})+\log p_{\rm
prior}(\hat{\alpha},\hat{d})-\frac{n_{H}}{2}\log N,$ (23)

where $\hat{\alpha}$ and $\hat{d}$ are the maximum likelihood values of the
parameters within each hypothesis, and $n_{H}$ is the number of parameters for
the hypothesis ($n_{H}=2$ for $H=3$, and $n_{H}=1$ otherwise). We remind the
reader that, by construction, $\hat{\alpha}=0$ for the long tail hypothesis,
$H=2$, and $\hat{d}=0$ for the short tailed hypothesis, $H=1$.

We determine the regions of the $N,S_{0},K_{1},K_{2},Q_{1}$ space, where one
of the three $\mathcal{L}_{H}$ dominates, and plot the slice of this phase
diagram in Fig. 3. Specifically, in the Figure, we vary the total number of
coincidences, $\Delta=N-K_{1}$, and the number of states with coincidences,
that is, the number of states with more than two counts, $K_{2}$. By sampling
many distributions, we empirically observe that the value $Q_{1}\sim
0.6(\Delta-K_{2})/2$ is when the rest of the $\Delta-K_{2}$ counts are
uniformly dispersed, and $Q_{1}$ tends to zero when the rest of the counts are
concentrated in a single state. Note that the maximum value $Q_{1}$ can take
is $Q_{\text{max}}=\frac{\Delta-K_{2}}{2}$. For this reason, we choose the
intermediate representative value $Q_{1}=0.3Q_{\max}=0.3\frac{\Delta-
K_{2}}{2}$.

To simplify the presentation, we plot the winning tail hypothesis as a
function of $\Delta/N$ and $K_{2}/K_{1}$. Normalized in this way, the diagram
is constrained to a square of size 1, as $0\leq\Delta/N,K_{2}/K_{1}\leq 1$. In
addition, $K_{2}\leq\Delta$, which means that the upper left corner is not
accessible. The ratio $\Delta/N$ determines how common are the coincidences,
and the ratio $K_{2}/K_{1}$ describes whether the coincidences in the data are
concentrates in a few states, or dispersed over many states (see Figure 3b).

Figure 3a show that the exponential tail hypothesis dominates when there are
many coincidences, $\Delta/N\sim 1$, or when the coincidences are dispersed,
that is $K_{2}/K_{1}\sim 1$ or $K_{2}/\Delta\sim 1$. Both cases can be
explained as corresponding to distributions that are relatively uniform on
some fixed number of states, and have zero probability elsewhere. A long tail
only dominates when the fraction of coincidences has an intermediate value,
but the coincidences are highly concentrated, $K_{2}/K_{1}\ll 1$. In other
words, in this case, there are dominant states, but a lot of samples still
fall outside of them. For other values of $\Delta/N$ and $K_{2}/K_{1}$, the
mixed tail hypothesis dominates.

Equipped with this picture of which tail hypothesis is selected by the PYM
estimator as a function of data statistics, we now can calculate how the
estimator corrects the ML entropy value $S_{0}$ for different data statistics.
Integrating the mean posterior entropy $\langle
S|\mathcal{K},\alpha,d\rangle_{a}$, Eq. (22), over our approximation of the
posterior, $p_{a}(\alpha,d|{\mathcal{K}})$, which we obtain by exponentiating
Eq. (20), we get the approximate PYM estimator $\hat{S}_{PYM,a}$. The Maximum
Likelihood estimate $S_{0}$ enters linearly in the posterior mean entropy, Eq.
(22). Thus we write

$\langle S|\mathcal{K},\alpha,d\rangle_{a}=b_{\alpha,d}\,S_{0}+\delta
S_{\alpha,d},$ (24)

where $b_{\alpha,d}$ and $\delta S_{\alpha,d}$ can be read off from Eq. (22).
Performing the integral over the approximate posterior, this becomes:

$\hat{S}=\delta S+b\,S_{0},$ (25)

where $\delta S$ and $b$ are averages of the corresponding $\alpha$\- and
$d$-dependent quantities. Thus independent of the Maximum Likelihood entropy
value, within our approximation, the PYM estimator obtains the entropy
estimate by decreasing the ML contribution from the well-sampled head of the
distribution and adding an offset that comes from the low probability tail.
This is similar to so-called partition-based entropy estimators, [17, 12, 25,
26], which divide the state space into sub-spaces, estimate entropy in each
sub-space, and then add the estimates weighted by the probability of being in
a corresponding sub-space. However, here this partitioning arises naturally
from the Bayesian framework within our approximations.

Figure 4: Corrections to entropy estimation as a function of determining data
statistics. We break down the final estimation for entropy in two parts, as
$\hat{S}=\delta
S(\nicefrac{{\Delta}}{{N}},\nicefrac{{K_{2}}}{{K_{1}}})+b(\nicefrac{{\Delta}}{{N}},\nicefrac{{K_{2}}}{{K_{1}}})\,S_{0}$,
where $\delta S$ is the additive correction and $b$ is scaling factor or
weight for the Maximum Likelihood estimate. Well-sampled distributions are
located in the upper-right corner where $\delta S=0$ and $b=1$. As in the
previous plots, we leave $Q_{1}=0.3\,(\Delta-K_{2})/2$. a: Additive correction
to entropy. b: Scaling correction to entropy.

Both the scale factor and the offset depend on the dominant $\alpha$ and $d$
contributing to the estimator, and hence on the usual statistics of the data,
$\Delta$, $K_{1}$, $K_{2}$, and $Q_{1}$. Specifically, we numerically observe
that the value of $b$ obtained from Eq. (25) satisfies

$b=\langle N/(\alpha+N)\rangle\leq 1,$ (26)

where the average is over the product of the approximate posterior obtained by
exponentiating Eq. (20) and the prior $p(\gamma)=e^{-7\gamma/100}$ with
$\gamma$ defined in Eq. 15. Note that $\alpha$ is a measure of how much
probability is concentrated in the tail. Thus the ratio $N/(\alpha+N)$
approximates the overall weight of the the well-sampled head of the
distribution, requiring to decrease the contribution to the entropy from the
head by this factor. This matches our assertion that the aPYM estimator is a
partition-based estimator, separating the head from the tail.

In Figure 4 we show results of numerical estimation of the offset $\delta S$
and the scaling factor $b$ as a function of the fraction of coincidences,
$\Delta/N$, and the dispersion of coincidences, $K_{2}/K_{1}$. As in the
previous case, we keep $Q_{1}=0.3Q_{\rm max}$. We also set $N=10^{4}$. Figure
4(a) shows that the additive term grows when the fraction of coincidences
$\Delta/N$ decreases, and when $K_{2}/K_{1}$ is small, so that coincidences
are concentrated. Both of these cases correspond to a lot of mass in the tail
(see corresponding long tail region in Figure 3(a). The largest values of
$\delta S$ occur along the boundary strip $(\Delta/N,K_{2}/K_{1}\ll 1)$ and
the boundary $K_{2}=\Delta$. Panel b shows that the scaling factor $b$ is
close to 1 in most areas, except near the boundary edge $K_{2}=\Delta$. Along
this boundary, the scaling factor becomes the largest when the number of
coincidences decreases, $\Delta/N\ll 1$. Figure 4 clearly highlights when
Bayesian corrections to the ML estimation of entropy are essential: regions of
few and concentrated coincidences.

##  V Discussion

The major finding of this work is an excellent approximation for the PYM
estimator, one of the best Bayesian entropy estimators, and its various
relatives (such as NSB). The approximation simplifies the numerics
considerably. Crucially, the approximation also shows that the output of the
PYM entropy estimator depends on just a few statistics of the data, namely the
maximum likelihood (ML) entropy estimate, the fraction of coincidences
$\Delta/N$, and the dispersion of coincidences $K_{1}/K_{2}$, and $Q_{1}$. We
showed that that workflow of the estimator can be interpreted as first
estimating the parameters $d$ and $\alpha$ based on the aforementioned
statistics, and with them the tail structure and the total weight of the tail.
Then the estimator rescales the ML entropy estimate by the weight of the well-
sampled head of the distribution, and adds to it the estimated entropy of the
tail. The phase diagrams of which tail structure the estimator selects, Fig.
3, and how it corrects the ML estimate, Fig. 4, illustrate these points.

Early work of Ma [15] showed that when states are equiprobable, in the under-
sampled regime, the coincidences in counts can help with the inference of the
entropy of a system. Later Nemenman [16] showed that in the severely under-
sampled regime ($K_{1}$ close to $N$), entropy estimation depends on the
number of coincidences $K_{1}$. Further, he pointed out how reliable entropy
estimates may be obtained by partitioning the overall state space of the
variable into sub-spaces with similar sampling properties [26]. Here we extend
these results to the whole regime where entropy estimation is challenging for
multinomial observations, $\exp(S/2)<N<\exp(S)$, by approximating the more
general PYM estimator. Our identification of the small set of statistics,
which define the output of the estimator, lifts the veil from its inner
workings, allowing for a simple, semi-analytical estimation procedure. In
particular, this allows us to predict if a particular estimator will be biased
simply by looking at the values of the select statistics of the data.

How to match a priori assumptions about the underlying distributions to the
data to allow for an unbiased estimation of quantities of interest—such as
entropy [17, 10] or the mutual information [27]— is an open problem [28]. It
requires understanding the relation between the a priori assumptions and the
data features that influence the inference. In this work, we build such a link
for entropy estimation, and we hope that similar links might exist for other
difficult estimation problems.

###### Acknowledgements.

This work was funded, in part, by the Simons Investigator award (IN), the
Simons-Emory International Consortium (AR and IN), NSF Grant 1822677 (DGH and
IN), NIH Grant 2R01NS084844 (AR, IN), the International Physics of Living
Systems Network (NSF Grant 1806833, AR). IN would like to acknowledge
hospitality of the Aspen Center for Physics, funded in part by NSF Grant
1607611.

## References

  * Shannon [1948] C. E. Shannon, A mathematical theory of communication, Bell System Technical Journal 27, 379 (1948).
  * Cover and Thomas [2012] T. M. Cover and J. A. Thomas,  _Elements of information theory_ (John Wiley & Sons, 2012).
  * Strong  _et al._ [1998] S. P. Strong, R. Koberle, R. R. d. R. van Steveninck, and W. Bialek, Entropy and information in neural spike trains, Physical Review Letters 80, 197 (1998).
  * Paninski [2003] L. Paninski, Estimation of entropy and mutual information, Neural Computation 15, 1191 (2003).
  * Miller [1955] G. Miller, Note on the bias of information estimates, Information theory in psychology: Problems and methods (1955).
  * Grassberger [2003] P. Grassberger, Entropy estimates from insufficient samplings, arXiv preprint physics/0307138 (2003).
  * Berry II  _et al._ [2013] M. J. Berry II, G. Tkačik, J. Dubuis, O. Marre, and R. A. da Silveira, A simple method for estimating the entropy of neural activity, Journal of Statistical Mechanics: Theory and Experiment 2013, P03015 (2013).
  * Wolpert and Wolf [1995] D. H. Wolpert and D. R. Wolf, Estimating functions of probability distributions from a finite set of samples, Physical Review E 52, 6841 (1995).
  * Nemenman  _et al._ [2002] I. Nemenman, F. Shafee, and W. Bialek, Entropy and inference, revisited, in  _Advances in neural information processing systems_ (2002) pp. 471–478.
  * Archer  _et al._ [2014] E. Archer, I. M. Park, and J. W. Pillow, Bayesian entropy estimation for countable discrete distributions, The Journal of Machine Learning Research 15, 2833 (2014).
  * Chao and Shen [2003] A. Chao and T.-J. Shen, Nonparametric estimation of shannon’s index of diversity when there are unseen species in sample, Environmental and ecological statistics 10, 429 (2003).
  * Chao  _et al._ [2013] A. Chao, Y. Wang, and L. Jost, Entropy and the species accumulation curve: a novel entropy estimator via discovery rates of new species, Methods in Ecology and Evolution 4, 1091 (2013).
  * Cerquetti [2019] A. Cerquetti, Exact good-turing characterization of the two-parameter poisson-dirichlet superpopulation model, arXiv preprint arXiv:1901.09665 (2019).
  * [14] A. Antos and I. Kontoyiannis, Estimating the entropy of discrete distributions, in  _Proceedings. 2001 IEEE International Symposium on Information Theory (IEEE Cat. No.01CH37252)_, p. 45.
  * Ma [1981] S.-k. Ma, Calculation of entropy from data of motion, Journal of Statistical Physics 26, 221 (1981).
  * Nemenman [2011] I. Nemenman, Coincidences and estimation of entropies of random variables with large cardinalities, Entropy 13, 2013 (2011).
  * Nemenman  _et al._ [2004] I. Nemenman, W. Bialek, and R. d. R. van Steveninck, Entropy and information in neural spike trains: Progress on the sampling problem, Physical Review E 69, 056111 (2004).
  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , ninth dover printing, tenth gpo printing ed. (Dover, New York City, 1964).
  * Minka [2000] T. Minka, Estimating a dirichlet distribution (2000).
  * Pitman  _et al._ [1997] J. Pitman, M. Yor, _et al._ , The two-parameter poisson-dirichlet distribution derived from a stable subordinator, The Annals of Probability 25, 855 (1997).
  * Ferguson [1973] T. S. Ferguson, A bayesian analysis of some nonparametric problems, The Annals of Statistics 1, 209 (1973).
  * Ishwaran and James [2001] H. Ishwaran and L. F. James, Gibbs sampling methods for stick-breaking priors, Journal of the American Statistical Association 96, 161 (2001).
  * Ishwaran and James [2003] H. Ishwaran and L. F. James, Generalized weighted chinese restaurant processes for species sampling mixture models, Statistica Sinica 13, 1211 (2003).
  * Schwarz [1978] G. Schwarz, Estimating the dimension of a model, The Annals of Statistics 6, 461 (1978).
  * Srivastava  _et al._ [2017] K. H. Srivastava, C. M. Holmes, M. Vellema, A. R. Pack, C. P. H. Elemans, I. Nemenman, and S. J. Sober, Motor control by precisely timed spike patterns, Proceedings of the National Academy of Sciences 114, 1171 (2017).
  * Nemenman  _et al._ [2015] I. Nemenman, M. E. Wall, and C. E. Strauss, Of fishes and birthdays: Efficient estimation of polymer configurational entropies, arXiv preprint arXiv:1502.02364 (2015).
  * Hernández and Samengo [2019] D. G. Hernández and I. Samengo, Estimating the mutual information between two discrete, asymmetric variables with limited samples, Entropy 21, 623 (2019).
  * Hernández and Samengo [2022] D. G. Hernández and I. Samengo, Inferring a property of a large system from a small number of samples, Entropy 24, 125 (2022).

##  VI Appendix

###  VI.1 Marginal likelihood approximation for a Pitman-Yor process

In this Appendix we show how to approximate the marginal posterior of a
Pitman-Yor process in the regime $K_{1}\lesssim N\leq\exp(S)$. We start by
manipulating each term in the logarithm of the evidence $\mathcal{L}=\log
p(\mathbf{n}|\alpha,d)$ from Eq. (17),

$\mathcal{L}(\mathbf{n}|\alpha,d)=\sum_{l=1}^{K_{1}-1}\log(\alpha+ld)+\\\
\sum_{i=1}^{K_{1}}\log\Gamma(n_{i}-d)-K_{1}\log\Gamma(1-d)\\\
+\log\Gamma(1+\alpha)-\log\Gamma(N+\alpha).$ (27)

To simplify the first term in Eq. (27), we rewrite it in terms of coincidences
$K_{1}$ as follows:

$I_{1}=\displaystyle\sum_{l=1}^{K_{1}-1}\log(\alpha+ld)\\\
=\sum_{l=1}^{K_{1}-1}\left[\log
d+\log\left(\frac{\alpha}{d}+l\right)\right]\displaystyle=(K_{1}-1)\log d\\\
+\sum_{l=1}^{K_{1}-1}\left[\log\Gamma\left(\frac{\alpha}{d}+l+1\right)-\log\Gamma\left(\frac{\alpha}{d}+l\right)\right]\\\
\displaystyle=(K_{1}-1)\log
d+\log\Gamma\left(\frac{\alpha}{d}+K_{1}\right)-\log\Gamma\left(\frac{\alpha}{d}+1\right).$
(28)

In order to rewrite the rest of the terms of Eq. (27) in terms of various
coincidence statistics, we use the identity Eq. (19). Joining the second and
third terms in Eq. (27) and rewriting them in terms of count multiplicities
yields

$\displaystyle\sum_{i=1}^{K_{1}}\log\Gamma(n_{i}-d)-K_{1}\log\Gamma(1-d)$
$\displaystyle=-K_{2}\log\Gamma(1-d)+\sum_{m=2}(K_{m}-K_{m+1})\log\Gamma(m-d)$
$\displaystyle=\sum_{m=2}K_{m}\left[\log\Gamma(m-d)-\log\Gamma(m-1-d)\right]$
$\displaystyle=\sum_{m=2}K_{m}\log(m-1-d)$
$\displaystyle=K_{2}\log(1-d)+Q(d),$ (29)

where

$Q(d)=\sum_{m=3}^{m_{f}}K_{m}\log(m-1-d).$ (30)

where $m_{f}$ denotes the largest occupancy of any state in the sample. Since
the domain of $0\leq d<1$ is small, $Q(d)$ is approximately linearly varying
with $d$, so that we can expand it around $d=0$:

$Q(d)=Q(0)-\sum_{j=1}\left[\sum_{m=3}\frac{K_{m}}{(m-1)^{j}}\right]\frac{d^{j}}{j}\\\
\approx Q(0)-\left[\sum_{m=3}\frac{K_{m}}{m-1}\right]d\\\
-\frac{1}{2}\left[\sum_{m=3}\frac{K_{m}}{(m-1)^{2}}\right]d^{2}+\mathcal{O}(Q_{3}),\\\
=Q_{0}-Q_{1}d-\frac{1}{2}Q_{2}d^{2}+\mathcal{O}(Q_{3}).$ (31)

where

$Q_{j}=\sum_{m=3}\frac{K_{m}}{(m-1)^{j}}$ (32)

for $j\geq 1$. As $d$ approaches $1$, the term $K_{2}\log(1-d)$ goes to
infinity, which renders any error in the Taylor expansion of $Q(d)$
irrelevant. This makes the approximations above useable even if we ignore
$\mathcal{O}(d^{2})$ terms.

Putting all of the approximations above together, the ensuing approximate
logarithm of the evidence $\mathcal{L}(\mathbf{n}|\alpha,d)$ is

$\mathcal{L}(\mathbf{n}|\alpha,d)\approx(K_{1}-1)\log
d+\log\Gamma\left(\frac{\alpha}{d}+K_{1}\right)\\\
-\log\Gamma\left(\frac{\alpha}{d}+1\right)+\log\Gamma(1+\alpha)-\log\Gamma(N+\alpha)\\\
+K_{2}\log(1-d)-Q_{1}d+\mathcal{O}\left(d^{2}\sum_{m=3}\frac{K_{m}}{(m-1)^{2}}\right),$
(33)

up to an additive constant. This is Eq. (20) in the main text.

###  VI.2 Maximum likelihood Entropy in terms of coincidences

To relate the conditional entropy, Eq. (13), to the Maximum Likelihood entropy
estimator $S_{0}$, we need to rewrite the latter in terms coincidences.
Utilizing the identity Eq. (19), we write

$N\left[S_{0}-\log N\right]=-\sum_{i}n_{i}\log n_{i}=\\\
-\sum_{m=2}(K_{m}-K_{m+1})m\log m\\\ =-K_{2}(2\log 2)\\\
-\sum_{m=3}K_{m}\big{[}m\log m-(m-1)\log(m-1)\big{]}.$ (34)

Rewriting the expression in brackets as

$m\log m-(m-1)\log(m-1)=1+\psi(m)+\mathcal{O}(m^{-2}).$ (35)

and plugging this into Eq. (34), we finally obtain,

$N\left[S_{0}-\log N\right]=-K_{2}\log 4-(N-K_{1}-K_{2})-\\\
\sum_{m=3}K_{m}\psi(m)+\mathcal{O}\left(\sum_{m}K_{m}/m^{2}\right).$ (36)

###  VI.3 Mean posterior entropy approximation for the Pitman-Yor Process

Similar to Appendix VI.1, here we approximate the posterior entropy, Eq. 13,
in the limit of small $d$. To simplify the notation, we use the shorthand
$S=\langle S|\bm{n},\alpha,d\rangle$ in this Appendix. Rearranging Eq. (13),
we obtain

$\displaystyle(\alpha+N)\left[S-\psi(N+\alpha+1)\right]=\\\
-\alpha\,\psi(1-d)-K_{1}\,d\,\psi(1-d)-\sum_{i}(n_{i}-d)\psi(n_{i}+1-d).$ (37)

We now again use Eq. (19) and a Taylor expansion in small $d$ to rewrite the
last term on the right hand side of Eq. (37):

$\displaystyle K_{1}\,d\,\psi(1-d)-\sum_{i}(n_{i}-d)\psi(n_{i}+1-d)$
$\displaystyle=K_{1}\,d\,\psi(1-d)-\sum_{m=1}(K_{m}-K_{m+1})(m-d)\psi(m+1-d)$
$\displaystyle=-\sum_{m=1}K_{m}\left[(m-d)\psi(m+1-d)-(m-1-d)\psi(m-d)\right]$
$\displaystyle=-\sum_{m=1}K_{m}\left[1+\psi(m-d)\right]$
$\displaystyle=-\sum_{m=1}K_{m}-\sum_{m=1}K_{m}\psi(m-d)$
$\displaystyle=-N-K_{1}\psi(1-d)-K_{2}\psi(2-d)-\sum_{m=3}K_{m}\psi(m-d).$
(38)

where we used $\psi(m+1-d)=\left(\psi(m-d)+\frac{1}{m-d}\right).$

Since $m\geq 3$, we can Taylor expand the sum in this last term around $d=0$
to obtain

$\sum_{m=3}K_{m}\psi(m-d)\approx\sum_{m=3}K_{m}\psi(m)+d\,\sum_{m=3}K_{m}\psi^{\prime}(m)\\\
+\mathcal{O}(d^{2}\sum_{m}K_{m}\psi^{\prime\prime}(m)).$ (39)

Now using the relations $\psi^{\prime}(m)=\frac{1}{m-1}+\mathcal{O}(m^{-2})$
and the expression for $\sum_{m=3}K_{m}\psi(m)$ in Eq. (36), we rewrite Eq.
(39) as

$\sum_{m=3}K_{m}\psi(m-d)\\\ \approx K_{2}\log
4+(N-K_{1}-K_{2})-N\left[S_{0}-\log N\right]\\\
+d\,\sum_{m=3}\frac{K_{m}}{m-1}+\mathcal{O}(d^{2},\sum_{m=3}K_{m}/m^{2}),$
(40)

where $\mathcal{O}(d^{2},\sum_{m=3}K_{m}/m^{2})$ means that we kept terms that
are at most linear in $d$ and whose summands are at most proportional to
$\sum_{m=3}K_{m}/m$. Plugging these approximation in Eq. (38) and noticing
that $Q_{1}=\sum_{m=3}\frac{K_{m}}{m-1}$, we obtain

$(\alpha+N)\left[S-\psi(N+\alpha+1)\right]=\\\ N(S_{0}-\log
N)-\alpha\,\psi(1-d)+K_{1}\left[-1-\,\psi(1-d)\right]\\\
+K_{2}\left[-1-\psi(2-d)+\log
4\right]-Q_{1}\,d+\mathcal{O}(d^{2},\sum_{m=3}K_{m}/m^{2}),$ (41)

which after isolating $S$ becomes Eq. (22) of the main text.


# Modifying quantum optical states by zero-photon subtraction

C.M. Nunn  J.D. Franson  T.B. Pittman  Department of Physics, University of
Maryland Baltimore County, Baltimore, MD 21250, USA

###### Abstract

The process of single-photon subtraction (SPS) is known to dramatically alter
the properties of certain quantum optical states. Somewhat surprisingly,
subtracting zero photons can also modify quantum states and has practical
applications in quantum communication. Here we experimentally investigate
zero-photon subtraction (ZPS) using a wide variety of input states and
conditional measurements based on actively detecting zero photons in one
output port of a variable beamsplitter. We find that SPS and ZPS can exhibit
complementary behavior depending on the photon statistics of the input states,
and highlight deeper connections with Mandel’s $Q$-parameter for classifying
quantum states.

##  I Introduction

In experimental quantum optics, the bosonic annihilation operator $\hat{a}$
can be realized through the process of “single-photon subtraction” (SPS) [1,
2]. When the input state contains a definite number of photons $\ket{n}$, this
operation transforms the state as $\hat{a}\ket{n}\rightarrow\sqrt{n}\ket{n-1}$
in the usual way, corresponding to the simple removal of one photon from the
state. However, when the input is a superposition of different number states,
the SPS process can lead to counterintuitive results [3, 4]. For example,
consider the input state
$\ket{\psi}_{in}=\frac{1}{\sqrt{2}}\left(\ket{1}+\ket{5}\right)$, which has a
mean number of photons $\langle\hat{n}\rangle=3$ (where the number operator
$\hat{n}\equiv\hat{a}^{\dagger}\hat{a}$). Applying $\hat{a}$ to this state
leads to
$\ket{\psi}_{out}=\frac{1}{\sqrt{6}}\left(\ket{0}+\sqrt{5}\ket{4}\right)$,
which has $\langle\hat{n}\rangle=3.\bar{3}$. In this sense, subtracting a
single photon from the state has actually increased the mean number of photons
[5, 6].

Figure 1: An implementation of “zero-photon subtraction” (ZPS) via conditional
measurements on a beamsplitter. A superposition state $\ket{\psi}_{in}$ with
expected photon number $\langle\hat{n}\rangle_{in}$ is prepared in the input
mode of a beamsplitter with reflectance $R$. In contrast to single-photon
subtraction (SPS), ZPS requires heralding on the detection of zero photons in
the reflected mode [7]. Heralding on zero photons yields the attenuated state
$\ket{\psi}_{out}$ with reduced mean photon number
$\langle\hat{n}\rangle_{out}<\langle\hat{n}\rangle_{in}$. The degree of
attenuation depends on both $R$ and the photon number statistics of the input
state.

In a similarly counterintuitive way, subtracting zero photons from a state can
actually decrease the mean number of photons. Figure 1 shows an implementation
of this “zero-photon subtraction” (ZPS) process using a beamsplitter with
reflectance $R$, and conditional measurements. A pulsed input state
$\ket{\psi}_{in}$ passes through the beamsplitter, and the transmitted output
$\ket{\psi}_{out}$ is heralded by the successful detection of zero photons in
the reflected mode. Despite no photons being physically removed from the
system, ZPS results in
$\langle\hat{n}\rangle_{out}<\langle\hat{n}\rangle_{in}$ for all but pure Fock
states [8]. Importantly, the ZPS process in Fig. 1 can be used to implement a
probabilistic noiseless attenuation protocol that is useful for quantum
communications [9, 10, 11].

As highlighted by the structure of Fig. 1, the key difference between SPS and
ZPS is a heralding signal based on the detection of one vs. zero photons,
respectively. While SPS has been experimentally studied extensively [12], ZPS
has only been briefly observed for super-Poissonian (i.e., thermal) states,
and with fixed values of beamsplitter reflectance $R$ [13, 14, 15, 16, 17, 18,
19]. In this paper, we systematically study ZPS for examples of super-
Poissonian, sub-Poissonian, and coherent state inputs, all as a function of
beamsplitter reflectance ranging from $R=0\rightarrow 1$. The observed trends
in attenuation demonstrate some complementary aspects of ZPS and SPS that
depend on the photon number distributions, and highlight the role of losses
and detector efficiency when heralding on zero photons in ZPS.

The remainder of the paper is structured as follows: in Section II we provide
a detailed theoretical background for ZPS, and introduce an experimentally
accessible parameter $K$ that can be used to quantify the degree of
attenuation. In Section III we describe our experimental system, which uses
(1) a conventional pulsed parametric down-conversion (PDC) source to produce
the desired input states [20], (2) a variable evanescent-mode fiber coupler to
continuously vary $R$ [21], and (3) the ability to actively herald on zero
photons using commercial single-photon detectors [7]. In Section IV we analyze
and discuss the experimental results, and briefly describe classical analogues
that provide some additional insight into the observed attenuation effects.
Finally, we summarize our study and conclude in Section V.

##  II Zero-Photon Subtraction

The process of zero-photon subtraction (ZPS) illustrated in Figure 1 was first
proposed as a method of noiseless attenuation by Mičuda et al. [9]. This
transformation can be defined by its action on Fock states $\ket{n}\rightarrow
t^{n}\ket{n}$ with beamsplitter transmittance $T=|t|^{2}$. For an arbitrary
input state $\hat{\rho}=\sum_{m,n=0}^{\infty}\rho_{mn}\ket{m}\bra{n}$,
noiseless attenuation yields the following expected photon number in the
output [8]:

$\langle\hat{n}\rangle_{out}=\frac{\sum_{n}n\rho_{nn}T^{n}}{\sum_{n}\rho_{nn}T^{n}}$
(1)

When $T=1$, we regain the expected photon number of the original state with no
attenuation, $\langle\hat{n}\rangle_{in}=\sum_{n}n\rho_{nn}$. Remarkably, when
$T<1$ it can be seen that
$\langle\hat{n}\rangle_{out}<\langle\hat{n}\rangle_{in}$ for all but pure Fock
states by differentiating Eq. 1 with respect to $T$:

$\frac{d\langle\hat{n}\rangle_{out}}{dT}=\frac{1}{T}\left[\frac{\sum_{n}(n-\langle\hat{n}\rangle_{out})^{2}\rho_{nn}T^{n}}{\sum_{n}\rho_{nn}T^{n}}\right]\equiv\frac{\langle(\Delta
n)^{2}\rangle_{out}}{T}\geq 0,$ (2)

and seeing that $\langle\hat{n}\rangle_{out}$ increases monotonically on the
interval $T\in(0,1]$. Here $\langle(\Delta n)^{2}\rangle_{out}$ is the photon
number variance of the transformed state.

Equation 2 is analogous to the result derived by Ueda et al. for stationary
fields [3], and it suggests that the degree of attenuation is closely related
to the photon number statistics of the input state. Experimentally, it is
convenient to quantify the degree of attenuation as a function of reflectance
$R$ with the following ratio:

$K(R)\equiv\frac{\langle\hat{n}\rangle_{out}}{(1-R)\langle\hat{n}\rangle_{in}}$
(3)

The denominator $(1-R)\langle\hat{n}\rangle_{in}$ simply corresponds to
ordinary attenuation by a beamsplitter, in which a fraction $T=1-R$ of the
photons are transmitted on average. Thus, $K(R)$ compares the mean photon
number of the heralded ZPS state $\langle\hat{n}\rangle_{out}$ to that of the
“ordinary” output state with no conditional measurements.

The relative attenuation function $K(R)$ contains information about the photon
number distribution and higher-order correlations. Most importantly, one can
derive from Eqns. 1 to 3 that:

$\frac{dK}{dR}\Bigr{|}_{R=0}=1-\frac{\langle(\Delta
n)^{2}\rangle_{in}}{\langle\hat{n}\rangle_{in}}\equiv-Q_{in}$ (4)

where $Q_{in}$ is Mandel’s $Q$-parameter for the input state [22].

The above result highlights an important connection between ZPS and typical
SPS. In the limit of low beamsplitter reflectance $R$, SPS is equivalent to
the annihilation operator $\hat{a}$ 111Historically, “photon subtraction”
refers to application of the annihilation operator, predating its realization
by conditional measurements in the output of a beamsplitter (BS) [12].
However, the term often encompasses the more general BS transformations which
relax the weak reflectance requirement (e.g., [33, 34]), which we follow here
in our discussion of ZPS., and increases the mean photon number of some states
as demonstrated in Sec. I. More precisely, this so-called “photon excess” is
given exactly by the $Q$-parameter, such that
$Q_{in}=\langle\hat{n}\rangle_{out}-\langle\hat{n}\rangle_{in}$ [4]. Thus, the
mean photon number of super-Poissonian states $(Q>0)$ counterintuitively
increases after performing SPS with a weakly reflecting beamsplitter. Equation
4 links this property of SPS to the behavior of ZPS in the same regime of
$R\ll 1$. For ZPS, the $Q$-parameter determines the initial slope $dK/dR$, and
thus deviations from $K=1$ as $R$ increases from zero. We can therefore say
super-Poissonian states exhibit a complementary “photon deficit” $(K<1)$ after
ZPS in this regime, such that the mean photon number is reduced below that of
ordinary attenuation.

In the same way that SPS has unique consequences for sub-, super- and
Poissonian states [6], it is also natural to investigate these three classes
of states for ZPS. Our experiment will examine the following cases: (1)
coherent states $\ket{\alpha}$, which possess Poissonian statistics; (2) the
single-mode squeezed vacuum state (SMSV) $\ket{\xi}$, which is super-
Poissonian; and (3) a single-photon Fock state $\ket{1}$, which is sub-
Poissonian. As detailed in Section III, the experimentally prepared single-
photon state is actually a mixture that includes the vacuum term,
$\hat{\rho}_{1}=(1-\beta)\ket{0}\bra{0}+\beta\ket{1}\bra{1}$. We can calculate
the expected relative attenuation for each of the three input states:

$\displaystyle K^{(\alpha)}(R)$ $\displaystyle=1$ (5) $\displaystyle
K^{(\xi)}(R)$ $\displaystyle\approx 1-R$ (6) $\displaystyle
K^{(\hat{\rho}_{1})}(R)$ $\displaystyle=\frac{1}{1-\beta R}$ (7)

where the approximation for the SMSV $\ket{\xi}$ in Eq. 6 holds for weak
squeezing.

##  III Experiment

Figure 2: Zero-photon subtraction (ZPS) experiment, in four stages. (1) Input
State Preparation– One of three input states is generated, $\ket{\alpha}$,
$\ket{\xi}$ or $\hat{\rho}_{1}$. Coherent states $\ket{\alpha}$ are produced
directly by an ultrafast pulsed laser (100 MHz, 780 nm). These pulses also
undergo second harmonic generation (SHG) to serve as 390 nm pump pulses for
Type-I parametric down-conversion (PDC) using a $\beta$-barium borate (BBO)
crystal. The resulting photon pairs are coupled into a Hong-Ou-Mandel (HOM)
interferometer to produce either $\ket{\xi}$ or $\hat{\rho}_{1}$ as described
in the text. (2) Attenuation– The input state enters a fiber-based variable
beamsplitter (VBS) with reflectance $R$. (3) Measurement– Each VBS output is
measured by single-photon detectors $D_{1}$ and $D_{2}$, with overall channel
efficiencies $\eta_{1}$ and $\eta_{2}$. The heralding detector $D_{1}$ can be
translated across the mode of interest by displacement $\Delta x$. Detection
events and the $D_{ref}$ reference signal are recorded by time-to-digital
converters (TDCs). (4) Post-selection– Time tags from all detections are used
to measure $D_{2}$ counting rates with and without post-selecting on “no-
click” events at $D_{1}$. Abbreviations: DM– dichroic mirror used to isolate
UV pump pulses; L– various lenses; $\Delta t$– glass wedge time delay; IF–
narrowband interference filters centered near 780 nm.

The full ZPS experiment is shown in Figure 2. As summarized in the first
panel, one of each type of input state (sub-, super- and Poissonian) is
prepared with a combination of standard techniques in quantum optics.

In the case of Poissonian statistics, coherent states $\ket{\alpha}$ are
prepared with a mode-locked fiber laser (Menlo Systems C-Fiber 780), which
generates a train of ultrashort pulses with a repetition rate of 100 MHz and a
center wavelength of 780 nm. These pulses are coupled into a single-mode fiber
and attenuated for use as ZPS input states.

Super-Poissonian SMSV states $\ket{\xi}$ are prepared with parametric down-
conversion (PDC) and a Hong-Ou-Mandel (HOM) interferometer [24]. The 780 nm
pulse train is first frequency doubled and used as a pump for Type-I PDC using
a $\beta$-barium borate (BBO) crystal. The resulting photon pairs are coupled
into single-mode fibers, and then combined in a 50-50 fiber coupler serving as
the HOM interferometer. The relative time delay $\Delta t$ between photons is
controlled with a pair of translating glass wedges before one of the fibers.
When $\Delta t=0$, interference ideally produces two disentangled SMSV states
in the HOM outputs [25]. With our low pump power, the PDC photon pair
production rate of $\sim 10^{-4}$ per pulse ensures we are in the weak
squeezing limit where Eq. 6 holds.

Using the same setup with noncollinear PDC and the HOM interferometer, we can
also generate heralded single-photon states with sub-Poissonian statistics.
First, a large time delay is introduced in one input of the interferometer,
eliminating HOM interference. Next, a single-photon detector $D_{0}$ without
photon-number resolution (non-PNR) is coupled to one output. When $D_{0}$
detects exactly one photon with a “click,” the twin photon is heralded in the
other mode (offset by the delay $\pm\Delta t$). Alternatively, a “click” could
result from two photons hitting $D_{0}$, heralding zero photons in the output.
The ideal result is a mixture
$\hat{\rho}_{1}=(1-\beta)\ket{0}\bra{0}+\beta\ket{1}\bra{1}$, with
$\beta=2/3$.

ZPS is performed at a variable beamsplitter (VBS), implemented with a tunable
fiber coupler [21]. The input state $\ket{\psi}_{in}$ (i.e., $\ket{\alpha}$,
$\ket{\xi}$ or $\hat{\rho}_{1}$) enters one input of the VBS, and the two
outputs are routed to detection channels. Each channel includes a free-space
U-bench with 25-nm-bandwidth rectangular bandpass filters centered near 780
nm, then coupled into multimode fibers and directed to single-photon counting
modules (SPCMs) $D_{1}$ and $D_{2}$ (silicon avalanche photodiodes, Excelitas
SPCM-AQ4C). The auxiliary heralding detector $D_{0}$ has a similar channel not
shown in Figure 2, with a more narrow 10-nm-bandwidth filter to increase
heralding efficiency [26, 27]. To serve as a universal clock for all “click”
and “no-click” events, all detection signals are recorded alongside a 100 MHz
mode-locking reference signal from an additional detector $D_{ref}$, using
time-to-digital converters (TDCs) with 81 ps timebin resolution (IDQuantique,
model ID801). All detection events are stored as time tags and processed using
the techniques described in Ref. [7].

The counting statistics of ZPS states are observed by post-selecting on “no-
click” events, in which $D_{ref}$ registers a pulse but the heralding detector
$D_{1}$ measures zero photons. After a 20-second exposure, the mean counting
rate at $D_{2}$ is calculated with and without this post-selection. Then the
$D_{2}$ dark count rate ($\sim$80 Hz, after filtering [28]) is subtracted from
each of these values, and their ratio is taken to obtain $K$. This is repeated
for multiple values of VBS reflectance $R$, revealing the behavior of $K(R)$
for each state.

Each stage of the experiment introduces losses which must be taken into
account for our analysis. Returning to Figure 1, we can group all losses into
three distinct channels: the input mode of the main beamsplitter (VBS); the
reflected auxiliary mode, containing heralding detector $D_{1}$; and the
transmitted output mode, containing the photon-counting detector $D_{2}$.
Input losses are primarily due to coupling free-space photon pairs from the
PDC source into single-mode fibers, as well as fiber connector losses at the
VBS. The fiber-coupling efficiency and connector transmission are denoted
$\kappa_{\text{PDC}}$ and $\kappa_{f}$, respectively. Additional losses after
the VBS are contained in the effective detector efficiencies $\eta_{1}$ and
$\eta_{2}$, illustrated in the third panel of Fig. 2.

As defined in Eq. 3, $K(R)$ is unaffected by losses in the output mode with
detector $D_{2}$. This can differ for non-PNR detectors as shown in the
Appendix, but these effects are negligible in our experiment. However, losses
in the heralding mode introduce unwanted noise that alters our counting
statistics [7]. Similarly, input losses introduce noise that alter the photon
statistics of the initial states. Even so, the resulting mixed states can be
analyzed with the same measurement of $K$, which only depends on diagonal
terms $\rho_{nn}$ in a full description of the state. Consequently, we can
modify our equations of $K(R)$ to include all losses (see Appendix for
details):

$\displaystyle K_{exp}^{(\alpha)}(R)$ $\displaystyle=1$ (8) $\displaystyle
K_{exp}^{(\xi)}(R)$ $\displaystyle\approx
1-\kappa_{\text{PDC}}\kappa_{f}R\eta_{1}$ (9) $\displaystyle
K_{exp}^{(\hat{\rho}_{1})}(R)$
$\displaystyle=\frac{1}{1-\beta\kappa_{f}R\eta_{1}}$ (10)

Results for the coherent state $\ket{\alpha}$ remain unchanged from Eq. 5 to
Eq. 8. In the case of the SMSV (Eq. 9), the probability of multipair emission
is negligible, and so the HOM interferometer output is dominated by zero- and
two-photon terms. Loss before ($\kappa_{\text{PDC}}$) and after the
interferometer ($\kappa_{f}$) introduce a significant single-photon component,
but the altered statistics remain super-Poissonian. For the state
$\hat{\rho}_{1}$, the existing single-photon term is similarly reduced by
$\kappa_{f}$ but remains sub-Poissonian. The initial single-photon probability
$\beta$ is also degraded by dark counts at $D_{0}$ and interferometer losses,
lowering it from the ideal value of $2/3$. In all cases, finite heralding
efficiency $\eta_{1}$ has the same effect on the measured value of $K$ as the
overall input losses.

To experimentally determine coupling values and detector efficiencies in our
system, and to align the apparatus for input state preparation, we first
perform a series of standard HOM tests [24] and channel loss measurements. We
bypass $D_{0}$ and the VBS in Figure 2 and perform a coincidence measurement
with the $D_{1}$ and $D_{2}$ channels connected directly to the HOM
interferometer outputs. We observe a HOM dip with 98% visibility with this
arrangement. Approximate Klyshko efficiencies [29], apart from the
interferometer coupling efficiency of $\kappa_{\text{PDC}}\approx 0.50$, are
found to be $\eta_{1}\approx 0.32$ and $\eta_{2}\approx 0.28$. This is
consistent with nominal SPCM detector efficiencies of $\sim 50\%$ at 780 nm
and U-bench transmission of $\sim 65\%$ and $\sim 60\%$. The values of
$\kappa_{f}\approx 0.86$ and $\beta\approx 0.38$ are determined in the
analysis of the main experiment.

##  IV Results and Discussion

Our main results are shown in Figure 3. The relative attenuation $K$ induced
by ZPS is shown as a function of reflectance $R$ for all three input states.
Each state exhibits very distinct behavior in agreement with equations 8-10.
For the coherent state shown in Fig. 3(a), $K=1$ for all $R$, indicating that
ZPS counting rates are identical to those of ordinary attenuation. This case
provides a benchmark for our experiment, and can be understood by the well-
known fact that a coherent state entering a beamsplitter produces two
uncorrelated coherent states in the outputs [6, 20, 25]. Consequently,
conditional measurements like those in ZPS have no effect on the output state
in this case.

Figure 3: Experimental measurements of relative attenuation $K(R)$ for the
three cases of (a) a coherent state $\ket{\alpha}$, (b) a SMSV state
$\ket{\xi}$, and (c) a heralded single photon state $\hat{\rho}_{1}$. The
measured data points (red circles) in plots (a)-(c) show distinct trends that
depend on the photon number statistics of the given state: the benchmark case
$\ket{\alpha}$ shows $K=1$ for all reflectance $R$, while $\ket{\xi}$ and
$\hat{\rho}_{1}$ have a negative and positive slope, respectively. In panel
(b), a fit to the data using Eq. 9 (black dashed curve) gives the value
$(\kappa_{\text{PDC}})(\kappa_{f})(\eta_{1})=0.14$, corresponding to an
overall loss of 86%. In panel (c), a similar fit to the data using Eq. 10
corresponds to an overall loss of 73% with single-photon probability
$\beta\approx 0.38$. In both panels (b) and (c), the blue (dot-dashed) and
green (dotted) theoretical curves show more pronounced effects of ZPS that
would be observed for 50% and 0% overall loss, respectively. For comparison,
panel (d) shows theoretical attenuation of an ideal Fock state $\ket{n}$
calculated for the same overall loss values as in panel (c).

In Fig. 3(b), measurements of $K$ for the SMSV state $\ket{\xi}$ are shown.
The data trend exhibits a negative initial slope in accordance with Eq. 4
$(Q>0)$. This “photon deficit” $K<1$ increases linearly with reflectance $R$.
This attenuation, however, is limited by heralding efficiency and losses. A
fit to the data using Eq. 9 shows that as $R\rightarrow 1$,
$K^{(\xi)}_{min}=0.861\pm 0.003$. This is consistent with the product of
efficiencies $(\kappa_{\text{PDC}})(\kappa_{f})(\eta_{1})\approx 0.14$ (i.e.,
overall loss of 86%). For comparison, the two theoretical curves in Fig. 3(b)
show the stronger attenuation that would be achieved with overall losses of
only 50% and 0%.

The heralded single photon case $\hat{\rho}_{1}$ in Fig. 3(c) displays
essentially opposite behavior. The sub-Poissonian statistics $(Q<0)$ determine
a positive initial slope such that $K>1$, and this trend continues for all
values of $R$. Note that by our definition of relative attenuation, this does
not indicate $\langle\hat{n}\rangle_{out}>\langle\hat{n}\rangle_{in}$, and is
much different from the “photon excess” observed for super-Poissonian states
after SPS. The state is still attenuated relative to the input, and this can
be seen by comparing the observed $K$ to the theoretically predicted values
for ideal Fock states in Fig. 3(d), for which
$\langle\hat{n}\rangle_{out}=\langle\hat{n}\rangle_{in}$. In contrast to the
states $\ket{\alpha}$ and $\ket{\xi}$, however, $K>1$ indicates that the
degree of heralded attenuation by ZPS is weaker than that of ordinary
attenuation. With previously determined values of $\kappa_{f}\approx 0.86$ and
$\eta_{1}\approx 0.32$, which combine to give an overall loss of 73%, a fit of
the data in panel (c) to Eq. 10 indicates a single photon probability of
$\beta\approx 0.38$ for our initial state. Two theoretical curves with
$\beta=0.38$ and improved overall losses of 50% and 0% show more extreme
deviations from $K=1$.

Interestingly, the theoretical curves in Figure 3(b)-(d) show as detector
inefficiency and losses increase, the ZPS statistics for both sub- and super-
Poissonian states converge toward the Poissonian case $K=1$. Losses before or
after the VBS play an identical role in reshaping $K(R)$ in Eqs. 8-10, and so
we can explain this in two ways. First, the photon number distributions of
$\ket{\xi}$ and $\hat{\rho}_{1}$ tend toward Poissonian statistics after
ordinary attenuation, i.e., loss [30]. As losses before the beamsplitter
increase, the observed $K(R)$ values should therefore tend to unity and
resemble those of the coherent state. Alternatively, losses before the
heralding detector $D_{1}$ reduce our ability to distinguish “true” vacuum in
the reflected mode from one or more photons [7]. As effective efficiency
decreases, “no-click” events herald a mixture of the desired ZPS state with
unwanted noise, becoming identical to ordinary attenuation in the zero-
efficiency limit.

Figure 4: Relative attenuation $K$ for the states $\ket{\xi}$ and
$\hat{\rho}_{1}$ at $R\approx 0.5$, measured as the heralding detection
channel $D_{1}$ is moved out of the mode of interest by a distance $\Delta x$,
which modulates heralding efficiency $\eta_{1}$. As $\eta_{1}$ decreases down
to zero (right axis), the relative attenuation values (left axis) converge to
$K=1$ as ZPS becomes increasingly ineffective. Here, heralding efficiency is
normalized to its maximum value of $\eta_{1}\approx 0.32$.

Figure 4 demonstrates this effect by studying $K(R)$ as detector efficiency
$\eta_{1}$ decreases. Here, the multimode fiber launcher coupled to heralding
detector $D_{1}$ is scanned a distance $\Delta x$ out of the mode of interest
within the U-bench (see inset and Fig. 2). The blue curve shows the degree of
spatial mode overlap as the fiber is moved, measured separately as the
fraction of power coupled from an auxiliary source. As detector mode overlap
decreases, it can be seen that the relative attenuations $K^{(\xi)}$ and
$K^{(\hat{\rho}_{1})}$ (measured at $R\approx 0.5$) converge to the Poissonian
value of $K=1$. Although we continue to herald on zero photons in $D_{1}$,
these measurements carry less information about the system and ZPS becomes
increasingly ineffective.

The results presented in Figs. 3 and 4 show a strong connection between the
effects of ZPS and the photon number distributions of initial states.
Importantly, however, measuring $K(R)$ yields no information about phase or
coherence between number states. Consequently, as has been argued for the
“photon excess” after SPS [5, 6], the effects of ZPS observed here can also be
replicated using counting statistics of classical particles and probabilistic
subtraction. For ZPS, probability of success drops exponentially with higher
numbers of photons, thus shifting the mean of the number distribution
downward. Additionally, the “photon excess” or “photon deficit” exhibited by
thermal light undergoing SPS or ZPS can also be understood as intensity
fluctuations of the classical electromagnetic field (i.e., correlated
intensities at detectors $D_{1}$ and $D_{2}$).

In this sense, the “photon deficit” observed in ZPS, much like the “photon
excess” of SPS, is not a purely quantum mechanical effect. Nonetheless, it is
important to note that the degree of attenuation measured after ZPS can reveal
existing nonclassicality of the input states. For example, our observations
confirmed the nonclassical [20] sub-Poissonian statistics of the heralded
single photon $\hat{\rho}_{1}$. Furthermore, quantum state tomography or some
other phase-dependent measurement would reveal that ZPS not only attenuates
quantum states, but does so noiselessly (i.e., preserves their coherence)
[31]. This property is exactly what makes ZPS promising for applications in
quantum communication [9, 10, 11].

##  V Conclusions

In summary, we have experimentally demonstrated that the zero-photon
subtraction (ZPS) process of Figure 1 can reduce the mean photon number of
quantum optical states, despite no photons being removed from the system. Our
experiment tested the effects of ZPS on three unique classes of input states
(sub-, super- and Poissonian) using a beamsplitter with variable reflectance
$R$. By studying the relative attenuation ratio $K$ as a function of $R$, the
observed trends reveal a connection to Mandel’s $Q$-parameter in the regime of
$R\ll 1$. More precisely, the initial slope of $K(R)$ near $R=0$ is equal to
$-Q_{in}$, resulting in distinct behavior for each input state. Consequently,
(sub-) super-Poissonian states will be attenuated (less) more by ZPS than by
ordinary attenuation with a weakly reflecting beamsplitter. These ZPS effects
are complementary to the effects of typical single-photon subtraction (SPS) in
the same regime $R\ll 1$. Most notably, super-Poissonian states that exhibit a
“photon excess” after SPS will also exhibit a unique “photon deficit” after
ZPS.

These observations were made possible by actively heralding on the detection
of zero photons with a single-photon detector [7]. We further confirmed the
need for high efficiency in the heralding mode by measuring the convergence of
non-Poissonian attenuation $K$ to the benchmark Poissonian case as losses
increased.

Although not revealed by the photon counting measurements reported here, ZPS
can preserve the coherence of quantum states [31], making it useful for
quantum communications as a noiseless attenuator [9, 10, 11]. Our results
provide further insight into the nature of this transformation and its
relationship to other techniques in quantum state engineering by conditional
measurements.

###### Acknowledgements.

We would like to thank S. U. Shringarpure for many valuable discussions on
this topic. This work was supported by the National Science Foundation under
Grant No. PHY-2013464.

*

## Appendix A

This Appendix provides further calculations regarding ZPS with realistic
single-photon detectors.

To account for imperfect heralding on zero with finite efficiency $\eta_{1}$,
we can derive an alternate version of Eq. 1:

$\displaystyle\langle\hat{n}\rangle_{out}$
$\displaystyle=\frac{\text{tr}\left\\{\hat{B}\hat{\rho}\hat{B}^{\dagger}\hat{\Pi}^{(NC)}_{1}\hat{n}_{2}\right\\}}{\text{tr}\left\\{\hat{B}\hat{\rho}\hat{B}^{\dagger}\hat{\Pi}^{(NC)}_{1}\right\\}}$
(11)
$\displaystyle=\frac{1-R}{1-R\eta_{1}}\frac{\sum_{n}n\rho_{nn}(1-R\eta_{1})^{n}}{\sum_{n}\rho_{nn}(1-R\eta_{1})^{n}}$

where the subscripts $1$ and $2$ indicate detection channels $D_{1}$ and
$D_{2}$, $\hat{B}$ is the unitary beam-splitter operator [20], and we have
used the standard POVMs for non-PNR detectors for “click” $(C)$ and “no-click”
$(NC)$ events [32]:

$\displaystyle\hat{\Pi}^{(C)}_{i}$ $\displaystyle=\mathbbm{1}-\Pi^{(NC)}_{i}$
$\displaystyle\hat{\Pi}^{(NC)}_{i}$
$\displaystyle=\sum_{n}(1-\eta_{i})^{n}\ket{n}\bra{n}$ (12)

for $i=1,2$. Substitution into the definition of $K$ (Eq. 3) yields the same
expression, but with the replacement $R\rightarrow R\eta_{1}$.

To account for losses before the attenuator, we introduce a preceding
beamsplitter with transmittance $\kappa$, where the reflected mode is lost to
the environment. Tracing over the outputs of both beamsplitters as in 11, we
once again find an expression where losses are included with a simple
replacement:

$\displaystyle K_{exp}(R)=K_{ideal}(\kappa R\eta_{1})$ (13)

The above relationship gives us Eqs. 8-10. Mathematically, this means the
experimental $K$ and the ideal lossless case have the same general behavior as
a function of $R$, except the effective domain is limited from $R\in[0,1)$ to
$R_{exp}\in[0,\kappa\eta_{1})$.

Finally, we can also account for the lack of PNR capability in the output
detector by replacing $\hat{n}_{2}$ in Eq. 11 with $\hat{\Pi}^{(C)}_{2}$. The
effective relative attenuation becomes:

$K_{\text{click}}(R)=\frac{\sum_{n}\rho_{nn}\left[(1-\kappa
R\eta_{1})^{n}-(1-\kappa
R\eta_{1}-T\eta_{2})^{n}\right]}{\left[\sum_{n}\rho_{nn}(1-\kappa
R\eta_{1})^{n}\right]\left[1-\sum_{n}\rho_{nn}(1-T\eta_{2})^{n}\right]}$ (14)

By taking the low-efficiency limit $\eta_{2}\rightarrow 0$ and applying
L’Hopital’s rule once, we find that $K_{\text{click}}(R)$ converges to $K(R)$.
For low mean photon numbers, this converges quickly enough to make the
approximation $K_{\text{click}}\approx K$, and the difference is negligible
for our experimental conditions.

## References

  * Wenger  _et al._ [2004] J. Wenger, R. Tualle-Brouri, and P. Grangier, Non-gaussian statistics from individual pulses of squeezed light, Phys. Rev. Lett. 92, 153601 (2004).
  * Ourjoumtsev  _et al._ [2006] A. Ourjoumtsev, R. Tualle-Brouri, J. Laurat, and P. Grangier, Generating optical schrödinger kittens for quantum information processing, Science 312, 83 (2006).
  * Ueda  _et al._ [1990] M. Ueda, N. Imoto, and T. Ogawa, Quantum theory for continuous photodetection processes, Phys. Rev. A 41, 3891 (1990).
  * Mizrahi and Dodonov [2002] S. S. Mizrahi and V. V. Dodonov, Creating quanta with an annihilation operator, J. Phys. A 35, 8847 (2002).
  * Parigi  _et al._ [2007] V. Parigi, A. Zavatta, M. Kim, and M. Bellini, Probing quantum commutation rules by addition and subtraction of single photons to/from a light field, Science 317, 1890 (2007).
  * Zavatta  _et al._ [2008] A. Zavatta, V. Parigi, M. S. Kim, and M. Bellini, Subtracting photons from arbitrary light fields: experimental test of coherent state invariance by single-photon annihilation,  10, 123006 (2008).
  * Nunn  _et al._ [2021] C. M. Nunn, J. D. Franson, and T. B. Pittman, Heralding on the detection of zero photons, Phys. Rev. A 104, 033717 (2021).
  * Gagatsos  _et al._ [2014] C. N. Gagatsos, J. Fiurášek, A. Zavatta, M. Bellini, and N. J. Cerf, Heralded noiseless amplification and attenuation of non-gaussian states of light, Phys. Rev. A 89, 062311 (2014).
  * Mičuda  _et al._ [2012] M. Mičuda, I. Straka, M. Miková, M. Dušek, N. J. Cerf, J. Fiurášek, and M. Ježek, Noiseless loss suppression in quantum optical communication, Phys. Rev. Lett. 109, 180503 (2012).
  * Brewster  _et al._ [2017] R. A. Brewster, I. C. Nodurft, T. B. Pittman, and J. D. Franson, Noiseless attenuation using an optical parametric amplifier, Phys. Rev. A 96, 042307 (2017).
  * Guo  _et al._ [2020] Y. Guo, J. Ding, Y. Mao, W. Ye, Q. Liao, and D. Huang, Quantum catalysis-based discrete modulation continuous variable quantum key distribution with eight states, Phys. Lett. A 384, 126340 (2020).
  * Bellini and Zavatta [2010] M. Bellini and A. Zavatta, Chapter 2 - manipulating light states by single-photon addition and subtraction, in  _Progress in Optics_, Progress in Optics, Vol. 55, edited by E. Wolf (Elsevier, 2010) pp. 41–83.
  * Allevi  _et al._ [2010] A. Allevi, A. Andreoni, M. Bondani, M. G. Genoni, and S. Olivares, Reliable source of conditional states from single-mode pulsed thermal fields by multiple-photon subtraction, Phys. Rev. A 82, 013816 (2010).
  * Zhai  _et al._ [2013] Y. Zhai, F. E. Becerra, B. L. Glebov, J. Wen, A. E. Lita, B. Calkins, T. Gerrits, J. Fan, S. W. Nam, and A. Migdall, Photon-number-resolved detection of photon-subtracted thermal light, Opt. Lett. 38, 2171 (2013).
  * Vidrighin  _et al._ [2016] M. D. Vidrighin, O. Dahlsten, M. Barbieri, M. S. Kim, V. Vedral, and I. A. Walmsley, Photonic maxwell’s demon, Phys. Rev. Lett. 116, 050401 (2016).
  * Bogdanov  _et al._ [2017] Y. I. Bogdanov, K. G. Katamadze, G. V. Avosopiants, L. V. Belinsky, N. A. Bogdanova, A. A. Kalinkin, and S. P. Kulik, Multiphoton subtracted thermal states: Description, preparation, and reconstruction, Phys. Rev. A 96, 063803 (2017).
  * Hloušek  _et al._ [2017] J. Hloušek, M. Ježek, and R. Filip, Work and information from thermal states after subtraction of energy quanta, Sci. Rep. 7, 13046 (2017).
  * Magaña-Loaiza  _et al._ [2019] O. S. Magaña-Loaiza, R. d. J. León-Montiel, A. Perez-Leija, A. B. U’Ren, C. You, K. Busch, A. E. Lita, S. W. Nam, R. P. Mirin, and T. Gerrits, Multiphoton quantum-state engineering using conditional measurements, npj Quantum Inf. 5, 80 (2019).
  * Katamadze  _et al._ [2020] K. G. Katamadze, G. V. Avosopiants, N. A. Bogdanova, Y. I. Bogdanov, and S. P. Kulik, Multimode thermal states with multiphoton subtraction: Study of the photon-number distribution in the selected subsystem, Phys. Rev. A 101, 013811 (2020).
  * Gerry and Knight [2004] C. C. Gerry and P. L. Knight,  _Introductory quantum optics_ (Cambridge University Press, Cambridge, 2004).
  * Digonnet and Shaw [1982] M. Digonnet and H. Shaw, Analysis of a tunable single mode optical fiber coupler, IEEE Transactions on Microwave Theory and Techniques 30, 592 (1982).
  * Mandel [1979] L. Mandel, Sub-poissonian photon statistics in resonance fluorescence, Opt. Lett. 4, 205 (1979).
  * Note [1] Historically, “photon subtraction” refers to application of the annihilation operator, predating its realization by conditional measurements in the output of a beamsplitter (BS) [12]. However, the term often encompasses the more general BS transformations which relax the weak reflectance requirement (e.g., [33, 34]), which we follow here in our discussion of ZPS.
  * Hong  _et al._ [1987] C. K. Hong, Z. Y. Ou, and L. Mandel, Measurement of subpicosecond time intervals between two photons by interference, Phys. Rev. Lett. 59, 2044 (1987).
  * Kim  _et al._ [2002] M. S. Kim, W. Son, V. Bužek, and P. L. Knight, Entanglement by a beam splitter: Nonclassicality as a prerequisite for entanglement, Phys. Rev. A 65, 032323 (2002).
  * Pittman  _et al._ [2005] T. B. Pittman, B. C. Jacobs, and J. D. Franson, Heralding single photons from pulsed parametric down-conversion, Opt. Commun. 246, 545 (2005).
  * Bovino  _et al._ [2003] F. A. Bovino, P. Varisco, A. Maria Colla, G. Castagnoli, G. Di Giuseppe, and A. V. Sergienko, Effective fiber-coupling of entangled photons for quantum communication, Opt. Commun. 227, 343 (2003).
  * Wahl  _et al._ [2020] M. Wahl, T. Röhlicke, S. Kulisch, S. Rohilla, B. Krämer, and A. C. Hocke, Photon arrival time tagging with many channels, sub-nanosecond deadtime, very high throughput, and fiber optic remote synchronization, Rev. Sci. Inst. 91, 013108 (2020).
  * Klyshko [1980] D. N. Klyshko, Use of two-photon light for absolute calibration of photoelectric detectors, Sov. J. Quantum Electron. 10, 1112 (1980).
  * Hu  _et al._ [2007] Y. Hu, X. Peng, T. Li, and H. Guo, On the poisson approximation to photon distribution for faint lasers, Phys. Lett. A 367, 173 (2007).
  * Shringarpure  _et al._ [2021] S. U. Shringarpure, C. M. Nunn, T. B. Pittman, and J. D. Franson, Coherence of quantum states after noiseless attenuation (2021), arXiv:2107.04582 .
  * Stevens [2013] M. J. Stevens, Chapter 2 - photon statistics, measurements, and measurements tools, in  _Single-Photon Generation and Detection_, Experimental Methods in the Physical Sciences, Vol. 45, edited by A. Migdall, S. V. Polyakov, J. Fan, and J. C. Bienfang (Academic Press, 2013) pp. 25–68.
  * Olivares  _et al._ [2003] S. Olivares, M. G. A. Paris, and R. Bonifacio, Teleportation improvement by inconclusive photon subtraction, Phys. Rev. A 67, 032314 (2003).
  * Kim  _et al._ [2005] M. S. Kim, E. Park, P. L. Knight, and H. Jeong, Nonclassicality of a photon-subtracted gaussian field, Phys. Rev. A 71, 043805 (2005).