# Emergence of the Born rule in quantum optics

Brian R. La Cour  Applied Research Laboratories, The University of Texas at
Austin, P.O. Box 8029, Austin, TX 78713-8029<EMAIL_ADDRESS>Morgan
C. Williamson

###### Abstract

The Born rule provides a fundamental connection between theory and observation
in quantum mechanics, yet its origin remains a mystery. We consider this
problem within the context of quantum optics using only classical physics and
the assumption of a quantum electrodynamic vacuum that is real rather than
virtual. The connection to observation is made via classical intensity
threshold detectors that are used as a simple, deterministic model of photon
detection. By following standard experimental conventions of data analysis on
discrete detection events, we show that this model is capable of reproducing
several observed phenomena thought to be uniquely quantum in nature, thus
providing greater elucidation of the quantum-classical boundary.

##  1 Introduction

Since the appearance of Bell’s inequality, it has become apparent that local
hidden variable models cannot be compatible with the complete mathematical
formalism of quantum mechanics [1, 2, 3, 4]. Indeed, recent loophole-free
experiments appear to be consistent with this conclusion [5, 6, 7, 8].
Nevertheless, there remains the open question of which observed phenomena, in
particular, are truly quantum in nature and have no classical analogue. This
question of elucidating the quantum-classical boundary is of practical
importance, as many new and emerging technologies, such as quantum computing,
quantum communication, and quantum sensing, rely upon this distinction for
their efficacy and security [9].

The field of quantum optics would seem to be a good place to explore this
question, as the systems of interest are relatively simple to describe in
terms of discrete field modes, while the important light-matter interactions
may be restricted to the physics of photodetection devices. One of the more
curious aspects of quantum optics is the concept of the vacuum or zero-point
field (ZPF). In quantum electrodynamics (QED), a vacuum state is defined
simply to be the lowest energy state of a given field mode [10]. The number of
photons in this state is taken to be zero, yet its energy is nonzero, giving
rise to the notion of “virtual” photons. Although the quantum vacuum is viewed
as being only virtual, its effects are quite real. Phenomena such as the
Casimir force, van der Waals attraction, Lamb shifts, and spontaneous emission
are all believed to have their origin in the quantum vacuum [11].

The prominence of vacuum states in quantum optics suggests that they may be
useful in developing a physical theory that explores the quantum-classical
boundary. In this work, we shall proceed by supposing that the quantum vacuum
of QED is real, not virtual. In doing so, we shall abandon all formal
reference to quantum theory and consider a world governed solely by classical
physics, albeit one in which the presence of a reified vacuum field is
unavoidable. Our connection to quantum theory will lie solely in the demand
that the statistical description of the real vacuum field match that of the
virtual one. Our goal in doing so will be to explore which observed quantum
phenomena can be explained under this supposition. In particular, we shall
explore in this work the emergence of the Born rule as a statistical
prediction that is applicable only within a certain regime of validity and
application.

Several previous attempts have been made to derive the Born rule from first
principles [12]. Max Born, in his original 1926 paper, considered the problem
of perturbative scattering and suggested that the resulting energy may be
interpreted as a statistical average if the scattering amplitudes, when
properly squared, are interpreted as probabilities [13]. Gleason provided the
first attempt at a mathematical derivation of the Born rule but relied on an
assumed association of Hermitian operators with measurement observables [14].
David Deutsch, in 1999, went further to argue that elementary decision theory
may be used to deduce the Born rule as a necessary consequence of the other
quantum axioms [15]. This argument has since been criticized to be circular,
as it requires the assumption of an agent with a particular predilection for
$L_{2}$ norms [16]. Zurek has suggested decoherence as an explanation of the
Born rule [17], although this view has been criticized as well to be
insufficient [18]. More recently, Masanes et al. have claimed to derive the
Born rule by assuming, among other things, that measurements consist of well-
defined trials and always produce one of a pre-defined set of outcomes [19].
While seemingly innocuous, this assumption does not always hold in real,
experimental settings where, for example, photons are detected at random times
or, more often, not at all. An interesting result from Allahverdyan et al.
provides a derivation of the Born rule from the dynamical law of quantum
mechanics with the context of spin systems [20].

Working within the confines of the formalism does not seem a promising
approach to deriving physical laws. What these and other attempts to derive
the Born rule lack is any attempt to model the actual physics of measurement.
This paper seeks to address that point by considering a deterministic model of
measurement together with a reified quantum vacuum.

A reified quantum vacuum is the premise behind the theory of stochastic
electrodynamics (SED), and we adopt a similar outlook here [21]. Previous work
in SED considered the statistical behavior of physical systems immersed in the
zero-point field. These included classical descriptions of the quantum
harmonic oscillator ground state as well as spontaneous parametric
downconversion [22]. Although these efforts were successful insofar as they
predicted probability density functions identical to the corresponding quantum
Wigner function, they failed to fully appreciate the critical role of
measurement and experimental procedure in the observation of quantum
phenomena. In particular, the role of post-selection and its relation to
contextuality has received little attention within SED.

To address this deficiency, we shall consider here a local, deterministic
model of photon detection wherein the only random variables determining the
outcome of a measurement are those associated with the relevant vacuum states
incident upon the device. This approach differs from previous work in
stochastic optics, an offshoot of SED focused on quantum optics, wherein the
intensity of incident waves above a given threshold determines only the
_probability_ of an outcome, leaving the actual realization to be determined
by yet another, implicit, hidden variable [23]. Our approach uses a
deterministic amplitude threshold crossing scheme to define detection events
and is similar to the work of other researchers in this regard [24, 25, 26]. A
key difference from previous work is the use of post-selection and the
examination of asymptotic behavior to approximate ideal quantum predictions.

The structure of the paper is as follows. In section 2 we describe the
mathematical model used to describe the reified vacuum field and use the
single-mode approximation to make the correct correspondence with quantum
optics. The connection to observation is made in section 3, where we describe
a deterministic model of quantum measurement using amplitude threshold
detection. From this, the Born rule is shown to arise as an emergent and
approximate property of the model in the presence of measurements. Finally, in
section 4 we consider the general problem of transformations of multiple
vacuum modes under linear optics to arrive at a model approximating single-
photon, multi-mode quantum states. Conclusions are summarized in section 5.

##  2 The Reified Vacuum Field

###  2.1 Continuum Description

Any classical electric field may be written in terms of a continuum of plane
wave modes. Thus, the electric field at a point $\boldsymbol{x}$ and time $t$
may be written, in Gaussian units, as

$\boldsymbol{E}(\boldsymbol{x},t)=\frac{1}{2\pi}\int\sqrt{\mathcal{E}(\boldsymbol{k})}\,\boldsymbol{a}(\boldsymbol{k})\,e^{i\boldsymbol{k}\cdot\boldsymbol{x}-i\omega(\boldsymbol{k})t}\,d\boldsymbol{k}+\mathrm{c.c.}\;,$
(1)

where $\mathcal{E}(\boldsymbol{k})\geq 0$ is a scale factor related to the
modal energy for wave vector $\boldsymbol{k}\in\mathbb{R}^{3}$,
$\boldsymbol{a}(\boldsymbol{k})\in\mathbb{C}^{3}$ gives the field direction
and phase, and $\omega(\boldsymbol{k})\geq 0$ is the angular frequency. For a
classical vacuum, $\omega(\boldsymbol{k})=\|\boldsymbol{k}\|c$, where
$\|\boldsymbol{k}\|$ is the magnitude of $\boldsymbol{k}$ and $c$ is the speed
of light. The term “c.c.” represents the complex conjugate of the term to the
left. The magnetic field is similarly described, with
$\boldsymbol{a}(\boldsymbol{k})$ replaced by
$\boldsymbol{k}\times\boldsymbol{a}(\boldsymbol{k})$, so that specifying
$\mathcal{E}(\boldsymbol{k})$, $\boldsymbol{a}(\boldsymbol{k})$, and
$\omega(\boldsymbol{k})$ for all $\boldsymbol{k}\in\mathbb{R}^{3}$ provides a
complete description of the electromagnetic field. Without loss of generality,
we shall take $\boldsymbol{a}(\boldsymbol{k})$ to be stochastic, while
$\mathcal{E}(\boldsymbol{k})$ is assumed fixed.

For convenience, we may decompose $\boldsymbol{a}(\boldsymbol{k})$ into
orthogonal polarization modes. For each wave vector $\boldsymbol{k}$, let
$\hat{\boldsymbol{e}}_{0}(\boldsymbol{k})$ and
$\hat{\boldsymbol{e}}_{1}(\boldsymbol{k})$ be any two orthogonal polarization
vectors (i.e., complex vectors such that
$\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{e}}_{\mu}(\boldsymbol{k})=0$, where
$\hat{\boldsymbol{k}}=\boldsymbol{k}/\|\boldsymbol{k}\|$, and
$\hat{\boldsymbol{e}}_{\mu}(\boldsymbol{k})^{*}\cdot\hat{\boldsymbol{e}}_{\nu}(\boldsymbol{k})=\delta_{\mu,\nu}$
for all $\mu,\nu\in\\{0,1\\}$). We may then write

$\boldsymbol{a}(\boldsymbol{k})=a_{0}(\boldsymbol{k})\,\hat{\boldsymbol{e}}_{0}(\boldsymbol{k})+a_{1}(\boldsymbol{k})\,\hat{\boldsymbol{e}}_{1}(\boldsymbol{k})\;,$
(2)

where
$a_{\mu}(\boldsymbol{k})=\hat{\boldsymbol{e}}_{\mu}(\boldsymbol{k})^{*}\cdot\boldsymbol{a}(\boldsymbol{k})\in\mathbb{C}$.
Note that the choice of polarization vectors is arbitrary, and may vary with
$\boldsymbol{k}$, but is otherwise taken to be fixed and nonrandom.

We now turn to the correspondence with quantum theory. Consistency with
quantum electrodynamics will require that, at zero temperature,

$\mathcal{E}(\boldsymbol{k})=\mathcal{E}_{0}(\omega(\boldsymbol{k})):=\tfrac{1}{2}\hbar\omega(\boldsymbol{k})\;,$
(3)

where we have now introduced $\hbar$, Planck’s constant divided by $2\pi$, as
setting the fundamental scale of the vacuum field. For nonzero temperatures,
$\mathcal{E}_{0}(\omega)$ is replaced by the expression

$\mathcal{E}_{T}(\omega)=\tfrac{1}{2}\hbar\omega+\frac{\hbar\omega}{e^{\hbar\omega/k_{B}T}-1}\;,$
(4)

where $\omega>0$, $k_{B}$ is Boltzmann’s constant, and $T>0$ is the absolute
temperature. Since the density of states is given by
$\omega^{2}/(\pi^{2}c^{3})$, the spectral energy density is

$\rho_{T}(\omega)=\frac{\hbar\omega^{3}}{\pi^{2}c^{3}}\left(\tfrac{1}{2}+\frac{1}{e^{\hbar\omega/k_{B}T}-1}\right)\;,$
(5)

which corresponds to Planck’s “second quantum theory” of blackbody radiation,
with a zero-point energy term included [27]. Note also that, at zero
temperature, $\rho_{0}(\omega)=\hbar\omega^{3}/(2\pi^{2}c^{3})$ is Lorentz
invariant, owing to the cubic dependence on frequency, so the spectral energy
density is the same in all inertial reference frames [28].

The stochastic nature of the field is described entirely by
$\boldsymbol{a}(\cdot)$, and consistency with QED requires that it be a
complex Gaussian random vector field such that, for any choice of polarization
vectors, $\mathsf{E}[a_{\mu}(\boldsymbol{k})]=0$ and

$\displaystyle\mathsf{E}[a_{\mu}(\boldsymbol{k})\,a_{\nu}(\boldsymbol{k}^{\prime})^{*}]$
$\displaystyle=\delta_{\mu,\nu}\,\delta(\boldsymbol{k}-\boldsymbol{k}^{\prime})$
(6a)
$\displaystyle\mathsf{E}[a_{\mu}(\boldsymbol{k})\,a_{\nu}(\boldsymbol{k}^{\prime})]$
$\displaystyle=0\;,$ (6b)

where $\mathsf{E}[\cdot]$ denotes an expectation value [29]. More generally,
the $n$-point correlations of the field are given by

$\mathsf{E}\left[\prod_{i=1}^{n}a_{\mu_{i}}(\boldsymbol{k}_{i})\,a_{\nu_{i}}(\boldsymbol{k}_{i}^{\prime})^{*}\right]=\prod_{i=1}^{n}\prod_{j=i}^{n}\delta_{\mu_{i},\nu_{j}}\,\delta(\boldsymbol{k}_{i}-\boldsymbol{k}_{j}^{\prime})\;,$
(7)

with all other combinations giving a zero expectation value. Of course, this
mathematical correlation structure is only an idealization; on _some_ spatio-
temporal scale, the field must surely be correlated. We would furthermore
expect that the statistical character of the field, its scale and
correlations, might also change over time and space. Nevertheless, we shall
proceed with this modest idealization of the zero-point field, as it will
provide a useful model for the quantum vacuum.

###  2.2 Discrete-Mode Approximation

One can approximate the continuum of wave vector modes by a set of closely
spaced discrete modes in a notional box. Given a cube of length $L>0$, we
define a set $K$ of discrete-mode wave vectors as follows:

$K=\left\\{\frac{2\pi n_{1}}{L}\hat{\boldsymbol{x}}+\frac{2\pi
n_{2}}{L}\hat{\boldsymbol{y}}+\frac{2\pi
n_{3}}{L}\hat{\boldsymbol{z}}\;:\;n_{1},n_{2},n_{3}\in\mathbb{Z}\right\\}\;.$
(8)

The continuum wave vector space may now be decomposed into notional discrete
cells $C(\boldsymbol{k})=\boldsymbol{k}+[-\frac{\pi}{L},\frac{\pi}{L})^{3}$,
each centered on a wave vector $\boldsymbol{k}\in K$. Since the cells are
disjoint and their union comprises all of $\mathbb{R}^{3}$, we may rewrite
equation (1) as follows:

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\boldsymbol{E}(\boldsymbol{x},t)=\frac{1}{2\pi}\sum_{\boldsymbol{k}\in
K}\sum_{\mu}\int_{C(\boldsymbol{k})}\sqrt{\mathcal{E}_{0}(\omega(\boldsymbol{k}^{\prime}))}\,a_{\mu}(\boldsymbol{k}^{\prime})\,\hat{\boldsymbol{e}}_{\mu}(\boldsymbol{k}^{\prime})\,e^{i\boldsymbol{k}^{\prime}\cdot\boldsymbol{x}-i\omega(\boldsymbol{k}^{\prime})t}\,d\boldsymbol{k}^{\prime}+\mathrm{c.c.}$
(9)

Furthermore, if the cells are small (i.e., $L$ is large), we may make the
approximation

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\boldsymbol{E}(\boldsymbol{x},t)\approx\frac{1}{2\pi}\sum_{\boldsymbol{k}\in
K}\sum_{\mu}\sqrt{\mathcal{E}_{0}(\omega(\boldsymbol{k}))}\,\hat{\boldsymbol{e}}_{\mu}(\boldsymbol{k})\,e^{i\boldsymbol{k}\cdot\boldsymbol{x}-i\omega(\boldsymbol{k})t}\int_{C(\boldsymbol{k})}a_{\mu}(\boldsymbol{k}^{\prime})\,d\boldsymbol{k}^{\prime}+\mathrm{c.c.}$
(10)

This last integral yields a complex Gaussian random variable with zero mean
and a variance of $\Delta\boldsymbol{k}=(2\pi/L)^{3}=8\pi^{3}/V$ corresponding
to the volume of each cell [30]. We may therefore write

$\int_{C(\boldsymbol{k})}a_{\mu}(\boldsymbol{k}^{\prime})\,d\boldsymbol{k}^{\prime}=z_{\mu,\boldsymbol{k}}\sqrt{\Delta\boldsymbol{k}}\;,$
(11)

where $z_{\mu,\boldsymbol{k}}$ is a standard complex Gaussian random variable
(i.e., a complex Gaussian random variable such that
$\mathsf{E}[z_{\mu,\boldsymbol{k}}]=0$,
$\mathsf{E}[|z_{\mu,\boldsymbol{k}}|^{2}]=1$, and
$\mathsf{E}[z_{\mu,\boldsymbol{k}}^{2}]=0$). Equivalently, we may write
$z_{\mu,\boldsymbol{k}}$ in the form $z_{\mu,\boldsymbol{k}}=(x+iy)/\sqrt{2}$,
where $x,y$ are independent, real-valued standard normal random variables.

We note that all discrete modes differing in either wave vector or
polarization are independent since, for
$\boldsymbol{k},\boldsymbol{k}^{\prime}\in K$,

$\begin{split}\mathsf{E}[z_{\mu,\boldsymbol{k}}\,z_{\nu,\boldsymbol{k^{\prime}}}^{*}]&=\frac{1}{\Delta\boldsymbol{k}}\int_{C(\boldsymbol{k})}\int_{C(\boldsymbol{k}^{\prime})}\mathsf{E}[a_{\mu}(\boldsymbol{k}^{\prime\prime})^{*}a_{\nu}(\boldsymbol{k}^{\prime\prime\prime})]\,d\boldsymbol{k}^{\prime\prime}d\boldsymbol{k}^{\prime\prime\prime}\\\
&=\frac{1}{\Delta\boldsymbol{k}}\int_{C(\boldsymbol{k})}\int_{C(\boldsymbol{k}^{\prime})}\delta_{\mu,\nu}\,\delta(\boldsymbol{k}^{\prime\prime\prime}-\boldsymbol{k}^{\prime\prime})\,d\boldsymbol{k}^{\prime\prime}d\boldsymbol{k}^{\prime\prime\prime}\\\
&=\delta_{\mu,\nu}\,\frac{1}{\Delta\boldsymbol{k}}\int_{C(\boldsymbol{k})\cap
C(\boldsymbol{k}^{\prime})}d\boldsymbol{k}^{\prime\prime}\\\
&=\delta_{\mu,\nu}\,\delta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}\;.\end{split}$
(12)

We shall chiefly be concerned with descriptions in terms of the discrete-mode
approximation, as this affords the clearest correspondence with quantum
optics. In particular, the lowering operator $\hat{a}_{\mu,\boldsymbol{k}}$
for discrete mode $(\mu,\boldsymbol{k})$ may be associated with the random
variable $z_{\mu,\boldsymbol{k}}/\sqrt{2}$ in the sense that the vacuum
expectation of the symmetrized number operator equals the variance of the
corresponding random variable. To see this, observe that

$\bra{0}\tfrac{1}{2}\Bigl{(}\hat{a}_{\mu,\boldsymbol{k}}^{\dagger}\hat{a}_{\mu,\boldsymbol{k}}+\hat{a}_{\mu,\boldsymbol{k}}\hat{a}_{\mu,\boldsymbol{k}}^{\dagger}\Bigr{)}\ket{0}=\tfrac{1}{2}(0+1)=\tfrac{1}{2}\;,$
(13)

and, similarly,

$\mathsf{E}\left[\tfrac{1}{2}\left(\frac{z_{\mu,\boldsymbol{k}}^{*}}{\sqrt{2}}\frac{z_{\mu,\boldsymbol{k}}}{\sqrt{2}}+\frac{z_{\mu,\boldsymbol{k}}}{\sqrt{2}}\frac{z_{\mu,\boldsymbol{k}}^{*}}{\sqrt{2}}\right)\right]=\mathsf{E}\left[\left|\frac{z_{\mu}(\boldsymbol{k})}{\sqrt{2}}\right|^{2}\right]=\tfrac{1}{2}\;.$
(14)

Note that, since $z_{\mu,\boldsymbol{k}}$ and $z_{\mu,\boldsymbol{k}}^{*}$
commute, whereas $\hat{a}_{\mu,\boldsymbol{k}}$ and
$\hat{a}_{\mu,\boldsymbol{k}}^{\dagger}$ do not, symmetrization of the
operators is important to achieve the correct correspondence.

The connection to quantum optics can be further elucidated by examining the
modal energy. Quantum mechanically, the energy of the vacuum state is given by
the expectation value of the Hamilitonian
$\hat{H}=\frac{1}{2}(\hat{a}_{\mu,\boldsymbol{k}}^{\dagger}\hat{a}_{\mu,\boldsymbol{k}}+\hat{a}_{\mu,\boldsymbol{k}}\hat{a}_{\mu,\boldsymbol{k}}^{\dagger})\,\hbar\omega(\boldsymbol{k})$,
which is simply the symmetrized number operator scaled by
$\hbar\omega(\boldsymbol{k})$. Thus,
$\bra{0}\hat{H}\ket{0}=\frac{1}{2}\hbar\omega(\boldsymbol{k})$ is identified
as the average energy per vacuum mode.

To find the corresponding classical value, we begin by computing the energy
density of the electromagnetic field for the selected mode, as given by

$u(\boldsymbol{x},t)=\frac{1}{8\pi}\Bigl{[}\|\Delta\boldsymbol{E}(\boldsymbol{x},t)\|^{2}+\|\Delta\boldsymbol{B}(\boldsymbol{x},t)\|^{2}\Bigr{]}\;,$
(15)

where the single-mode electric field is

$\Delta\boldsymbol{E}(\boldsymbol{x},t)=\frac{1}{2\pi}\sqrt{\mathcal{E}_{0}(\omega)}\,z_{\mu,\boldsymbol{k}}\,\hat{\boldsymbol{e}}_{\mu}(\boldsymbol{k})\,e^{i\boldsymbol{k}\cdot\boldsymbol{x}-i\omega
t}\,\sqrt{\Delta\boldsymbol{k}}+\mathrm{c.c.}$ (16)

and the single-mode magnetic field is

$\Delta\boldsymbol{B}(\boldsymbol{x},t)=\frac{1}{2\pi}\sqrt{\mathcal{E}_{0}(\omega)}\,z_{\mu,\boldsymbol{k}}\,[\hat{\boldsymbol{k}}\times\hat{\boldsymbol{e}}_{\mu}(\boldsymbol{k})]\,e^{i\boldsymbol{k}\cdot\boldsymbol{x}-i\omega
t}\,\sqrt{\Delta\boldsymbol{k}}+\mathrm{c.c.}\;.$ (17)

Using the fact that, for any complex vector $\boldsymbol{v}$,
$\|\boldsymbol{v}+\boldsymbol{v}^{*}\|^{2}=2(\boldsymbol{v}^{*}\cdot\boldsymbol{v}+\mathrm{Re}[\boldsymbol{v}\cdot\boldsymbol{v}])$,
we find that

$\|\Delta\boldsymbol{E}(\boldsymbol{x},t)\|^{2}=\frac{2\pi\hbar\omega(\boldsymbol{k})}{V}\,\left(|z_{\mu,\boldsymbol{k}}|^{2}+\mathrm{Re}\left[z_{\mu,\boldsymbol{k}}^{2}e^{i2(\boldsymbol{k}\cdot\boldsymbol{x}-\omega(\boldsymbol{k})t)}\right]\right)\;,$
(18)

and, since
$\|\Delta\boldsymbol{E}(\boldsymbol{x},t)\|^{2}=\|\Delta\boldsymbol{B}(\boldsymbol{x},t)\|^{2}$,
we conclude that

$u(\boldsymbol{x},t)=\frac{\hbar\omega(\boldsymbol{k})}{2V}\,\left(|z_{\mu,\boldsymbol{k}}|^{2}+\mathrm{Re}\left[z_{\mu,\boldsymbol{k}}^{2}e^{i2(\boldsymbol{k}\cdot\boldsymbol{x}-\omega(\boldsymbol{k})t)}\right]\right)\;.$
(19)

Now consider the time average of $u(\boldsymbol{x},t)$, a spatially
independent random variable given by

$\bar{u}=\frac{\omega(\boldsymbol{k})}{2\pi}\int_{0}^{2\pi/\omega(\boldsymbol{k})}u(\boldsymbol{x},t)\,dt=\frac{\hbar\omega(\boldsymbol{k})}{2V}\,|z_{\mu,\boldsymbol{k}}|^{2}\;.$
(20)

The expectation value of this time average over realizations of the ZPF is
therefore

$\mathsf{E}[\bar{u}]=\tfrac{1}{2}\hbar\omega(\boldsymbol{k})/V\;.$ (21)

This result matches the quantum mechanical prediction if one integrates over a
box of volume $V$ to find the total expected energy. Of course, this volume is
only notional and arises as an artifact of our discrete-mode approximation. It
describes the degree to which the single-mode approximation is valid rather
than any physical volume. For, say, a conical beam with a small half-angle of
$\Delta\theta$ and a filtered bandwidth of $\Delta\omega$, we have
$\Delta\boldsymbol{k}=\pi\Delta\theta^{2}\Delta\omega/c$. Thus, as the beam is
narrowed, the notional volume increases and the energy density decreases
proportionally. An equivalent, and perhaps more physically meaningful,
interpretation of the quantum mechanical energy, then, might be that the
quantity $\bra{0}\hat{H}\ket{0}/V$ gives the expected energy density of a
single vacuum mode when the wave vector is filtered and collimated to a
resolution of $\Delta\boldsymbol{k}=8\pi^{3}/V$. In this context, we may there
identify a correspondence between the Hamiltonian operator $\hat{H}$ and the
time-averaged classical electromagnetic energy $H=\bar{u}V$.

###  2.3 Coherent States

In quantum optics, coherent states are considered the closest analogue to a
classical state. Previous work in SED has identified coherent states as
arising from, for example, classical driven harmonic oscillators coupled to
the ZPF [31]. Here we shall consider an optical analogue in which we add a
classical plane wave to a single mode of the ZPF.

Recall that, previously, we had defined the ZPF to be of the form

$\boldsymbol{E}(\boldsymbol{x},t)=\frac{1}{2\pi}\int\sqrt{\mathcal{E}_{0}(\omega(\boldsymbol{k}))}\,\boldsymbol{a}(\boldsymbol{k})\,e^{i\boldsymbol{k}\cdot\boldsymbol{x}-i\omega(\boldsymbol{k})t}\,d\boldsymbol{k}+\mathrm{c.c.}\;.$
(22)

We now add to this a plane wave with wave vector $\boldsymbol{k}_{0}$ and
polarization $\hat{\boldsymbol{e}}_{0}$ of the form

$\boldsymbol{F}(\boldsymbol{x},t)=E_{0}\,\hat{\boldsymbol{e}}_{0}\,e^{i\boldsymbol{k}_{0}\cdot\boldsymbol{x}-i\omega_{0}t}+E_{0}^{*}\,\hat{\boldsymbol{e}}_{0}^{*}\,e^{-i\boldsymbol{k}_{0}\cdot\boldsymbol{x}-i\omega_{0}t}\;,$
(23)

where $E_{0}\in\mathbb{C}$ is a complex number representing the amplitude and
phase of the external plane wave. The total electric field is now

$\boldsymbol{F}(\boldsymbol{x},t)+\boldsymbol{E}(\boldsymbol{x},t)\\\
=\int\left[E_{0}\,\hat{\boldsymbol{e}}_{0}\,\delta(\boldsymbol{k}-\boldsymbol{k}_{0})+\frac{1}{2\pi}\sqrt{\mathcal{E}_{0}(\omega(\boldsymbol{k}))}\,\boldsymbol{a}(\boldsymbol{k})\right]e^{i(\boldsymbol{k}\cdot\boldsymbol{x}-\omega
t)}\,d\boldsymbol{k}\ +\mathrm{c.c.}$ (24)

In the single-mode approximation with $\Delta\boldsymbol{k}=8\pi^{3}/V$, the
total field becomes

$\boldsymbol{F}(\boldsymbol{x},t)+\Delta\boldsymbol{E}(\boldsymbol{x},t)\approx\left[E_{0}+\frac{1}{2\pi}\sqrt{\mathcal{E}_{0}(\omega_{0})}\,z\,\sqrt{\Delta\boldsymbol{k}}\right]\hat{\boldsymbol{e}}_{0}\,e^{i(\boldsymbol{k}_{0}\cdot\boldsymbol{x}-\omega_{0}t)}+\mathrm{c.c.}\;$
(25)

where $z$ is a standard complex Gaussian random variable. For reasons that
will soon become apparent, we shall express $E_{0}$ in the form

$E_{0}=\alpha\sqrt{\frac{2\pi\hbar\omega_{0}}{V}}\;,$ (26)

where $\alpha\in\mathbb{C}$ is a complex number that will later be identified
as the coherent state parameter. The combined field in the single-mode
approximation may now be written

$\boldsymbol{F}(\boldsymbol{x},t)+\Delta\boldsymbol{E}(\boldsymbol{x},t)=\sqrt{\frac{2\pi\hbar\omega_{0}}{V}}\left(\alpha+\frac{z}{\sqrt{2}}\right)\hat{\boldsymbol{e}}_{0}\,e^{i(\boldsymbol{k}_{0}\cdot\boldsymbol{x}-\omega_{0}t)}+\mathrm{c.c.}$
(27)

The energy density of the corresponding electromagnetic field is

$\begin{split}u(\boldsymbol{x},t)&=\frac{1}{8\pi}\|\boldsymbol{F}(\boldsymbol{x},t)+\Delta\boldsymbol{E}(\boldsymbol{x},t)\|^{2}\times
2\\\
&=\frac{1}{4\pi}\left[\|\boldsymbol{F}(\boldsymbol{x},t)\|^{2}+\|\Delta\boldsymbol{E}(\boldsymbol{x},t)\|^{2}+2\boldsymbol{F}(\boldsymbol{x},t)\cdot\Delta\boldsymbol{E}(\boldsymbol{x},t)\right]\;,\end{split}$
(28)

where

$\|\boldsymbol{F}(\boldsymbol{x},t)\|^{2}=\frac{4\pi\hbar\omega_{0}}{V}\left(|\alpha|^{2}+\mathrm{Re}\left[\alpha^{2}e^{i2(\boldsymbol{k}_{0}\cdot\boldsymbol{x}-\omega_{0}t)}\right]\right)$
(29)
$\|\Delta\boldsymbol{E}(\boldsymbol{x},t)\|^{2}=\frac{4\pi\hbar\omega_{0}}{V}\,\frac{1}{2}\left(|z|^{2}+\mathrm{Re}\left[z^{2}e^{i2(\boldsymbol{k}_{0}\cdot\boldsymbol{x}-\omega_{0}t)}\right]\right)$
(30)
$\boldsymbol{F}\cdot\Delta\boldsymbol{E}(\boldsymbol{x},t)=\frac{4\pi\hbar\omega_{0}}{V}\,\frac{1}{\sqrt{2}}\mathrm{Re}\left[\alpha
z^{*}+\alpha
z\,e^{i2(\boldsymbol{k}_{0}\cdot\boldsymbol{x}-\omega_{0}t)}\right]\;.$ (31)

Taking the time average of $u(\boldsymbol{x},t)$ gives

$\bar{u}=\frac{\hbar\omega_{0}}{V}\left(|\alpha|^{2}+|z|^{2}+\sqrt{2}\mathrm{Re}[\alpha
z^{*}]\right)=\left|\alpha+\frac{z}{\sqrt{2}}\right|^{2}\frac{\hbar\omega_{0}}{V}\;,$
(32)

and the expectation value of $\bar{u}$ over realizations of the ZPF is

$\mathsf{E}[\bar{u}]=\left(|\alpha|^{2}+\tfrac{1}{2}\right)\frac{\hbar\omega_{0}}{V}\;.$
(33)

This result matches the familiar energy density
$\bra{\alpha}\hat{H}\ket{\alpha}/V$ of a quantum optical coherent state
$\ket{\alpha}$.

For general thermal states, $\mathcal{E}_{0}$ is replaced by
$\mathcal{E}_{T}$, as defined in equation (4), and, hence, $z/\sqrt{2}$ is
replaced by the scaled quantity $\sigma z$, where

$\sigma=\sqrt{\mathcal{E}_{T}(\omega_{0})/(\hbar\omega_{0})}\;.$ (34)

In this case, the average energy density becomes
$(|\alpha|^{2}+\sigma^{2})\hbar\omega_{0}/V$. Note that nonzero temperatures
merely have the effect of rescaling the ZPF for the given mode. At high
temperatures ($\sigma\gg|\alpha|$), the coherent state becomes
indistinguishable from thermal noise. Conversely, at large amplitudes
($|\alpha|\gg\sigma$), the coherent state becomes indistinguishable from a
classical plane wave of fixed amplitude and phase.

##  3 Amplitude Threshold Detection

We have described a mathematical model for the QED vacuum in terms of a
stochastic electromagnetic field. To make the important connection to
observation and discrete detection events, we now introduce a simple
deterministic model of photon detection based on amplitude threshold crossings
and motivated by the observed behavior of real detectors.

Suppose that we have, to arbitrary precision, isolated a single angular
frequency $\omega_{0}$, polarization mode $\hat{\boldsymbol{e}}_{0}$, and wave
vector mode $\boldsymbol{k}_{0}$ of the vacuum in the discrete-mode
approximation with wave vector resolution $\Delta\boldsymbol{k}$. For the
vacuum and coherent states, the energy density $u(\boldsymbol{x},t)$ at
position $\boldsymbol{x}$ and time $t$ varies sinusoidally in time and space.
We imagine a detection device that reacts slowly enough as to be sensitive
only to the time average, $\bar{u}$, of the energy density and note that this
averaging eliminates both the temporal and spatial dependence of the energy
density. Using a time average is justified by the fact that a typical period
of light is orders of magnitude shorter than the corresponding lag time for
the photoelectric effect [32].

Now, although $\bar{u}$ is constant across space and time, it varies from one
vacuum field realization to another due to the presence of the random variable
$z$. We now ask whether this time-averaged energy density falls above some
threshold $\Gamma^{2}\geq 0$. Such an outcome will be deemed a detection event
or “click” of a detector, and the probability of such an event occurring will
be denoted $\Pr[\bar{u}>\Gamma^{2}]$. Note that the vacuum realization $z$ is
the only source of randomness in determining this probability, and, in the
single-mode limit, the coherence time of the given vacuum mode is infinite.

###  3.1 Dark Counts

The time-averaged energy density of the vacuum is given by equation (32) with
$\alpha=0$ and, being the sum of two squared independent normal distributions,
has an exponential distribution with a mean of $\frac{1}{2}\hbar\omega_{0}/V$.
The probability of a detection event is therefore

$\Pr\left[\bar{u}>\Gamma^{2}\right]=\exp(-2V\Gamma^{2}/\hbar\omega_{0})\;.$
(35)

Since we have assumed $V$ is large, we may take $\Gamma^{2}$ to be comparably
small. In particular, we shall adopt the _single-mode limit_ , analogous to
the thermodynamic limit, in which

$\lim_{V\to\infty}\frac{V\Gamma^{2}}{\hbar\omega}=\gamma^{2}$ (36)

for some $\gamma\geq 0$. (Note that $\gamma$ may be specific to a particular
polarization, frequency and wave vector resolution.) In the single-mode limit,
the probability of a detection event is $\exp(-2\gamma^{2})$, which we
interpret as the probability of a dark count for the vacuum state at zero
temperature. In a thermal state ($T>0$), we replace $\frac{1}{2}\hbar\omega$
with $\sigma^{2}\hbar\omega$, so the probability of a dark count becomes
$\exp(-\gamma^{2}/\sigma^{2})$. This, again, becomes an effective rescaling of
the detection threshold, so there is no loss of generality in supposing $T=0$.

The prediction of a nonzero dark count rate at zero temperature is, strictly
speaking, at variance with quantum mechanical predictions. Even under ideal
conditions, our model predicts a nonzero probability of a vacuum detection
event; quantum mechanically this probability should be exactly zero. However,
even at extremely low temperatures, nonzero dark count rates are
experimentally observed [33].

For an explicit, albeit notional, example of a physical detection mechanism,
one may consider a classical charged particle in a bifurcating harmonic
potential. Such a potential has the quadratic form
$\frac{1}{2}m\omega^{2}x^{2}$ for mass $m$ and displacement $x$ for
$|x|\leq\ell$. For $|x|>\ell$, the potential is strongly repulsive and the
particle quickly accelerates away, thereby creating an observable event. Since
the trapped particle behaves as a high-Q linear filter, its behavior will
closely match that of the resonant vacuum mode. If the polarization is linear
and aligned with the displacement of the potential, the particle’s motion will
bifurcate and run away if the modal amplitude is sufficiently high.

Despite some similarities, the adoption of a threshold detection scheme for
modeling photon detection should not be construed as a semi-classical
treatment, as we are still completely within the confines of classical
physics. Although we have adopted a very simple model of single-photon
detection, these general qualitative observations are expected to hold in a
more detailed physical model. In what follows, we shall make no further
reference to the particular physical mechanism used for detection and will
instead focus on the more abstract notion of threshold detection in the
single-mode limit.

###  3.2 Emergence of the Born Rule

For coherent states, a detection event in the single-mode limit may be written

$\left|\alpha+\frac{z}{\sqrt{2}}\right|^{2}>\gamma^{2}\;.$ (37)

With detection events so defined, we may identify the complex amplitude $a$,
given by

$a=\alpha+\frac{z}{\sqrt{2}}\;,$ (38)

and note that $4|a|^{2}$ follows a non-central $\chi^{2}$ distribution with
two degrees of freedom and a noncentrality parameter of $4|\alpha|^{2}$. Thus,
the cumulative distribution function (cdf) of $|a|^{2}$ is given by the
expression [34]

$\Pr\Bigl{[}|a|^{2}\leq\gamma^{2}\Bigr{]}=1-Q_{1}\left(2|\alpha|,2\gamma\right)\;,$
(39)

where $Q_{1}(\cdot,\cdot)$ is the Marcum Q-function, defined by

$Q_{1}(\mu,\nu)=\int_{\nu}^{\infty}x\,e^{-(x^{2}+\mu^{2})/2}\,I_{0}(\mu
x)\,dx\;,$ (40)

and $I_{0}(\cdot)$ is the zeroth-order modified Bessel function of the first
kind.

The probability distribution for $|a|^{2}$ can be related to the familiar
Poisson distribution of photon number in coherent states as follows. For
integers $k\geq 1$, the moments $\mathsf{E}[|a|^{2k}]$ all exist and,
therefore, uniquely determine the probability distribution of $|a|^{2}$. We
observe that the number operator $\hat{n}=\hat{a}^{\dagger}\hat{a}$ provides a
quantum mechanical analog of $|a|^{2}$ in the sense that
$\bra{\alpha}\mathcal{S}(\hat{n}^{k})\ket{\alpha}=\mathsf{E}[|a|^{2k}]$, where
$\mathcal{S}(\hat{n}^{k})$ is the symmetrized form (or Weyl ordering) of
$\hat{n}^{k}$, by the optical equivalence theorem [35]. Furthermore,
$\mathcal{S}(\hat{n}^{k})$ can be written as a degree-$k$ polynomial
$\mathcal{P}_{k}(\hat{n})$ in $\hat{n}$; for example,
$\mathcal{P}_{1}(\hat{n})=\hat{n}+\frac{1}{2}$,
$\mathcal{P}_{2}(\hat{n})=\hat{n}^{2}+\hat{n}+\frac{1}{2}$, etc. [36]. In
terms of the photon number basis, then, we may write the moments as

$\mathsf{E}[|a|^{2k}]=\bra{\alpha}\mathcal{S}(\hat{n}^{k})\ket{\alpha}=\sum_{n=0}^{\infty}\mathcal{P}_{k}(n)\,|\braket{n}{\alpha}|^{2}\;,$
(41)

where $|\braket{n}{\alpha}|^{2}=e^{-|\alpha|^{2}}|\alpha|^{2n}/n!$ is the
probability associated with photon number $n$. The discrete Poisson
distribution for photon number could therefore be interpreted as a
calculational device, a mathematical artifice, so to speak, for determining
the distribution of the continuous random variable $|a|^{2}$. The begs the
question of whether threshold exceedances may be interpreted as photon
detection events.

The probability of such a detection event is given by

$\begin{split}\Pr\left[|a|^{2}>\gamma^{2}\right]&=Q_{1}\left(2|\alpha|,2\gamma\right)=\int_{2\gamma}^{\infty}xe^{-(x^{2}+4|\alpha|^{2})/2}\,\frac{1}{\pi}\int_{0}^{\pi}e^{2|\alpha|x\cos\theta}d\theta\,dx\;.\end{split}$
(42)

We now note that, to fourth order in $|\alpha|$,

$\begin{split}\Pr\left[|a|^{2}>\gamma^{2}\right]&=\int_{2\gamma}^{\infty}xe^{-x^{2}/2}\left[1+(x^{2}-2)|\alpha|^{2}+\tfrac{1}{4}(x^{4}-8x^{2}+8)|\alpha|^{4}\right]dx+\mathcal{O}(|\alpha|^{6})\\\
&=e^{-2\gamma^{2}}\left(1+4\gamma^{2}|\alpha|^{2}+4\gamma^{2}(\gamma^{2}-1)|\alpha|^{4}\right)+\mathcal{O}(|\alpha|^{6})\;.\end{split}$
(43)

The presence of $|\alpha|^{2}$ is the lowest-order approximation is the first
indication of the emergence of the Born rule, although the correspondence is
subtle and requires some discussion.

According to quantum mechanics, the probability of observing $n$ photons given
a coherent state $\ket{\alpha}$ is
$p_{n}=|\braket{n}{\alpha}|^{2}=e^{-|\alpha|^{2}}|\alpha|^{2n}/n!$. Hence, the
probability of observing no photons at all is $p_{0}=e^{-|\alpha|^{2}}$, while
the probability of observing at least one photon is
$1-p_{0}=1-e^{-|\alpha|^{2}}\approx|\alpha|^{2}$ for $|\alpha|\ll 1$.
According to equation (42), for $\alpha=0$ we have
$\Pr[|a|^{2}>\gamma^{2}]=e^{-2\gamma^{2}}$, which we interpret as the dark
count probability of the vacuum state. For $\alpha\neq 0$, equation (43) will
be a good approximation for $|\alpha|^{2}\ll 1/(4\gamma^{2})$. Furthermore,
for $\gamma^{2}\gg\frac{1}{2}$ we will have low dark counts. So, for
$|\alpha|^{2}\ll 1/(4\gamma^{2})\ll\frac{1}{2}$ (i.e., $|\alpha|$ small and
$\gamma$ large) we expect to be in the near-single-photon regime. However,
taking both $|\alpha|$ to be small and $\gamma$ to be large does not
necessarily provide the best agreement with the Born rule, as we shall see.

Figure 1 shows an example using $\alpha=0.707\cos\theta$, $\gamma=1$, and
$N=10^{4}$ random realizations. Examining $N\Pr[|a|^{2}>\gamma^{2}]$ as a
function of $\theta$, we observe a near-perfect sinusoidal pattern with a
period of $\pi$ that has a minimum of $Ne^{-2\gamma^{2}}\approx 1353$ and a
maximum of $NQ_{1}(1.414,2)\approx 3942$. Subtracting the dark counts and
renormalizing by the resulting maximum value, as one normally does in
practice, gives a good approximation to the $\cos^{2}\theta$ probability law
one would expect for an application of the Born rule to single-photon
detection. Furthermore, reducing the magnitude of $\alpha$, and of course
ignoring the many non-detection events, gives arbitrarily good agreement. (If
$\alpha$ is identically zero we will have a constant dark count rate which,
when subtracted out, gives the quantum mechanical prediction of zero.)

Figure 1: Plot of simulated counts (blue circles) for $N=10^{4}$ trials versus
the linear polarization angle $\theta$ for $\alpha=0.707\cos\theta$ and
$\gamma=1$. The detection probability (thick solid black line) given by
equation (42) and its approximation (thin dashed red line) given by equation
(43), both scaled by $N$, are shown as well.

Now, for a general coherent state $\ket{\alpha}$, quantum mechanics does not
actually predict a probability of $\cos^{2}\theta$, as our detector only
indicates the presence of one _or more_ photons. The actual predicted
probability is $1-e^{-|\alpha|^{2}}$, which is only approximately sinusoidal.
Comparing this to equation (42), suitably normalized, we observe a subtle
difference. For $\alpha=\cos\theta$ and $\gamma=1$, our model predicts a
slightly lower probability than the quantum prediction of
$1-e^{-|\alpha|^{2}}$. (See figure 2.) For $\gamma=0.5$ it is slightly higher.
Treating $\gamma$ as an adjustable parameter, then, allows for an arbitrarily
good fit.

Figure 2: Plot of predicted probabilities for a coherent state with
$\alpha=\cos\theta$ and $\gamma=1$. The black solid line is the baseline
$\cos^{2}\theta$ prediction. The dashed red line is our prediction based on a
normalized version of equation (42). The dotted blue line is the quantum
mechanical prediction for detecting one or more photons.

The Poissonian nature of detection events is also an important characteristic
of coherent light. Strictly speaking, the single-mode approximation we have
made entails an infinite coherence time, so the temporal distribution of
events is not well defined within the present model. Under realistic
conditions, however, there will be some non-zero bandwidth associated with any
given mode and, hence, an associated coherence time $\tau$ inversely
proportional to this bandwidth. For laser light, this is typically on the
order of nanoseconds to microseconds. We may therefore suppose that over a
time $T\gg\tau$, there are about $N=T/\tau$ independent detection
opportunities. Under our model, the probability of a detection for each such
opportunity, after subtracting for dark counts and assuming $\gamma|\alpha|\ll
1$, is about $p=4\gamma^{2}e^{-2\gamma^{2}}|\alpha|^{2}$. Now, the number of
detections out of these $N$ trials is binomially distributed with a mean of
$Np$ and a variance of $Np(1-p)$. If $N$ is large and $p$ is small (e.g.,
$|\alpha|\ll 1$ for fixed $\gamma$), this distribution becomes approximately
Poissonian. The average number of counts, $Np$, may be interpreted as the
product of an incident photon rate $|\alpha|^{2}/\tau$, a detection efficiency
$4\gamma^{2}e^{-2\gamma^{2}}$, and an observation time $T$. Deviations from
Poissonian behavior are to be expected when $|\alpha|$ is large, the detection
efficiency is high, or the observation time is short.

A further comparison to experimental observations can be made by considering
detection efficiency. For Poisson-distributed photon statistics,
experimentalists often use a parametric model of the form

$p=1-(1-\delta)e^{-\eta|\alpha|^{2}}\;,$ (44)

where $p$ is the probability of a detection event, $\delta$ is the dark count
probability, and $\eta\in[0,1]$ is the detection efficiency [37, 38]. Our
model conforms with this general expression in the small $|\alpha|$ limit if
we take $\delta=e^{-2\gamma^{2}}$ and

$\eta=\frac{4\gamma^{2}\,e^{-2\gamma^{2}}}{1-e^{-2\gamma^{2}}}\;.$ (45)

Note that, in this interpretation, the effective detection efficiency
increases as the threshold $\gamma$ is decreased, attaining near unit
efficiency for $\gamma\approx 0.8$; however, for values much lower than this
the efficiency is over unity and this interpretation is no longer valid.

Finally, another important quantity in experimental quantum optics is the
interferometric visibility, which measures the degree of coherence in the
prepared state. This may be defined as the ratio of the difference in maximum
and minimum probabilities to their sum, which in our case is

$\mathcal{V}=\frac{Q_{1}(2|\alpha|,2\gamma)-e^{-2\gamma^{2}}}{Q_{1}(2|\alpha|,2\gamma)+e^{-2\gamma^{2}}}\;.$
(46)

Taking $\gamma$ to be large, with $\alpha$ fixed, therefore gives a fringe
visibility arbitrarily close to unity. For larger values of $|\alpha|$,
corresponding more closely to the classical regime, the convergence to unity
occurs more rapidly. We illustrate this in figure 3, plotting visibility as a
function of the threshold for different values of $\alpha$. It is important to
note that the visibility described here is in terms of the probability of
detection, $\Pr[|a|^{2}>\gamma^{2}]$, not the intensity, $|a|^{2}$, which is
random, nor the expected intensity, $\mathsf{E}[\,|a|^{2}]$, which would give
a visibility of one half. This point is important for a proper comparison with
quantum mechanics, which predicts visibilities as high as one for actual
measured counts, not classical intensities.

Figure 3: Plot of the interferometric visibility as a function of the
detection threshold $\gamma$. The lower blue dotted curve is for $\alpha=0.5$,
the middle black solid curve is for $\alpha=1$, and the upper red dashed curve
is for $\alpha=1.5$.

In summary, we expect to find good agreement with quantum mechanical
predictions when $|\alpha|$ is small and $\gamma$ is large, with optimal but
imperfect agreement for suitable finite, nonzero choices of both. Parameter
regimes in which agreement with the Born rule is best, such as low values of
$|\alpha|$, may produce lower visibility. In general, no combination of
$\alpha$ and $\gamma$ can give both arbitrarily good agreement with the Born
rule and a visibility arbitrarily close to unity. The detailed tradeoffs
between these and other competing metrics are discussed further in section 4.

###  3.3 Dual-mode Detection

Previously, we considered measurements along a single polarization mode and
found that the associated probabilities follow the Born rule, albeit with a
threshold-dependent rescaling and fixed offset in accordance with equation
(44). Such measurements cannot distinguish between a missed detection and an
event that would have resulted in a detection in an orthogonal polarization.
Dual-modal detection provides an alternative method for comparing against the
Born rule that overcomes this deficiency.

Let $\hat{\boldsymbol{e}}_{H}$ and $\hat{\boldsymbol{e}}_{V}$ denote the
horizontal and vertical polarization unit vectors for given wave vector. A
linearly polarized coherent state for this wave vector may be described by the
complex amplitude vector

$\boldsymbol{a}=a_{H}\hat{\boldsymbol{e}}_{H}+a_{V}\hat{\boldsymbol{e}}_{V}=\left(\alpha\cos\theta+\frac{z_{H}}{\sqrt{2}}\right)\hat{\boldsymbol{e}}_{H}+\left(\alpha\sin\theta+\frac{z_{V}}{\sqrt{2}}\right)\hat{\boldsymbol{e}}_{V}\;,$
(47)

where $z_{H}$ and $z_{V}$ are independent standard complex Gaussian random
variables. Note that $\theta=0$ corresponds to a vacuum state in the
orthogonal polarization mode, which is always assumed to be present.

As a consequence of the independence of $z_{H}$ and $z_{V}$, the random
variables $|a_{H}|^{2}$ and $|a_{V}|^{2}$ are also independent, and their
joint cdf is given by the product of their marginal distributions. Let us
suppose a dual-mode detector that will register separate events if either
$|a_{H}|>\gamma$ or $|a_{V}|>\gamma$ are true. This would be the case if the
detector were, say, a pair of bifurcating harmonic oscillators oriented in the
horizontal and vertical polarization directions. Equivalently, we may consider
a polarizing beam splitter that separates the components to two single-mode
detectors. The probability of no detection occurring is then

$\begin{split}P_{0}&=\Pr\\!\left[|a_{H}|^{2}\leq\gamma^{2},|a_{V}|^{2}\leq\gamma^{2}\right]\\\
&=\Bigl{[}1-Q_{1}(2|\alpha\cos\theta|,2\gamma)\Bigr{]}\Bigl{[}1-Q_{1}(2|\alpha\sin\theta|,2\gamma)\Bigr{]}\end{split}$
(48)

Likewise, the probabilities for the three possible detection events are

$\displaystyle P_{H}$
$\displaystyle=Q_{1}(2|\alpha\cos\theta|,2\gamma)\Bigl{[}1-Q_{1}(2|\alpha\sin\theta|,2\gamma)\Bigr{]}$
(49) $\displaystyle P_{V}$
$\displaystyle=\Bigl{[}1-Q_{1}(2|\alpha\cos\theta|,2\gamma)\Bigr{]}Q_{1}(2|\alpha\sin\theta|,2\gamma)$
(50) $\displaystyle P_{HV}$
$\displaystyle=Q_{1}(2|\alpha\cos\theta|,2\gamma)\,Q_{1}(2|\alpha\sin\theta|,2\gamma)\;,$
(51)

where $P_{H}$ is the probability of a single detection of $H$, $P_{V}$ is the
probability of a single detection of $V$, and $P_{HV}$ is the probability of
both.

In actual experiments with coherent light it is common to reject events in
which there are two detections and, of course, ignore those with none. Out of
a notional, unknown number $N$ of independent trials, one measures
$S_{H}=NP_{H}$ single counts for $H$, $S_{V}=NP_{V}$ single counts for $V$,
and $NP_{HV}$ “accidental” coincidence counts. If we post-select on single
detection events, the conditional probability, $p_{H}$, of detecting $H$ is

$p_{H}=\frac{S_{H}}{S_{H}+S_{V}}=\frac{P_{H}}{P_{H}+P_{V}}\;.$ (52)

We may now compare $p_{H}$ with the Born rule prediction of $\cos^{2}\theta$.
An example is plotted in figure 4 for $|\alpha|^{2}=0.5$ and $\gamma=1$. The
agreement is perfect when $\theta=45^{\circ},135^{\circ}$ (diagonal and anti-
diagonal polarization, respectively), resulting in balanced probabilities and
a conditional probability of $\frac{1}{2}$. For other values of $\theta$, we
find $\frac{1}{2}(1-\mathcal{V})\leq p_{H}\leq\frac{1}{2}(1+\mathcal{V})$,
where

$\mathcal{V}=\frac{Q_{1}(2|\alpha|,2\gamma)(1-e^{-2\gamma^{2}})-[1-Q_{1}(2|\alpha|,2\gamma)]e^{-2\gamma^{2}}}{Q_{1}(2|\alpha|,2\gamma)(1-e^{-2\gamma^{2}})+[1-Q_{1}(2|\alpha|,2\gamma)]e^{-2\gamma^{2}}}$
(53)

is the visibility. For our particular case, $\mathcal{V}=0.61$, so $0.19\leq
p_{H}\leq 0.81$.

The maximum discrepancy arises when $\theta$ is $0^{\circ}$ or $90^{\circ}$.
For $\theta=0^{\circ}$, the polarization of the wave is horizontal, but we are
still not guaranteed an $H$ outcome, even conditionally, because the
probability of a “false” $V$ detection is still nonzero. Similarly, for
$\theta=90^{\circ}$, the polarization of the wave is vertical, but an $H$
outcome is still possible due to dark counts. In any realistic experiment,
such events will be unavoidable and are quantified by a visibility below
unity. Such anomalous events are effectively removed by renormalization,
resulting in a modified conditional probability of the form

$\hat{p}_{H}=\frac{1}{\mathcal{V}}\left(p_{H}-\frac{1}{2}\right)+\frac{1}{2}\;.$
(54)

This renormalized conditional probability gives excellent agreement with the
Born rule prediction, as shown in figure 4.

Figure 4: Plot of the conditional probability $p_{H}$ (blue dotted line) and
renormalized conditional probability $\hat{p}_{H}$ (red dashed line) against
the $\cos^{2}\theta$ Born rule prediction (black solid line), versus the
polarization angle $\theta$ for $|\alpha|^{2}=0.5$ and $\gamma=1$.

###  3.4 Particle-like Behavior

Consider a coherent state prepared in some polarization mode
$\hat{\boldsymbol{e}}_{0}$ and spatial mode $\boldsymbol{k}_{R}$ traveling to
the right that is incident upon a 50/50 beam splitter (BS). The outgoing beams
have orthogonal spatial modes of $\boldsymbol{k}_{R}$ and $\boldsymbol{k}_{D}$
traveling right and down, respectively, each with the same polarization mode.
The initial state may be described by the vector

$\boldsymbol{a}=\begin{bmatrix}\alpha+z_{R}/\sqrt{2}\\\
z_{D}/\sqrt{2}\end{bmatrix}=\alpha\begin{bmatrix}1\\\
0\end{bmatrix}+\frac{1}{\sqrt{2}}\begin{bmatrix}z_{R}\\\
z_{D}\end{bmatrix}\;,$ (55)

where $z_{R}$ and $z_{D}$ are independent standard complex Gaussian random
variables corresponding to the ZPF components of the two spatial modes. For
simplicity, we ignore the orthogonal polarization modes.

The beam splitter acts as a Hadamard gate $\boldsymbol{\mathsf{H}}$,
transforming $\boldsymbol{a}$ into

$\boldsymbol{a}^{\prime}=\boldsymbol{\mathsf{H}}\boldsymbol{a}=\frac{\alpha}{\sqrt{2}}\begin{bmatrix}1\\\
1\end{bmatrix}+\frac{1}{2}\begin{bmatrix}z_{R}+z_{D}\\\
z_{R}-z_{D}\end{bmatrix}\;.$ (56)

Note that $z_{R}^{\prime}=(z_{R}+z_{D})/\sqrt{2}$ and
$z_{D}^{\prime}=(z_{R}-z_{D})/\sqrt{2}$ are again independent standard complex
Gaussian random variables, so the noise term for $\boldsymbol{a}^{\prime}$ has
the same form as that for $\boldsymbol{a}$.

If we place single-mode detectors at each output port of the beam splitter,
there will be four possible outcomes with four corresponding probabilities: no
detections ($P_{0}$), a single detection for mode $\boldsymbol{k}_{R}$
($P_{R}$), a single detection for mode $\boldsymbol{k}_{D}$ ($P_{D}$), and
coincident detections on both modes ($P_{RD}$). These probabilities are as
follows:

$\displaystyle P_{0}$
$\displaystyle=\Bigl{[}1-Q_{1}(\sqrt{2}|\alpha|,2\gamma)\Bigr{]}^{2}$ (57)
$\displaystyle P_{R}$
$\displaystyle=P_{D}=\Bigl{[}1-Q_{1}(\sqrt{2}|\alpha|,2\gamma)\Bigr{]}Q_{1}(\sqrt{2}|\alpha|,2\gamma)$
(58) $\displaystyle P_{RD}$
$\displaystyle=Q_{1}(\sqrt{2}|\alpha|,2\gamma)^{2}\;.$ (59)

Note that $P_{RD}\geq P_{R}P_{D}$, since each detection event is independent
of the other. A similar result is found in the semiclassical treatment of
photon detection [39]. In the single-photon regime ($|\alpha|\ll 1$) one would
expect particle-like behavior, so coincident detections should be quite rare.
Quantum mechanically, the probability of a coincident detection for a true,
single-photon state would be exactly zero.

Experimentally, one counts the number of single-detection events,
$S_{R}=NP_{R}$ and $S_{D}=NP_{D}$, for transmitted and reflected light,
respectively, as well as the number of coincidences, $C=NP_{RD}$, where $N$ is
the nominal number of trials. The difficulty with such experiments is that $N$
is often unknown or perhaps unknowable. If $N$ is known, the ratio
$R=CN/(S_{R}S_{D})$, more commonly associated with the degree of second-order
temporal coherence $g^{(2)}(0)$, would be expected to have a value no less
than one, since

$R=\frac{CN}{S_{R}S_{D}}=\frac{P_{RD}}{P_{R}P_{D}}\geq 1\;.$ (60)

If $C=0$, as predicted by quantum mechanics, and $S_{R},S_{D}>0$, then $R=0$,
thereby violating the inequality. Early experiments of this sort were
performed by Grangier et al. using both a light-emitting diode (LED) [40] and
a heralded photon source [41]. The LED light source was turned on briefly
using a controlled electronic trigger, allowing $N$ to be know precisely.
Since the LED light was strongly attenuated, a value of $R$ near unity, and
consistent with the inequality $R\geq 1$, was measured, as one might expect.

In the case of the heralded photon source, $N$ was taken to be the number of
trigger events, $N_{t}$, each of which was taken to indicate the presence of a
single, heralded photon. Under this assumption, the experimenters obtained a
value of $R_{t}=CN_{t}/(S_{R}S_{D})$ significantly less than one. A value less
than one is generally considered to be evidence of photon antibunching. The
true value of $N$, however, could not be known and may well have been much
larger than $N_{t}$, in which case a value below unity would not be
surprising. A similar experiment, also using heralded events, was performed
recently by Thorn et al. using a modern parametric downconversion source and
avalanche photodiodes, with similar results [42].

For our model, the single-photon regime provides a good approximation to a
true, single-photon state, so long as we ignore non-detection events. Taking
$N$ to be $N_{d}=S_{R}+S_{D}+C$, the total number of detection events, we
obtain a result similar to heralding. From this we may compute the ratio

$R_{d}=\frac{CN_{d}}{S_{R}S_{D}}=\frac{P_{RD}(1-P_{0})}{P_{R}P_{D}}\;.$ (61)

This may equivalently be seen as replacing the absolute probabilities
$P_{R},P_{D},P_{RD}$ in the expression for $R$ with the _conditional_
probabilities
$p_{R}=P_{R}/(1-P_{0}),p_{D}=P_{D}/(1-P_{0}),p_{RD}=P_{RD}/(1-P_{0})$. Such
conditioning is similar to the experimental procedure of using heralding to
define the number of trials. It is now easy to show that $R_{d}$ can be less
than unity when either $|\alpha|\ll 1$ or $\gamma\gg 1$. As an example, figure
5 shows the values of $R$ and $R_{d}$ as a function of $|\alpha|$ for
$\gamma=1$. For this example, $R_{d}$ achieves a minimum value of about 0.34,
whereas quantum mechanics predicts $R=0$ for an ideal single-photon state. By
increasing the value of $\gamma$, however, this minimum can be made
arbitrarily small. For example, the observed value of $R_{d}=0.018$ in
reference [42] could be achieved with $\alpha=0.3$ and $\gamma=1.6$. In this
way, a purely classical model of light, when analyzed in a similar way, can
exhibit the same anomalous quantum behavior.

Figure 5: Plot of $R$ (dotted blue line) and $R_{d}$ (solid red line) as a
function of $|\alpha|$ for $\gamma=1$. Values below one indicate particle-like
behavior.

##  4 General Linear Transformations

###  4.1 Single-Photon, Four-Mode Entanglement

In quantum mechanics, a single photon can be entangled across multiple modes.
Similar behavior can be modeled classically. Consider a coherent state
prepared with polarization $\hat{\boldsymbol{e}}_{H}$ traveling to the right
and incident upon a 50/50 beam splitter (BS). The initial, four-mode state may
be written

$\boldsymbol{a}=\begin{bmatrix}\left(\alpha+\frac{z_{1}}{\sqrt{2}}\right)\hat{\boldsymbol{e}}_{H}+\frac{z_{2}}{\sqrt{2}}\hat{\boldsymbol{e}}_{V}\\\
\frac{z_{3}}{\sqrt{2}}\hat{\boldsymbol{e}}_{H}+\frac{z_{4}}{\sqrt{2}}\hat{\boldsymbol{e}}_{V}\end{bmatrix}=\alpha\begin{bmatrix}1\\\
0\\\ 0\\\ 0\end{bmatrix}+\frac{1}{\sqrt{2}}\begin{bmatrix}z_{RH}\\\ z_{RV}\\\
z_{DH}\\\ z_{DV}\end{bmatrix}$ (62)

where $z_{RH},z_{RV},z_{DH},z_{DV}$ are independent and identically
distributed (iid) standard complex Gaussian random variables arising from the
zero-point field and corresponding to the two spatial modes
($\boldsymbol{k}_{R}$ and $\boldsymbol{k}_{D}$) and polarization modes
($\hat{\boldsymbol{e}}_{H}$ and $\hat{\boldsymbol{e}}_{V}$).

After the beam splitter, the state becomes

$\boldsymbol{a}^{\prime}=(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})\boldsymbol{a}=\frac{\alpha}{\sqrt{2}}\begin{bmatrix}1\\\
0\\\ 1\\\ 0\end{bmatrix}+\frac{1}{2}\begin{bmatrix}z_{RH}+z_{DH}\\\
z_{RV}+z_{DV}\\\ z_{RH}-z_{DH}\\\ z_{RV}-z_{DV}\end{bmatrix}\;,$ (63)

where $\otimes$ is the Kronecker product and $\boldsymbol{\mathsf{I}}$ is the
$2\times 2$ identity. Finally, we may apply an $\boldsymbol{\mathsf{X}}$ gate
(i.e., a half-wave plate rotated $45^{\circ}$) on the downward mode to change
the polarization. This has the effect of performing a controlled NOT gate
$\boldsymbol{\mathsf{C}}$, with the spatial mode (i.e., the optical path) as
the control and the polarization as the target. The resulting state is now

$\boldsymbol{a}^{\prime\prime}=\boldsymbol{\mathsf{C}}(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})\boldsymbol{a}=\frac{\alpha}{\sqrt{2}}\begin{bmatrix}1\\\
0\\\ 0\\\ 1\end{bmatrix}+\frac{1}{2}\begin{bmatrix}z_{RH}+z_{DH}\\\
z_{RV}+z_{DV}\\\ z_{RV}-z_{DV}\\\ z_{RH}-z_{DH}\end{bmatrix}$ (64)

Note that the second term is again a vector of iid standard complex Gaussian
random variables. We may therefore rewrite $\boldsymbol{a}^{\prime\prime}$ as

$\boldsymbol{a}^{\prime\prime}=\frac{\alpha}{\sqrt{2}}\begin{bmatrix}1\\\ 0\\\
0\\\ 1\end{bmatrix}+\frac{1}{\sqrt{2}}\begin{bmatrix}z_{RH}^{\prime}\\\
z_{RV}^{\prime}\\\ z_{DH}^{\prime}\\\
z_{DV}^{\prime}\end{bmatrix}=\alpha\,\boldsymbol{\psi}+\frac{\boldsymbol{z}^{\prime}}{\sqrt{2}}\;,$
(65)

where
$\boldsymbol{z}^{\prime}=[z_{RH}^{\prime},\ldots,z_{DV}^{\prime}]^{\mathsf{T}}$
and $\boldsymbol{\psi}$ is a column vector of unit amplitude. The vector
$\boldsymbol{\psi}$ has the mathematical form of an entangled Bell state

$\boldsymbol{\psi}=\frac{\ket{R,H}+\ket{D,V}}{\sqrt{2}}\;,$ (66)

where
$\ket{R,H}=\ket{R}\otimes\ket{H}=[1,0]^{\mathsf{T}}\otimes[1,0]^{\mathsf{T}}$
and
$\ket{D,V}=\ket{D}\otimes\ket{V}=[0,1]^{\mathsf{T}}\otimes[0,1]^{\mathsf{T}}$.

To perform a measurement of all four modes, each spatial mode is put into a
dual-mode detector and threshold detection is performed. There are four
components and, so, 16 possible outcomes, including multiple detections. For
$|\alpha|\ll 1$, the most likely outcome is no detections at all, with single
detections being the next most likely outcome. At the opposite extreme, for
$|\alpha|\gg 1$ the most likely outcome is detection on all four modes. For
small values of $|\alpha|$, the probability of a single detection on either
$\ket{R,H}$ or $\ket{D,V}$ (both equally likely) is much more likely than a
single detection on $\ket{D,H}$ or $\ket{R,V}$.

Let $\Pr[R,H]=\Pr[D,V]$ and $\Pr[R,V]=\Pr[D,H]$ denote the probabilities for
single-detection events on each of the four modes. These will be given by

$\displaystyle\Pr[R,H]$
$\displaystyle=\Pr[D,V]=\left(1-e^{-2\gamma^{2}}\right)^{2}Q_{1}(\sqrt{2}|\alpha|,2\gamma)\Bigl{[}1-Q_{1}(\sqrt{2}|\alpha|,2\gamma)\Bigr{]}$
(67) $\displaystyle\Pr[R,V]$
$\displaystyle=\Pr[D,H]=e^{-2\gamma^{2}}\left(1-e^{-2\gamma^{2}}\right)\Bigl{[}1-Q_{1}(\sqrt{2}|\alpha|,2\gamma)\Bigr{]}^{2}$
(68)

These probabilities are illustrated in figure 6 for $\alpha=1$. We see that
the dominant modes peak in probability at threshold values somewhat greater
than $1$ but are relatively much larger than the other two modes. This
comports with the general behavior one would expect of a single-photon state
that is hyperentangled in spatial and polarization modes [43]. If we consider
_only_ single-mode detection events (i.e., detections on one spatial mode and
one polarization mode), then the conditional probability of each dominant mode
converges to 0.5, the ideal quantum prediction, when $\gamma$ is large.
Conversely, the conditional probability converges to a nonzero value, which is
dependent on $\alpha$, when $\gamma$ is small. Qualitatively similar behavior
is found when $|\alpha|$ is varied while holding $\gamma$ fixed. Thus, a
correspondence with quantum mechanical predictions is achieved, but only in
the limit of larger threshold values and only when one post-selects on single-
mode detection events. This asymptotic behavior is a result of the symmetry of
the Bell state and would not be expected more generally. Correlations between
the modes are purely a result of post-selection, as the modes themselves are
statistically independent. Of course, this model reproduces only _local_
correlations between single-photon modes and not the _nonlocal_ correlations
one expects from a multi-photon entangled state.

Figure 6: Plots of $\Pr[R,H]=\Pr[D,V]$ (red dashed line), $\Pr[R,V]=\Pr[D,H]$
(dotted blue line), and the conditional probability of mode $\ket{R,H}$ given
any single-mode detection (solid black line) versus $\gamma$ and for
$|\alpha|=1$.

###  4.2 Wave/Particle Duality

In quantum mechanics, photons can exhibit both particle- and wave-like
behavior. This, too, can be modeled classically. Consider, an initial quantum
state $\ket{R,H}$ that undergoes a transformation via a 50/50 beam splitter
and a phase shifter in the $\ket{D,H}$ mode. Using a pair of mirrors, the two
paths are recombined in a second beam splitter to form a Mach-Zehnder
interferometer. The two output ports are then measured with detectors. Quantum
mechanically, the final state (before measurement) is

$\begin{split}\boldsymbol{\psi}&=(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})(\boldsymbol{\mathsf{R}}_{\phi}\otimes\boldsymbol{\mathsf{I}})(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})\ket{R,H}\\\
&=\tfrac{1}{2}(1+e^{i\phi})\ket{R,H}\;+\;\tfrac{1}{2}(1-e^{i\phi})\ket{D,H}\;,\end{split}$
(69)

where $\boldsymbol{\mathsf{R}}_{\phi}$ is the phase shift gate

$\boldsymbol{\mathsf{R}}_{\phi}=\begin{pmatrix}1&0\\\
0&e^{i\phi}\end{pmatrix}\;.$ (70)

Accordingly, the probability of finding a photon in the $\ket{R,H}$ mode is
$\cos^{2}(\phi/2)$.

We can model the problem classically by starting with an initial coherent
state $\boldsymbol{a}=\alpha\ket{R,H}+\boldsymbol{z}/\sqrt{2}$ and
transforming it via the same linear operations into

$\boldsymbol{a}_{\rm
MZ}=(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})(\boldsymbol{\mathsf{R}}_{\phi}\otimes\boldsymbol{\mathsf{I}})(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})\boldsymbol{a}=\alpha\boldsymbol{\psi}+\frac{\boldsymbol{z}^{\prime}}{\sqrt{2}}\;.$
(71)

The conditional probability of a detection in mode $\ket{R,H}$, given a single
detection in either mode $\ket{R,H}$ or $\ket{D,H}$, is now found to be

$p_{\rm
MZ}(\phi)=\left(1+\frac{Q_{1}(|\alpha(1-e^{i\phi})|,2\gamma)}{Q_{1}(|\alpha(1+e^{i\phi})|,2\gamma)}\cdot\frac{1-Q_{1}(|\alpha(1+e^{i\phi})|,2\gamma)}{1-Q_{1}(|\alpha(1-e^{i\phi})|,2\gamma)}\right)^{-1}\;.$
(72)

The resulting interference pattern, as shown in figure 7, is similar to what
one would expect from classical light if one were observing intensities;
however, we are showing probabilities. The pattern also reflects the nonlocal
(i.e., spatially distributed) nature of the interferometer: light travels
along both arms and interferes only when recombined. In this way, our
classical model exhibits the wave-like nature of light in terms of discrete
detection events.

Figure 7 also affords an interesting comparison to experimental data.
Reference [44] describes an experimental realization of Wheeler’s delayed-
choice experiment using a nitrogen-vacancy center single-photon source. The
experimenters report an anti-correlation parameter of $R_{d}=0.12$ and a
visibility of $\mathcal{V}=94\%$, after subtracting dark counts, when the
second beam splitter is in place. An estimated 2600 photons were used for each
measurement sample, and the resulting data points were fit to a cosine wave
with a period of $2\pi$. We have followed a similar analysis procedure in
figure 7, using a set of fitted samples with a standard deviations of
$1/\sqrt{2600}$ and subtracting the expected dark count probability of
$e^{-2\gamma^{2}}$. Using parameter values of $\gamma=1.6$ and $\alpha=0.95$
gives $R_{d}=0.12$, $\mathcal{V}=94\%$, and a root-mean-square error of 0.04,
in excellent agreement with the experimental results. Although the two light
sources are quite different, our model is nevertheless capable of reproducing
similar observations.

Figure 7: Plot of the idealized Mach-Zehnder interference pattern given by
$\cos^{2}(\phi/2)$ (black solid line) and the interference pattern predicted
by $p_{\rm MZ}(\phi)$ of equation (72) (blue dashed line) for $\alpha=0.95$
and $\gamma=1.6$. The red circles are fitted samples analogous to those in
reference [44].

To recover the particle-like nature of the system, we may create a Wheeler
delayed-choice experiment by removing the final beam splitter The resulting
state is now

$\boldsymbol{a}_{\rm
DC}=(\boldsymbol{\mathsf{R}}_{\phi}\otimes\boldsymbol{\mathsf{I}})(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})\boldsymbol{a}=\frac{\alpha}{\sqrt{2}}\Bigl{(}\ket{R,H}+e^{i\phi}\ket{D,H}\Bigr{)}+\frac{\boldsymbol{z}^{\prime\prime}}{\sqrt{2}}\;.$
(73)

The conditional probability of a detection in mode $\ket{R,H}$ given a single
photon detection in either mode $\ket{R,H}$ or $\ket{D,H}$ is now
$\frac{1}{2}$, independent of $\phi$. In addition, the probability of a double
detection in both modes may be made arbitrarily small by decreasing $|\alpha|$
or, equivalently, increasing $\gamma$. This is the behavior one would expect
from a localized particle. Note that it does not matter when the choice to
remove the final beam splitter is made.

A similar result is obtained if we simply provide “which way” information by
marking one of the two arms with, say, a change in polarization. If we apply
an $\boldsymbol{\mathsf{X}}$ gate on the lower arm before the final beam
splitter, the resulting state will be

$\begin{split}\boldsymbol{a}_{WW}&=(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})\boldsymbol{\mathsf{C}}(\boldsymbol{\mathsf{R}}_{\phi}\otimes\boldsymbol{\mathsf{I}})(\boldsymbol{\mathsf{H}}\otimes\boldsymbol{\mathsf{I}})\boldsymbol{a}\\\
&=\frac{\alpha}{2}\Bigl{(}\ket{R,H}+e^{i\phi}\ket{R,V}+\ket{D,H}-e^{i\phi}\ket{D,V}\Bigr{)}+\frac{\boldsymbol{z}^{\prime\prime\prime}}{\sqrt{2}}\;.\end{split}$
(74)

The interference pattern is once again lost (i.e., the detection probabilities
are independent of $\phi$). Each of the four modal outcomes occurs now with
equal probability, with the likelihood of multiple detections again vanishing
as $|\alpha|$ is decreased or $\gamma$ is increased. Replacing the NOT gate
with a unitary gate that is close to the identity will result in a diminished
but still discernible interference pattern, so one may consider measuring the
path information weakly as well. So, if there is only partial which-way
information, the interference pattern will simply diminish by degrees.

Finally, if one considers the total number of single-detection events on
either mode, it is natural to suppose that this should be insensitive to
whether the final beam splitter is present or not. Under the present model,
this need not be the case. Although the total intensity is the same, since the
beam splitter constitutes a unitary transformation, the probability of a count
on either detector _with_ the beam splitter is

$\begin{split}P_{\rm
MZ}&=1-\left[1-Q_{1}\left(|\alpha(1+e^{i\phi})|,2\gamma\right)\right]\left[1-Q_{1}\left(|\alpha(1-e^{i\phi})|,2\gamma\right)\right]\;,\end{split}$
(75)

while the probability of this event _without_ the final beam splitter is

$\begin{split}P_{\rm
DC}&=1-\left[1-Q_{1}\left(\sqrt{2}|\alpha|,2\gamma\right)\right]\left[1-Q_{1}\left(\sqrt{2}|\alpha
e^{i\phi}|,2\gamma\right)\right]\;.\end{split}$ (76)

In general, these two probabilities are different; quantum mechanics predicts
that they should be the same. Nevertheless, for fixed $\gamma$ and
$|\alpha|\to 0$ we do find that $P_{\rm MZ}/P_{\rm DC}\to 1$, as one might
expect. This is consistent with previous observation that the approximation
$Q_{1}(2|\alpha|,2\gamma)\approx e^{-2\gamma}(1+4\gamma^{2}|\alpha|^{2})$ is
valid when $\gamma|\alpha|$ is small but not necessarily otherwise. The
subtlety of this relationship will elaborated upon further in the next
section.

###  4.3 General Multimodal States

The transformation of coherent light via a sequence of linear optical
components can be described, in general, by a $d\times d$ unitary matrix
$\boldsymbol{\mathsf{U}}$. Without loss of generality, we may suppose that the
initial state is of the form

$\boldsymbol{a}=\alpha\begin{bmatrix}1\\\ 0\\\ \vdots\\\
0\end{bmatrix}+\frac{1}{\sqrt{2}}\begin{bmatrix}z_{1}\\\ z_{2}\\\ \vdots\\\
z_{d}\end{bmatrix}=\alpha\,\boldsymbol{\psi}+\frac{\boldsymbol{z}}{\sqrt{2}}\;,$
(77)

where $\boldsymbol{z}$ is a $d$-dimensional vector of iid standard complex
Gaussian random variables and $\alpha\in\mathbb{C}$. Following the
transformation, the new state is

$\boldsymbol{a}^{\prime}=\alpha\,\boldsymbol{\mathsf{U}}\boldsymbol{\psi}+\frac{\boldsymbol{z}^{\prime}}{\sqrt{2}}\;,$
(78)

where $\boldsymbol{z}^{\prime}=\boldsymbol{\mathsf{U}}\boldsymbol{z}$ is,
again, a $d$-dimensional vector of iid standard complex Gaussian random
variables, owing to the unitarity of $\boldsymbol{\mathsf{U}}$.

Detection measurements on the $d$ modes will result in one of $2^{d}$ possible
outcomes. Let $(n_{1},\ldots,n_{d})\in\\{0,1\\}^{d}$ denote the outcome in
which mode 1 has $n_{1}$ detections, mode 2 has $n_{2}$, etc., and let
$P(n_{1},\ldots,n_{d})$ denote the probability of this outcome occurring.
Since the random variables $z^{\prime}_{1},\ldots,z^{\prime}_{d}$ are
statistically independent, this probability is given by

$P(n_{1},\ldots,n_{d})=\prod_{i=1}^{d}q_{i}^{n_{i}}(1-q_{i})^{1-n_{i}}\;,$
(79)

where $q_{i}=Q_{1}(2|\alpha\psi^{\prime}_{i}|,2\gamma)$ is the probability of
a threshold crossing event for mode $i$. (We assume, for simplicity, that all
detectors have the same threshold.)

We will be most particularly concerned with single-detection events (i.e.,
those for which $n_{1}+\cdots+n_{d}=1$), as these would be interpreted as
single-photon detections. Although such events occur with vanishingly small
probability as $|\alpha|$ becomes large (and low probability for $|\alpha|$
small), we may condition, via post-selection, on only such events and thereby
obtain a nonvanishing probability. Specifically, let $p_{i}$ denote the
probability that a single-detection event occurs on mode $i$, given that a
single-detection event occurs on any one mode. It follows that

$\begin{split}p_{i}&=q_{i}\prod_{j\neq
i}(1-q_{j})\left[\sum_{k=1}^{d}q_{k}\prod_{\ell\neq
k}(1-q_{\ell})\right]^{-1}=\frac{q_{i}}{1-q_{i}}\left[\sum_{k=1}^{d}\frac{q_{k}}{1-q_{k}}\right]^{-1}\;,\end{split}$
(80)

provided that $q_{k}\neq 1$ for all $k$. Note that, if there exists an $i$
such that $|\psi_{i}|>|\psi_{j}|$ for all $j\neq i$, then $p_{i}\to 1$ as
$|\alpha|\to\infty$. In other words, for bright light only the most probable
mode will have a single detection. For states with no unique maximum, the
asymptotic probability is spread uniformly amongst the maxima. The latter case
is consistent with the quantum mechanical predictions, while the former is
not. The right correspondence with quantum mechanics is then to be expected
for small or intermediate values of $|\alpha|$.

To examine the validity of our model, we performed linear quantum state
tomography (QST) on a random sample of pure states formed by applying unitary
matrices drawn from a Haar measure [45]. For a given transformed state
$\boldsymbol{\psi}$, QST was performed using a complete set of $d$-dimensional
Hermitian basis matrices
$\boldsymbol{\mathsf{B}}_{1},\ldots,\boldsymbol{\mathsf{B}}_{d^{2}}$ that are
orthonormal in the Hilbert-Schmidt inner product. Each
$\boldsymbol{\mathsf{B}}_{k}$ is diagonalized by a unitary matrix
$\boldsymbol{\mathsf{U}}_{k}$ such that
$(\boldsymbol{\mathsf{U}}_{k}^{\dagger}\boldsymbol{\mathsf{B}}_{k}\boldsymbol{\mathsf{U}}_{k})_{ij}=\beta_{ki}\delta_{ij}$.
To measure in this basis, we therefore transformed $\boldsymbol{\psi}$ to
$\boldsymbol{\psi}^{\prime}=\boldsymbol{\mathsf{U}}_{k}^{\dagger}\boldsymbol{\psi}$
and computed $p_{i}$ according to equation (80). The expectation value of
$\boldsymbol{\mathsf{B}}_{k}$ was taken to be
$p_{1}\beta_{k1}+\cdots+p_{d}\beta_{k1}$, so the inferred density matrix for
$\boldsymbol{\psi}$ is defined to be

$\boldsymbol{\mathsf{\rho}}=\sum_{k=1}^{d^{2}}\boldsymbol{\mathsf{B}}_{k}\,\sum_{i=1}^{d}p_{i}\,\beta_{ki}\;.$
(81)

Since the basis matrices are such that
$\mathrm{Tr}[\boldsymbol{\mathsf{B}}_{1}^{\dagger}\boldsymbol{\mathsf{B}}_{1}]=1$
and
$\mathrm{Tr}[\boldsymbol{\mathsf{B}}_{k}^{\dagger}\boldsymbol{\mathsf{B}}_{k}]=0$
for $k>1$, the trace
$\mathrm{Tr}[\boldsymbol{\mathsf{\rho}}\boldsymbol{\mathsf{B}}_{k}]$ gives the
expectation value of $\boldsymbol{\mathsf{B}}_{k}$, as defined above. This
allows us to identify $\boldsymbol{\mathsf{\rho}}$ as playing the role of a
quantum mechanical density operator.

With $\boldsymbol{\mathsf{\rho}}$ so computed, we examined the fidelity,
defined by the vector inner product

$F=\bra{\boldsymbol{\psi}}\boldsymbol{\mathsf{\rho}}\ket{\boldsymbol{\psi}}=\sum_{ij}\psi_{i}^{*}\rho_{ij}\psi_{j}\;,$
(82)

as a function of $d$, $|\alpha|$, and $\gamma$ over an ensemble of pure states
$\boldsymbol{\psi}$. In figure 8 we have plotted $F$ versus $|\alpha|$ for
$d=4$ and $\gamma=1$ for an ensemble of 30 randomly drawn pure states. For
this case, tensor products of the Pauli matrices were used for the orthonormal
basis. We observe that $F=1/d$ (corresponding to $p_{i}=1/d$) for
$|\alpha|=0$, as expected for pure vacuum noise. As $|\alpha|$ increases, $F$
increases monotonically to a value near unity. However, for sufficiently large
values of $|\alpha|$ the inferred density matrix acquires negative eigenvalues
and becomes invalid. For general quantum states, taking
$|\alpha|\sim\gamma\sqrt{d/2}$ tends to give near unity fidelity, albeit with
invalid density matrices. For “classical” states (i.e., those for which
$|\psi_{i}|=1$ for exactly one index $i$) the density matrix remains valid for
all $|\alpha|$ and the fidelity asymptotically approaches unity as
$|\alpha|\to\infty$. Qualitatively similar behavior is found when $\gamma$ is
varied while holding $\alpha$ fixed.

Figure 8: Plot of the fidelity of the QST solution versus the state amplitude
$|\alpha|$ for $\gamma=1$ over an ensemble of 30 pure states with $d=4$. The
red portion of the curves indicates where one or more eigenvalues in the
density matrix are negative.

Negative eigenvalues in density matrices obtained through linear state
tomography are a common occurrence in experimental quantum optics,
particularly for low-entropy, high fidelity states. Their presence might be
interpreted as an observed deviation from the Born rule, but they are more
commonly ascribed to mere “experimental inaccuracies and statistical
fluctuations” [46]. According to our model, such results are an inevitable
consequence of the parameter regime investigated and the data analysis methods
used to infer the quantum state. Since we compute the probabilities exactly to
perform state tomography, we may also conclude that the potential for negative
eigenvalues is an intrinsic property of the model and not one due to mere
sampling errors.

To address the problem of invalid density matrices obtained from linear state
tomography, maximum-likelihood estimation (MLE) methods are often used [47,
48, 49]. In this approach, one parameterizes a general, positive semi-definite
density matrix and estimates the parameters of this matrix using an
optimization scheme based on an assumption of Gaussian errors. By
construction, this approach always yields a valid density matrix. We
reexamined our results using the MLE-based state tomography tools provided by
the Kwiat Quantum Information Group [50]. The results of five randomly sampled
states with $d=4$ are shown in figure 9. For $|\alpha|$ less than or close to
unity, the MLE results agree with the previous linear tomography results.
However, for larger values of $|\alpha|$, the curves peak near unity and then
slowly decrease as we approach the classical regime of $|\alpha|\gg 1$.
Qualitatively similar behavior is found when $\gamma$ is varied while holding
$\alpha$ fixed. This shows that it is possible to infer density matrix
estimates from our model that are both valid and of high fidelity.

Figure 9: Plot of the fidelity of the QST solution using the MLE method versus
the state amplitude $|\alpha|$ for $\gamma=1$ over an ensemble of five pure
states with $d=4$.

The density matrix derived from QST may also be used to examine entanglement.
According to the Peres-Horodecki positive partial transpose (PPT) criterion, a
density operator that acts on a tensor product Hilbert space
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ will be separable with respect to
$\mathcal{H}_{A}$ and $\mathcal{H}_{B}$ if all the eigenvalues of its partial
transpose are positive [51]. In our case, $\otimes$ is the Kronecker product,
$\mathcal{H}_{A}=\mathbb{C}^{d_{A}}$, and
$\mathcal{H}_{B}=\mathbb{C}^{d_{B}}$, for $d_{A},d_{B}\in\mathbb{N}$. If we
write the density matrix $\boldsymbol{\mathsf{\rho}}$ as

$\boldsymbol{\mathsf{\rho}}=\sum_{ij}\sum_{k\ell}\rho_{ij;k\ell}\>\boldsymbol{e}^{A}_{i}(\boldsymbol{e}^{A}_{j})^{\dagger}\otimes\boldsymbol{e}^{B}_{k}(\boldsymbol{e}^{B}_{\ell})^{\dagger}\;,$
(83)

where $\boldsymbol{e}^{A}_{i},\boldsymbol{e}^{A}_{j}$ and
$\boldsymbol{e}^{B}_{k},\boldsymbol{e}^{B}_{\ell}$ are the standard unit
vectors in $\mathbb{C}^{d_{A}}$ and $\mathbb{C}^{d_{B}}$, respectively, then
the partial transpose with respect to $\mathcal{H}_{B}$ is

$\boldsymbol{\mathsf{\rho}}^{\mathsf{T}_{B}}=\sum_{ij}\sum_{k\ell}\rho_{ij;k\ell}\>\boldsymbol{e}^{A}_{i}(\boldsymbol{e}^{A}_{j})^{\dagger}\otimes\boldsymbol{e}^{B}_{\ell}(\boldsymbol{e}^{B}_{k})^{\dagger}\;.$
(84)

Negative eigenvalues of the partial transpose are a necessary, though not
sufficient, condition for the density matrix to be nonseparable (i.e.,
entangled). For certain cases, such as $d_{A}=d_{B}=2$, this condition is also
sufficient and therefore may be used as an entanglement witness [52]. In
figure 10 we have plotted the minimum eigenvalue of the partial transpose for
$d_{A}=d_{B}=2$ as a function of $|\alpha|$ for a maximally entangled Bell
state using a detection threshold of $\gamma=1$ and the aforementioned MLE
method to infer the quantum state. It is perhaps surprising that, although our
inferred density matrix does not have perfect fidelity, it is nevertheless
entangled (i.e., nonseparable), as witnessed by the negative eigenvalues of
the partial transpose for values of $|\alpha|$ above 0.6. The behavior for
large $|\alpha|$ shows an asymptotic approach to $-0.5$, the value predicted
by quantum mechanics for an ideal Bell state. Qualitatively similar behavior
is found when $\gamma$ is varied while holding $\alpha$ fixed.

Figure 10: Plot of the PPT entanglement witness for $d=4$ versus the state
amplitude $|\alpha|$ for a maximally entangled Bell state using a detection
threshold of $\gamma=1$.

In the above examples we have considered fidelity as a function of $|\alpha|$;
the dependency on both $|\alpha|$ and $\gamma$ is more subtle. In figure 11 we
have plotted the fidelity of the QST solution using the MLE method as a
function of both $|\alpha|$ and $\gamma$, averaged over an ensemble of 100
pure states with $d=4$. For any given $\gamma>0$, we see that there is a
unique value of $|\alpha|$ giving locally optimal fidelity. Generally, small
values of $|\alpha|$ and large values of $\gamma$ give good, albeit imperfect,
fidelity. The globally optimal fidelity is found to be about 0.98 and occurs
near $|\alpha|=1.2$ and $\gamma=1.5$. The fidelity for some states can be
higher and qualitatively different. For example, states with a high degree of
symmetry, such as the Bell states, can exhibit fidelities approaching unity
when either $|\alpha|$ or $\gamma$ is large.

Figure 11: Contour plot of fidelity using the MLE method versus the state
amplitude $|\alpha|$ and detection threshold $\gamma$ over an ensemble of 100
pure states with $d=4$. The black dot indicates a maximum fidelity of $0.98$
for $|\alpha|=1.2$ and $\gamma=1.5$. The black curve is a spline-interpolation
of the local maxima for each given value of $\gamma$.

In addition to fidelity, we also considered visibility as a metric of
agreement with quantum predictions. Visibility converges monotonically to one
as either $|\alpha|$ or $\gamma$ tends to infinity, but optimal fidelity
occurs only for finite values of these parameters, so optimally satisfying
both may not be possible. We have plotted the visibility as a function of
$|\alpha|$ and $\gamma$ in figure 12. The point of maximum fidelity, from
figure 11, corresponds to a visibility of 0.94. Taking $\gamma$ to be large
and $|\alpha$, suitably chosen, to be small, can achieve higher visibility
while still maintaining a high fidelity, but the two metrics cannot be made
arbitrarily close to unity. Experimentally, fidelities and visibilities as
high as 0.99 have been observed, though not necessarily in the same context. A
more detailed analysis of specific experiments would therefore be needed for a
proper comparison.

Figure 12: Plot of visibility versus the state amplitude $|\alpha|$ and
detection threshold $\gamma$ for $d=4$. The black dot at $|\alpha|=1.2$ and
$\gamma=1.5$ corresponds to a visibility of 0.94.

##  5 Conclusion

Assuming a classical zero-point field and deterministic threshold detectors,
we have shown that one is able to reproduce many of the experimentally
observed phenomena attributed to single photons and thought to be uniquely
quantum in nature. In so doing we have established that such phenomena do, in
fact, have classical analogues. Weak coherent light in combination with a
reified zero-point field considered in the single-mode regime are found to
give probabilistic outcomes that are in close agreement with the Born rule for
single-photon, multi-mode states when post-selection on single-detection
events is performed. This agreement was verified explicitly by performing
quantum state tomography and computing the fidelity of the resulting density
matrix. The model results are not always in perfect agreement with the
idealized quantum mechanical predictions, but they are largely consistent with
experimental observations and data analysis methods in the appropriate
parameter regimes. Deviations are, however, expected in regimes in which
either the amplitude of the light or the threshold of detection is large. The
best overall agreement appears to be for a set of parameter values in which
the amplitude is small and threshold is large. This model therefore provides a
local, realistic picture of wave/particle duality and single-photon
entanglement that is grounded in a physical and wholly classical model. A
similar classical description of homodyne measurements, temporal behavior, and
multi-photon entanglement are left for future work.

This work was supported by the U. S. Office of Naval Research under Grant No.
N00014-17-1-2107. The authors also would like to thank the reviewers for their
thoughtful comments and suggestions.

## References

  * [1] J. von Neumann.  Mathematische Grundlagen der Quantentheorie.  Springer, 1931.
  * [2] J. S. Bell.  On the Einstein Podolsky Rosen paradox.  Physics, 1:195, 1964.
  * [3] W. M. de Muynck.  Foundations of Quantum Mechanics, An Empiricist Approach.  Kluwer Academic Publishers, 2002.
  * [4] N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner.  Bell nonlocality.  Reviews of Modern Physics, 2014:419, 2014.  https://doi.org/10.1103/RevModPhys.86.419.
  * [5] B. Hensen et al.  Loophole-free bell inequality violation using electron spins separated by 1.3 kilometres.  Nature, 526:682, 2015.  https://doi.org/10.1038/nature15759.
  * [6] M. Giustina et al.  Significant-loophole-free test of Bell’s theorem with entangled photons.  Physical Review Letters, 115:250401, 2015.  https://doi.org/10.1103/PhysRevLett.115.250401.
  * [7] L. K. Shalm et al.  Strong loophole-free test of local realism.  Physical Review Letters, 115:250402, 2015.  https://doi.org/10.1103/PhysRevLett.115.250402.
  * [8] W. Rosenfeld et al.  Event-ready bell test using entangled atoms simultaneously closing detection and locality loopholes.  Physical Review Letters, 119:010402, 2017.  https://doi.org/10.1103/PhysRevLett.119.010402.
  * [9] P. Bierhorst et al.  Experimentally generated randomness certified by the impossibility of superluminal signals.  Nature, 556:223, 2018.  https://doi.org/10.1038/s41586-018-0019-0.
  * [10] V. B. Berestetskii, L. P. Pitaevskii, and E. M. Lifshitz.  Quantum Electrodynamics, volume 4.  Elsevier, 2nd edition, 1982.
  * [11] P. W. Milonni.  The Quantum Vacuum. An Introduction to Quantum Electrodynamics.  Academic Press, 1994.
  * [12] N. P. Landsman.  Compendium of Quantum Physics, chapter The Born rule and its interpretation.  Springer, 2008.
  * [13] M. Born.  Zur quantenmechanik der stoßvorgänge.  Zeitschrift für Physik, 37:863, 1926.  https://doi.org/10.1007/BF01397477.
  * [14] A. M. Gleason.  Measures on the closed subspaces of a Hilbert space.  Journal of Mathematical Mechanics, 6:885, 1957.  https://doi.org/10.1007/978-94-010-1795-4_7.
  * [15] D. Deutsch.  Quantum theory of probability and decisions.  Proceedings of the Royal Society A, 455:3129, 1999.  https://doi.org/10.1098/rspa.1999.0443.
  * [16] H. Barnum, C. M. Caves, J. Finkelstein, C. A. Fuchs, and R. Schack.  Quantum probability from decision theory?  Proceedings of The Royal Society, A456:1175, 2000.  https://doi.org/10.1098/rspa.2000.0557.
  * [17] W. H. Zurek.  Probabilities from entanglement, Born’s rule $p_{k}=|\psi_{k}|^{2}$ from envariance.  Physical Review A, 71:052105, 2005.  https://doi.org/10.1103/PhysRevA.71.052105.
  * [18] M. Schlosshauer and A. Fine.  On Zurek’s derivation of the Born rule.  Foundations of Physics, 35:197, 2005.  https://doi.org/10.1007/s10701-004-1941-6.
  * [19] L. Masanes, T. D. Galley, and M. P. Müller.  The measurement postulates of quantum mechanics are operationally redundant.  Nature Communications, 10:1361, 2019.  https://doi.org/10.1038/s41467-019-09348-x.
  * [20] A. E. Allahverdyan, R. Balian, and T. M. Nieuwenhuizen.  Understanding quantum measurement from the solution of dynamical models.  Physics Reports, 525:1, 2013.  https://doi.org/10.1016/j.physrep.2012.11.001.
  * [21] L. de la Peña and A. M. Cetto.  The Quantum Dice: An Introduction to Stochastic Electrodynamics.  Kluwer, 1995.
  * [22] A. Casado, T. W. Marshall, and E. Santos.  Parametric downconversion experiments in the Wigner representation.  Journal of the Optical Society of America B, 14:494, 1997\.  https://doi.org/10.1364/JOSAB.14.000494.
  * [23] T. W. Marshall and E. Santos.  Stochastic optics: A reaffirmation of the wave nature of light.  Foundations of Physics, 18:185, 1988.  https://doi.org/10.1007/BF01882931.
  * [24] G. Adenier.  Violation of Bell inequalities as a violation of fair sampling in threshold detectors.  In AIP Conference Proceedings 1101, page 8, 2009.  https://doi.org/10.1063/1.3109977.
  * [25] B. La Cour.  A locally deterministic, detector-based model of quantum measurement.  Foundations of Physics, 44:1059, 2014.  https://doi.org/10.1007/s10701-014-9829-6.
  * [26] A. Khrennikov.  Beyond Quantum.  Pan Stanford Publishing, 2014.
  * [27] M. Planck.  Eine neue strahlungshypothese.  Verhandlungen der Deutschen Physikalischen Gesellschaft, 13:138, 1911.
  * [28] T. W. Marshall.  Random electrodynamics.  Proceedings of the Royal Society, A276:475, 1963.  http://doi.org/10.1098/rspa.1963.0220.
  * [29] M. Ibison and B. Haisch.  Quantum and classical statistics of the electromagnetic zero-point field.  Physical Review A, 54:2737, 1996.  https://doi.org/10.1103/PhysRevA.54.2737.
  * [30] A. Lasota and M. Mackey.  Chaos, Fractals, and Noise: Stochastic Aspects of Dynamics.  Springer, 2nd edition, 1998.
  * [31] H. M. França and T. W. Marshall.  Excited states in stochastic electrodynamics.  Physical Review A, 38:3258, 1988.  https://doi.org/10.1103/PhysRevA.38.3258.
  * [32] M. Ossiander et al.  Absolute timing of the photoelectric effect.  Nature, 561:374, 2018.  https://doi.org/10.1038/s41586-018-0503-6.
  * [33] H. Shibata, K. Shimizu, H. Takesue, and Y. Tokura.  Ultimate low system dark count rate for superconducting nanowire single-photon detector.  Optics Letters, 40:3428, 2015.  https://doi.org/10.1364/OL.40.003428.
  * [34] N. L. Johnson, S. Kotz, and N. Balakrishnan.  Continuous Univariate Distributions.  John Wiley and Sons, 1994.
  * [35] K. Cahill and R. Glauber.  Density operators and quasiprobability distributions.  Physical Review, 177:1882, 1969.  https://doi.org/10.1103/PhysRev.177.1882.
  * [36] K. Fujii and T. Suzuki.  A new symmetric expression of Weyl ordering.  Modern Physics Letters A, 19:827, 2004.  https://doi.org/10.1142/S021773230401374X.
  * [37] B. F. Levine, D. G. Bethea, and J. C. Campbell.  Near room temperature 1.3um single photon counting with a ingaas avalanche photodiode.  Electronics Letters, 20:596, 1984.  https://doi.org/10.1049/el:19840411.
  * [38] J. Oh, C. Anonelli, M. Tur, and M. Brodsky.  Method for characterizing single photon detectors in saturation regime by cw laser.  Optics Express, 18:5906, 2010.  https://doi.org/10.1364/OE.18.005906.
  * [39] R. Loudon.  Non-classical effects in the statistical properties of light.  Reports on Progress in Physics, 43:913, 1980.  https://doi.org/10.1088/0034-4885/43/7/002.
  * [40] R. W. Boyd, S. G. Lukishova, and V. N. Zadkov, editors.  The First Single Photon Sources and Single Photon Interference Experiments, chapter 1.  Springer, 2019.
  * [41] P. Grangier, G. Roger, and A. Aspect.  Experimental evidence for a photon anticorrelation effect on a beam splitter: A new light on single-photon interferences.  Europhysics Letters, 1:173, 1986.  https://doi.org/10.1209/0295-5075/1/4/004.
  * [42] J. J. Thorn et al.  Observing the quantum behavior of light in an undergraduate laboratory.  American Journal of Physics, 79:1210, 2004.  https://doi.org/10.1119/1.1737397.
  * [43] P. Kwiat and H. Weinfurter.  Embedded Bell-state analysis.  Physical Review A, 58:R2623, 1998.  https://doi.org/10.1103/PhysRevA.58.R2623.
  * [44] V. Jacques, E. Wu, F. Grosshans, F. Treussart, P. Grangier, A. Aspect, and J.-F. Roch.  Experimental realization of Wheeler’s delayed-choice gedanken experiment.  Science, 315:966, 2007.  https://doi.org/10.1126/science.1136303.
  * [45] F. Mezzadri.  How to generate random matrices from the classical compact groups.  Notices of the AMS, 54:592, 2007.  http://www.ams.org/notices/200705/index.html.
  * [46] D. F. V. James, P. G. Kwiat, W. J. Munro, and A. G. White.  Measurement of qubits.  Physical Review A, 64:052312, 2001.  https://doi.org/10.1103/PhysRevA.64.052312.
  * [47] Z. Hradil.  Quantum-state estimation.  Physical Review A, 55:R1561(R), 1997.  https://doi.org/10.1103/PhysRevA.55.R1561.
  * [48] Z. Hradil, J. Summhammer, G. Badurek, and H. Rauch.  Reconstruction of the spin state.  Physical Review A, 62:014101, 2000.  https://doi.org/10.1103/PhysRevA.62.014101.
  * [49] J. B. Altepeter, D. F. V. James, and P. G. Kwiat.  Lecture Notes in Physics, chapter Quantum State Estimation.  Springer, Berlin, 2004.
  * [50] Kwiat Quantum Information Group.  Guide to Quantum State Tomography.  http://research.physics.illinois.edu/QI/Photonics/Tomography/ Accessed 20 December 2019.
  * [51] A. Peres.  Separability criterion for density matrices.  Physical Review Letters, 77:1413, 1996.  https://doi.org/10.1103/PhysRevLett.77.1413.
  * [52] P. Horodecki M. Horodecki and R. Horodecki.  Separability of mixed states: necessary and sufficient conditions.  Physics Letters A, 223:1, 1996.  https://doi.org/10.1016/S0375-9601(96)00706-2.