# Graviton detection and the quantization of gravity

Daniel Carney<EMAIL_ADDRESS>Physics Division, Lawrence Berkeley National
Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA  Valerie Domcke
<EMAIL_ADDRESS>Theoretical Physics Department, CERN, 1 Esplanade des
Particules, CH-1211 Geneva 23, Switzerland  Nicholas L. Rodd
<EMAIL_ADDRESS>Theoretical Physics Department, CERN, 1 Esplanade des
Particules, CH-1211 Geneva 23, Switzerland

###### Abstract

We revisit a question asked by Dyson: “Is a graviton detectable?” We
demonstrate that in both Dyson’s original sense and in a more modern
measurement-theoretic sense, it is possible to construct a detector sensitive
to single gravitons, and in fact a variety of existing and near-term
gravitational wave detectors can achieve this. However, while such a signal
would be consistent with the quantization of the gravitational field, we draw
on results from quantum optics to show how the same signal could just as well
be explained via classical gravitational waves. We outline the kind of
measurements that would be needed to demonstrate quantization of gravitational
radiation and explain why these are substantially more difficult than simply
counting graviton clicks or observing gravitational noise in an
interferometer, and likely impossible to perform in practice.

††preprint: CERN-TH-2023-155

A wide array of detector architectures which are sensitive to single photons
now exist [1, 2]. On the contrary, designing a detector capable of detecting
single gravitons may seem like an insurmountable challenge. Dyson, and
separately Rothman and Boughn, have argued that this would be fundamentally
impossible in certain architectures, but possible with others in a sense we
make precise below [3, 4, 5]. Supposing that one could construct such a
detector, an important question arises: would observation of a signal
attributable to single gravitons imply that the gravitational field is
quantized?

First, let us motivate the idea that single graviton detection is in principle
possible, following Dyson. The energy density in a gravitational wave is
$\rho=\tfrac{1}{4}M_{\rm
Pl}^{2}\langle\partial_{t}h^{\mu\nu}\partial_{t}h_{\mu\nu}\rangle$, with
$M_{\rm Pl}$ the reduced Planck mass [6].111We take $\hbar=c=1$ and signature
$(-,+,+,+)$ throughout. For a wave of strain $h$ and frequency $\omega$, this
corresponds to $\rho=\tfrac{1}{4}h^{2}\omega^{2}M_{\rm Pl}^{2}$. Dividing this
energy up into gravitons of energy $\omega$, we find that the number of quanta
per de Broglie volume $\lambda_{\rm dB}^{3}=f^{-3}$ is

$n\lambda_{\rm dB}^{3}=\frac{\pi h^{2}M_{\rm Pl}^{2}}{2f^{2}}\simeq 2\times
10^{35}\left(\frac{h}{10^{-22}}\right)^{2}\left(\frac{1~{}{\rm
kHz}}{f}\right)^{2}\\!,$ (1)

with $f=\omega/2\pi$ is the linear frequency. These parameters are benchmarked
to a LIGO event, where it is clear that the number density of gravitons in an
observable wave is astronomically large. However, it is possible to look for
gravitons at much higher frequencies, where this occupation density is
diluted.

Figure 1: The parameter space for gravitational radiation at high frequency,
and the region where the field – if quantized – would be better described as a
highly dilute state of graviton particles rather than a gravitational wave
[$n\lambda_{\rm dB}^{3}<1$ from Eq. (1)]. Instruments such as CAST are already
reaching into this parameter space.

For example, consider the CERN Axion Solar Telescope (CAST) [7]. This device
consists of a large permanent magnet $B\sim 10~{}{\rm T}$ with an X-ray
photodetector. The coupling
${\mathcal{L}}=\tfrac{1}{2}h_{\mu\nu}T^{\mu\nu}_{\rm EM}$ leads to linear
conversion between gravitational and electromagnetic radiation in the presence
of the magnetic field [8, 9, 10, 11, 12]. Given incident gravitational
radiation, a detected X-ray photon implies the absorption of an individual
graviton, under the assumption that the gravitational field is quantized.
Moreover, at X-ray frequencies $f\sim 10^{18}~{}{\rm Hz}$ and with CAST’s
nominal sensitivity $h\sim 10^{-27}$ [10], this detector is sensitive to
signals in the regime where the number quanta per de Broglie cell can become
less than one. In such a state, the distribution of gravitons is highly
dilute,222The analogue of Eq. (1) for dark matter is $n\lambda_{\rm
dB}^{3}=(2\pi)^{3}\rho/m^{4}v^{3}$, so that locally dark matter undergoes the
analogous transition at $m_{\scriptscriptstyle\textrm{DM}}\sim 10~{}{\rm eV}$.
For a review, see Ref. [13]. and in this sense CAST can detect individual
gravitons.333We know of no sufficiently strong sources for gravitational
radiation at this frequency (see App. A). In this paper, we will only be
concerned with questions of principle involving detector sensitivities.

In Fig. 1 we show that the region where $n\lambda_{\rm dB}^{3}<1$ covers a
large portion of the parameter space for high frequency gravitational
radiation, and that CAST and its successor IAXO [14] clearly cut into that
space. The remaining aspects of the figure largely follow Ref. [15], and we
refer there for details, see also Refs. [16, 17, 18, 11, 19, 20, 12, 21, 22,
23, 24, 25] for recent work. We emphasize that for transient signals, the
strain must be interpreted carefully as discussed in App. A; in particular,
the figure assumes a signal duration comparable to or longer than the
measurement time, which for CAST is roughly a year.

While an intensity detector like CAST could be sensitive to such single
graviton signals, it has also been suggested that linear gravitational wave
detectors, for example LIGO, could be sensitive to analogous “shot noise”
signals due to gravitational quantization in highly squeezed states [26, 27,
28, 29, 30]. Supposing that such a signal was observed in either case, a
natural question to ask is whether it would constitute experimental proof that
the gravitational field is quantized.444Assuming the events are not attributed
to a background event or an axion. For example, in CAST, a significant
background source is cosmic-ray muons interacting with the detector material
and generating fluorescent X-rays [7, 14]. The answer is no: classical
gravitational radiation can equally well explain the observed detector
response.

To understand why, let us review a classic argument from quantum optics,
following Ref. [31]. Consider a photodetector operating based on the
photoelectric effect, such as a collection of approximately free electrons in
a semiconductor. A fully quantized Hamiltonian describing a given detector
electron coupled to the incident electromagnetic field is

$\hat{H}_{\rm Q}=\hat{H}_{\rm
det}+\hat{H}_{A}-\frac{e}{m}\hat{{\mathbf{p}}}\cdot\hat{{\mathbf{A}}},$ (2)

where the first term describes a ground state and conduction band separated by
an energy gap $\Delta>0$, $\hat{{\mathbf{p}}}$ is the momentum operator of an
electron, $\hat{{\mathbf{A}}}$ is the gauge field operator, and $\hat{H}_{A}$
is its kinetic energy. The probability for a single electron to be excited
from the ground state into the conduction band in a short time $\delta t$
computed from Eq. (2) is,

$P(t)\delta t\simeq\eta\braket{I(t)}\Theta(\omega-\Delta)\delta t.$ (3)

Here, $\eta$ is a constant representing the detector efficiency, and we assume
the incident field is nearly monochromatic with frequency $\omega$ and a
slowly varying intensity $I(t)$. The result exhibits the hallmarks of the
photoelectric effect: zero emission for $\omega<\Delta$, and emission at a
rate proportional to the light intensity beyond this. Moreover, the individual
events can be resolved as discrete “clicks”, occuring when an individual
electron is excited.

In contrast, consider a semi-classical model, where the gauge field is _not_
quantized, but viewed as a classical field:

$\hat{H}_{\rm SC}=\hat{H}_{\rm
det}+H_{A}-\frac{e}{m}\hat{{\mathbf{p}}}\cdot{\mathbf{A}}.$ (4)

Here ${\mathbf{A}}$ is a time dependent c-number that in principle can be
stochastic, drawn from a classical distribution $P_{\rm cl}({\mathbf{A}})$.
The detector electrons remain quantized, so our experiment will continue to
register a series of discrete clicks. In fact, a straightforward calculation
in perturbation theory (again see Ref. [31] or App. B) yields precisely the
same result as Eq. (3) for the rate of electron excitations. In this latter
case, the bracket in $\braket{I(t)}$ is an average over the classical
distribution $P_{\rm cl}({\mathbf{A}})$.555The presence of the energy
conservation $\Theta(\omega-\Delta)$ may seem surprising. In the semi-
classical calculation, this arises from a factor $\left|\int_{0}^{\delta
t}dt\,e^{i(E+\Delta-\omega)t}\right|^{2}\sim\delta(E+\Delta-\omega)$ in the
perturbative transition amplitude for a given electron to be excited into a
conduction state of energy $E$ by an incident plane wave, ${\mathbf{A}}\propto
e^{-i\omega t}$. For $\omega\delta t\gg 1$, certainly the case for example in
optical detection with $\omega\gtrsim 10^{15}~{}{\rm Hz}$, the approximation
as a delta function is excellent. More details are provided in App. B. We
conclude that the photoelectric effect can be equally well explained by a
classical electromagnetic field or quantized photons.

This leads to a final question: what kind of data stream could we observe in
the detector output that can be explained by a quantum model of the radiation
like Eq. (2) but _not_ by a classical model like Eq. (4)? A famous example in
quantum optics is sub-Poisson counting statistics, in which the variance in
the observed photon or graviton count rate is smaller than the mean rate [31,
32]. Heuristically, if the radiation were quantized, one could send a stream
of photons with fixed separation between them, leading to zero variance
amongst the arrival times at the detector, and for an ideal detector, zero
variance in the measured response. As we review in Sec. I, sub-Poisson
counting statistics can be explained by a quantized radiation field prepared
in a variety of non-trivial quantum states (e.g., a squeezed state or number
eigenstate), but cannot be explained by a classical field. Thus, observation
of such a signal would constitute proof that the incident field is quantized.

However, as we will emphasize, the requirements for observing such a
quantization signature for gravitational radiation are very stringent. We
first need some source of non-classical graviton states, and then we need a
detector which can see the subtle quantum effects that differentiate this
state from a classical ensemble. In particular, simply counting gravitons, or
simply observing gravitational shot noise (including that produced by a highly
squeezed state [26, 27, 28, 29, 30]) is not enough to demonstrate
quantization. And even if a source for such states exists in nature, it is
_substantially_ more difficult to make a detector capable of teasing out the
subtle quantum parts of these signatures, and probably impossible in practice.

The rest of the paper is dedicated to making the above statements precise. We
begin in Sec. I by sketching a simple model of a gravitational radiation
detector for illustrative purposes. We then use this to discuss the various
quantization witnesses described above. In Sec. II we then discuss what kinds
of experimental parameters would be required to construct a detector sensitive
to the quantization measures. Finally in Sec. III we relate this to other
tests of quantization of gravity, for example involving the microwave
background or tabletop entanglement experiments.

##  I Graviton measurement theory

The gravitational analogue to our models of quantized and semiclassical
radiation coupled to a detector is immediate. We assume that spacetime is
nearly flat $g_{\mu\nu}\simeq\eta_{\mu\nu}$ near the detector and consider
small perturbations of the total metric $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$.
We will then consider the difference between a quantized gravitational field

$\hat{H}_{\rm Q}=\hat{H}_{\rm det}+\hat{H}_{h}-\frac{1}{2}\int
d^{3}{\mathbf{x}}\,\hat{h}_{\mu\nu}\hat{T}^{\mu\nu},$ (5)

where $\hat{T}^{\mu\nu}$ is some part of the matter stress tensor operator in
the detector (e.g., the mirrors in LIGO or the electromagnetic field in the
CAST tube), and the classical version

$\hat{H}_{\rm SC}=\hat{H}_{\rm det}+H_{h}-\frac{1}{2}\int
d^{3}{\mathbf{x}}\,h_{\mu\nu}\hat{T}^{\mu\nu},$ (6)

where again $h_{\mu\nu}$ is a classically random c-number field. The quantized
model in Eq. (5) should be, as usual, interpreted as an effective field
theory, valid at the low energy densities considered in this paper.

The goal is to determine how these hypotheses could be distinguished. We begin
in Sec. I.1 with a review of the theory of coherence and classicality
conditions in a radiation field. In Sec. I.2 we give a toy model for a
gravitational wave detector which exhibits the difference between an intensity
(number) detector and a linear (amplitude) detector. These two detectors offer
different observables which can be used to distinguish Eqs. (5) and (6), and
we study these in Sec. I.3 and I.4, respectively.

###  I.1 Classicality vs. quantization

The radiation fields can be decomposed as usual into modes, whether we are
dealing with classical or quantized radiation. We will find it simplest to
discretize the spectrum by placing everything in a large box of length $L$,
and take $L\to\infty$ at the end. In terms of positive and negative frequency
components,

$\displaystyle A_{\mu}(x)$
$\displaystyle=\sum_{{\mathbf{k}},s}\epsilon_{s,\mu}({\mathbf{k}})u_{{\mathbf{k}}}({\mathbf{x}})e^{i\omega_{{\mathbf{k}}}t}a_{{\mathbf{k}},s}^{\dagger}+{\rm
h.c.},$ (7) $\displaystyle h_{\mu\nu}(x)$ $\displaystyle=M_{\rm
Pl}^{-1}\sum_{{\mathbf{k}},s}\epsilon_{s,\mu\nu}({\mathbf{k}})v_{{\mathbf{k}}}({\mathbf{x}})e^{i\omega_{{\mathbf{k}}}t}b_{{\mathbf{k}},s}^{\dagger}+{\rm
h.c.},$

where the $\epsilon_{s}$ with $s=1,2$ form bases for the polarization vectors
and tensors, and the $u_{{\mathbf{k}}},\,v_{{\mathbf{k}}}$ form complete sets
of mode functions. In this section, everything will be treated symmetrically
between the electromagnetic and gravitational case. Since we are primarily
interested in gravity, we use the gravitational $b_{{\mathbf{k}},s}$ for the
Fourier components, but the same results all apply for electromagnetism $b\to
a$. See App. E for more details on these mode expansions.

The distinction between classical and quantum models comes from how we treat
the coefficients $b_{{\mathbf{k}},s}$. In the classical case, the
$b_{{\mathbf{k}},s}$ are c-number dynamical variables. A general state
consists of a probability distribution $P_{\rm cl}(\\{b_{{\mathbf{k}},s}\\})$,
which can for example include classical correlations between the different
modes. In the quantum case, the $b_{{\mathbf{k}},s}$ are promoted to
operators. A general state consists of a density matrix $\rho$, which can
include both classical correlations and quantum correlations (entanglement).
Our goal is to understand what kind of observations can distinguish between
these two cases.

For simplicity, we will focus on narrow-band signals in what follows; we
discuss generalizations in the appendices. Consider a single mode of the
radiation field of interest (gravitational or electromagnetic), with wave
vector ${\mathbf{k}}$, polarization $s$, and let $b=b_{{\mathbf{k}},s}$.666In
the case of a classical gravitational wave, $b$ translates to the strain $h$
as follows: the wave has fixed energy density $\rho=\omega^{2}M_{\rm
Pl}^{2}h^{2}/4=\omega\sum_{s}|b_{s}|^{2}/4L^{3}$, so for a single polarization
$b_{s}^{2}=\omega M_{\rm Pl}^{2}h_{s}^{2}L^{3}$. We can give a simple
description of the general quantum state of this mode which will make it easy
to compare to the classical case. Define the coherent states $\ket{\beta}$ for
this mode as usual, $b\ket{\beta}=\beta\ket{\beta}$, for any complex number
$\beta$. These form an overcomplete basis for the Hilbert space of the mode,
since they are non-orthogonal
$\braket{\beta}{\beta^{\prime}}=e^{-|\beta-\beta^{\prime}|^{2}/2}$. In fact,
one can express any density matrix of the mode as

$\rho_{\rm rad}=\int d\beta\,P(\beta)\ket{\beta}\bra{\beta}\\!,\ \ \int
d\beta\,P(\beta)=1.$ (8)

This representation of the density matrix is called the Glauber-Sudarshan (or
sometimes “Glauber P”) representation [33, 34]. The normalization condition on
$P(\beta)$ simply enforces ${\rm Tr}[\rho_{\rm rad}]=1$.

The representation Eq. (8) is highly suggestive: it looks like we can write
any state of the mode as a classical ensemble of coherent states. Indeed, in
the case that we _can_ write the radiation state as such an ensemble, one can
always reproduce a detector’s output with the purely classical radiation
model, namely an ensemble of classical plane waves with distribution
$P(\beta)$ [34, 35, 31]. This is true _even including_ observation of “vacuum
fluctuations” or shot noise, which can be explained by the quantum mechanical
nature of the detector even with a purely classical signal, as discussed in
the introduction.

However, in a general quantum state, the weight function $P(\beta)$ can take
negative values, or be a highly singular function.777By a theorem of Schwartz,
these singular $P$-functions must take negative values. For further details,
see App. C. In this case, $P(\beta)$ is not a classical probability
distribution. It is precisely in this case – when $P(\beta)$ fails to be a
classical distribution – that one can find observables in a detector output
that can be explained by a quantized field model like Eqs. (2) or (5), but
cannot explained by a classical field model like in Eqs. (4) or (6). We
therefore refer to such states as “non-classical.” Examples of such states
include number eigenstates, superpositions of coherent states such as
$\ket{\psi}\sim\ket{\beta_{1}}+\ket{\beta_{2}}$, and squeezed coherent states;
see App. C. Notice in particular that the vacuum is classical in this
definition, since it is the trivial ensemble consisting of just the zero-
amplitude coherent state.888As originally shown by Glauber [34], $P(\beta)$
obeys an analogue of the central limit theorem, implying that Gaussian
$P(\beta)$, which are positive definite, are highly generic. For example,
thermal radiation can be described by a Gaussian $P(\beta)$.

An illuminating example of such a non-classicality observable can be
constructed in the optics setting, where the radiation field can be
controlled. Consider an incoming beam with annihilation operator $b_{1}$
prepared in a coherent state $\ket{\beta}_{1}$. After passing this beam
through a 50-50 beam splitter with vacuum $\ket{0}_{2}$ sent through the other
port, the state of the two arms is

$\ket{\beta}\ket{0}\to\ket{\beta/\sqrt{2}}\ket{\beta/\sqrt{2}}\\!.$ (9)

For completeness, we provide the calculation in App. D. If we have a pair of
photodetectors, one for each outgoing beam, they will continue to each
register half of the intensity, even for arbitrarily small beam energy
$|\beta|^{2}\to 0$. On the other hand, consider something “non-classical” like
a beam of single-photon states, again with vacuum through the dark port:

$\ket{1}\ket{0}\to\frac{1}{\sqrt{2}}\left[\ket{1}\ket{0}+\ket{0}\ket{1}\right]\\!.$
(10)

Now the photodetectors will be perfectly anti-correlated, i.e., either one
detector registers the photon or the other does, at each detection event.
There is no classical random ensemble that can reproduce this result [34].
Note that the non-classical state produced an entangled two-body state, while
the coherent state did not; this is a general feature of non-classical states
[36]. While this example uses a beamsplitter and is thus not easy to
generalize to the gravitational setting, we will give analogous examples of
these kinds of measurements with realistic gravitational radiation detectors
in the next sections.

To summarize, in order to perform a measurement on a gravitational radiation
mode that would require a non-classical description and thus demonstrate
quantization of the radiation field, we have two basic requirements:

  1. 1.

Nature has to provide a non-classical state of the radiation mode; and

  2. 2.

We have to detect some signature of this non-classicality, for example,
correlation statistics as just described.

These are very stringent conditions. In this paper, we make no comments on the
possible sources of such states, although see [26, 27, 28, 29, 30] for some
ideas. Our focus will be on the second condition, which as we shall see in
what follows, entails extremely challenging detector parameters.

###  I.2 A simple detector model

In order to be precise about detecting gravitons and non-classicality
observables, we will need to discuss two important issues. The first is that
different detectors measure the field in different bases. As concrete
examples, CAST is an intensity detector (it clicks when a graviton is
absorbed), while lower-frequency axion haloscopes based on cavities, for
example ADMX [37] or HAYSTAC [38], could be operated as amplitude detectors
(they measure the amplitude of a cavity mode). Similarly, LIGO is an amplitude
detector. The second issue is that, in general, real detectors do not perform
projective measurements but more general, “weak” measurements on the incoming
radiation field. This is particularly important in the gravitational case,
where the incident field is extremely weakly coupled to the detector.

$a$$B_{0}$$\otimes$$b$${\mathbf{p}}_{0}$$a$$B_{0}$$\otimes$$b$$\omega_{0}$
Figure 2: Schematic implementation of our toy model detectors. The
gravitational mode $b$ is incident on a fixed magnetic field $B_{0}$. In this
field, $b$ can convert to an electromagnetic mode $a$, which can then be
detected. The caps at the end represent photon counters. Top: directly
counting the converted photons produces a gravitational intensity detector, as
in Eq. (16). Bottom: an external light beam at frequency $\omega_{0}$ and
known phase is interfered with the light coming out of the detector field.
This produces a homodyne (amplitude) measurement of the light amplitude $X$,
and thus an amplitude measurement of the gravitational signal, as in Eq. (20).

To make our discussion of both of these issues concrete, consider a simple
gravitational radiation detector consisting of a strong homogeneous magnetic
field $B_{0}$ along the $z$-axis, distributed over some spatial volume
$\ell^{3}$, as shown in Fig. 2. For conceptual simplicity, we will consider an
incoming “beam” of gravitational radiation along the $x$-axis with momentum
${\mathbf{k}}_{0}=k_{0}\hat{{\mathbf{x}}}$. The modes of the electromagnetic
field couple to the gravitational field through the Gertsenshtein effect
(i.e., the coupling ${\mathcal{L}}=\tfrac{1}{2}h_{\mu\nu}T^{\mu\nu}_{\rm EM}$)
as discussed in the introduction.999Note that this interaction is quite
general. In addition to its clear connection to devices like CAST, it can also
describe the basic detection mechanism of interferometers like LIGO. While
these are commonly described in terms of gravitational waves moving test
mirrors, they can equivalently be described in terms of gravitational waves
acting on the laser between the mirrors [39]. The interaction Hamiltonian, in
the Schrödinger picture, is

$\displaystyle V=$ $\displaystyle
B_{0}\int_{\ell^{3}}d^{3}{\mathbf{x}}\,\left(F_{yx}h_{yy}+F_{zx}h_{yz}\right)$
(11) $\displaystyle=$
$\displaystyle\frac{B_{0}}{2}\int_{\ell^{3}}d^{3}{\mathbf{x}}\,\left(\delta
B_{z}h_{+}-\delta B_{y}h_{\times}\right)e^{ik_{0}x}+\textrm{h.c.},$

where $\delta B_{y,z}$ denote the perturbations of the magnetic field strength
around the background field, and $h_{s}=b_{s}/M_{\rm Pl}\sqrt{\omega L^{3}}$
are the amplitudes of the individual polarizations, according to the
conventions we establish in App. E. We include this here to emphasize that
each polarization couples to a specific direction of the magnetic field
perturbations.

Suppose that we can monitor a particular mode ${\mathbf{p}}_{0},r_{0}$ of the
electromagnetic field, with annihilation operator
$a=a_{{\mathbf{p}}_{0},r_{0}}$ and frequency $\omega_{0}=|{\mathbf{p}}_{0}|$,
polarized with $B$ along the $y$-axis for definiteness. Again for simplicity,
we will assume that this mode is perfectly matched to the incoming
gravitational radiation ${\mathbf{p}}_{0}={\mathbf{k}}_{0}$. In this case, the
interaction Hamiltonian can be simplified to

$V=ig\big{(}a^{\dagger}b-ab^{\dagger}\big{)},\ \ \
g=\frac{B_{0}}{2\sqrt{2}M_{\rm Pl}}\sqrt{\frac{\ell^{3}}{L^{3}}}.$ (12)

The coupling $g$ has units of energy. The artificial spatial scale $L$
reflects the discretized radiation modes $b$; in the continuum limit these
scale like $L^{3/2}$ and one would essentially replace $L\to\delta k^{-1}$ in
$g$ to get the effective coupling to a radiation signal of bandwidth $\delta
k$. A more general and rigorous form of this discussion is provided in App. E.
Our simplifying assumptions here make this the “best case” scenario: we have a
narrowband, perfectly mode-matched, optimally orientated detector. We will see
that even in this case, detection of non-classicality witnesses is
extraordinarily difficult; in a more realistic case the difficulties will only
increase.

First we consider how to use this as an intensity detector. We assume that at
each time step of length $\delta t$, the electromagnetic mode $a$ is prepared
in its vacuum state $\ket{0}$, and we then measure the mode $a$ in its number
basis $\ket{n}$, for example with a photodiode. Suppose to begin with that the
gravitational radiation mode is in some pure state $\ket{\psi_{\rm rad}}$, so
that the initial total state is $\ket{0,\,\psi_{\rm rad}}$. The probability
that we see $n$ photons excited in the detector mode in the timestep $\delta
t$ is given by

$\displaystyle p(n)$ $\displaystyle=\sum_{f_{\rm rad}}|\braket{n,\,f_{\rm
rad}}{U}{0,\,\psi_{\rm rad}}|^{2}$ (13) $\displaystyle=\braket{\psi_{\rm
rad}}{K_{n}^{\dagger}K_{n}}{\psi_{\rm rad}}\\!,$

where the sum over the final state of the gravitational field $\ket{f_{\rm
rad}}$ arises because we do not measure it directly. Here

$U=\exp\left[-i\int_{0}^{\delta t}dt\,V_{I}(t)\right]$ (14)

is the interaction-picture evolution operator, which encodes the interaction
between the radiation and detector. We have used this to define the Krauss
operator

$K_{n}=\braket{n}{U}{0}\\!,$ (15)

which acts on the radiation field. These Krauss operators generalize
projective measurements, which can be recovered when the operators
$E_{n}=K_{n}^{\dagger}K_{n}$ are projectors.101010In measurement theory, the
operators $E_{n}=K_{n}^{\dagger}K_{n}$ form what is called a positive
operator-valued measure (POVM) [40]. They satisfy $\sum_{n}E_{n}=1$ and
therefore give an operator-valued measure on the space of outcomes, in this
case, values of $n$ that can be registered by the detector. If the radiation
is in a mixed state $\rho_{\rm rad}$, then Eq. (13) is generalized to
$p(n)={\rm tr}[E_{n}\rho_{\rm rad}]$. The gravitational coupling is so weak
that we are instead making a very weak, non-projective measurement. For
example, with our specific coupling Eq. (12), we have

$K_{n}^{\dagger}K_{n}=(1-g^{2}\delta
t^{2}b^{\dagger}b)\,\delta_{n,0}+g^{2}\delta t^{2}b^{\dagger}b\,\delta_{n,1},$
(16)

up to terms of ${\mathcal{O}}(g^{4}\delta t^{4})$. What this means is that
with overwhelming probability, the detector will have zero excited photons,
while the next most likely outcome, with probability $g^{2}\delta
t^{2}\braket{\psi_{\rm rad}}{b^{\dagger}b}{\psi_{\rm rad}}$, is that a single
detector photon is excited.

Essentially, this says that our detector is a graviton version of a
photodiode: it is excited with probability proportional to the incoming
intensity $I\sim\omega b^{\dagger}b$. Explicitly, in the model where $b$ is a
classical c-number, we find that the probability of a single detector
excitation is

$p(1)=g^{2}\delta
t^{2}|b|^{2}=\frac{1}{2}\frac{\rho}{\omega}\frac{B_{0}^{2}\ell^{3}}{M_{\rm
Pl}^{2}}\delta t^{2},$ (17)

using Eq. (12). Thus, the conversion probability is proportional to the number
density of incident gravitons, $\rho/\omega$. This expression is reminiscent
of the Gertsenshtein effect, and indeed we will show the connection is exact.

The same basic detector architecture can also be used to implement linear
(“amplitude,” “quadrature”) measurements on the gravitational radiation mode.
We define the electromagnetic field quadrature operators

$X=\frac{a+a^{\dagger}}{\sqrt{2}},\hskip
14.22636ptY=-i\frac{a-a^{\dagger}}{\sqrt{2}},$ (18)

which are canonically conjugate $[X,Y]=i$ and act like position and momentum
operators; these are sometimes called the amplitude and phase variables.
Suppose that instead of monitoring the electromagnetic mode $a$ with a
photocounting device, we monitor it with something like a homodyne
interferometer (see Fig. 2). This means that we projectively measure it in one
of the (continuous) quadrature bases, say the amplitude basis $\ket{X}$. Again
we assume that we do this repeatedly every time step $\delta t$.111111We are
making another important simplifying assumption here, which is that we know
the phase of the incoming signal. In a more realistic setting, say with an
astrophysical source, the unknown relative phase $\phi$ between the signal and
the detector mode would have to appear either in the coupling $g\to
ge^{i\phi}$ in Eq. (12) or in the measurement basis $X\to
X(\phi)=(ae^{i\phi}+a^{\dagger}e^{-i\phi})/\sqrt{2}$. As with our other
simplifying assumptions, this assumption will only make the task of detecting
field quantization _easier_ , and since as we will see this is already
impossibly difficult, we will continue to assume a phase-matched detector for
a simpler presentation. The probability distribution for the detector to be
measured with value $X$ is121212Note that in this equation $\ket{0}$ is still
the vacuum $n=0$ state, not the amplitude $X=0$ eigenstate.

$p(X)=\braket{\psi_{\rm rad}}{K_{X}^{\dagger}K_{X}}{\psi_{\rm rad}}\\!,\ \ \
K_{X}=\braket{X}{U}{0}.$ (19)

Again using our specific interaction Hamiltonian Eq. (12), we find to linear
order in $g\delta t$,

$K_{X}^{\dagger}K_{X}=f_{0}(X)+g\delta tf_{1}(X)[b+b^{\dagger}],$ (20)

where

$\displaystyle f_{0}(X)$ $\displaystyle=\left|\braket{X}{0}\right|^{2},$ (21)
$\displaystyle f_{1}(X)$ $\displaystyle=\braket{0}{X}\braket{X}{1}\\!,$

involve the position-space harmonic oscillator wavefunctions $\braket{X}{n}$,
and we used the fact that $f_{1}(X)$ is real. We see that the detector can in
principle be excited to any value of $X$. The probability distribution on
these detector $X$ outcomes is shifted by an amount proportional to the
_gravitational radiation field_ quadrature $X_{\rm
b}=(b+b^{\dagger})/\sqrt{2}$. For example, with an incident coherent state
$\ket{\psi_{\rm rad}}=\ket{\beta}$, the average value on the detector is
$\overline{X}=\sqrt{2}g\delta t{\rm Re}(\beta)\propto h$, where $h$ is the
usual strain observable. Thus this constitutes a gravitational amplitude
detector.

We emphasize that a gravitational radiation intensity detector is one that
measures in the number basis. In this sense one can “detect single gravitons”
simply by building such a detector. As we will explain next, traditional
quantum optics experiments demonstrating quantization of the photon have used
this kind of detector, operated in a regime where the incident light is
intense (i.e., $n\lambda_{\rm dB}^{3}\gg 1$). Therefore, to detect a
quantization measure, it is _not_ necessary to be detecting in the dilute beam
limit $n\lambda_{\rm dB}^{3}\ll 1$, contrary to the implicit suggestion in
Dyson’s original work [3]. Conversely, detecting a signal in this dilute limit
does not imply or require a quantized description of the signal field.

###  I.3 Intensity measurements

We now turn to intensity measurements on the gravitational radiation and how
they can be used to test the quantization of the gravitational field. What we
will do is calculate the most general signal possible with the classical model
in Eq. (6), then perform the analogous computation with the quantized model of
Eq. (5), and then show some examples of intensity data that can be explained
in the quantized but not the classical model. Following our discussion above,
we will model the detector data as a series of clicks: in each interval
$\delta t$, either $0$ or $1$ click occurs. Over a total integration time
$T\gg\delta t$, we can then ask for the probability of observing $N$ clicks,
$p(N)$.

In the classical case, Eqs. (4) or (6), the general state of the radiation
mode is described by a probability distribution $P_{\rm cl}(b)$ for the
Fourier amplitude. First consider an incoming state with a fixed value of $b$.
Assuming that the coupling $g$ in Eq. (12) is weak $g\delta t\ll 1$, the only
measurement outcomes in each timestep $\delta t$ with reasonable probability
are zero clicks (with $p(0)=1-g^{2}\delta t^{2}|b|^{2}$), and one click, with
$p(1)=g^{2}\delta t^{2}|b|^{2}$ from Eq. (17). The average number of clicks in
time $T$ is then

$\overline{N}=p(1)\frac{T}{\delta t}=\eta\phi T,$ (22)

where we have defined

$\phi=\frac{\rho\ell^{2}}{\omega}=\frac{|b|^{2}\ell^{2}}{4L^{3}},\ \ \
\eta=4g^{2}\frac{L^{3}}{\ell^{3}}\ell\delta t=\frac{B_{0}^{2}}{2M_{\rm
Pl}^{2}}\ell\delta t,$ (23)

which represent the incoming flux $\phi$ and detector efficiency $\eta$. (In
the quantum picture, $\phi$ corresponds to the number of gravitons incident on
the detector per unit time, whereas the efficiency determines the fraction of
those gravitons that are absorbed.)

Within a long time interval $T\gg\delta t$, we imagine collecting $N$ total
clicks. Because we are rapidly resetting the detector at each time step, the
individual detection events are independent, and so will be Poisson
distributed, and thereby fully described by the mean in Eq. (22). Thus we have

$p(N|\,b)=\frac{1}{N!}[\eta\phi T]^{N}e^{-\eta\phi T}.$ (24)

The variance in the number of clicks is then given by the usual Poisson
formula,

$\Delta N^{2}=\overline{N}.$ (25)

Notice, in particular, that the signal $b$ here is completely deterministic;
this variance represents randomness in the detector purely due to the
detector’s quantum fluctuations.

In the most general classical model, $b$ is not fixed, but instead is a random
variable drawn from a distribution $P_{\rm cl}(b)$. As a result, the flux is
stochastic, and therefore so too is the expected number of events in a small
time interval, $\eta\phi(b)\delta t$, where we have made explicit the fact
that the flux $\phi=\phi(b)$ depends on $b$ through Eq. (23). Nevertheless,
the probability for observing $N$ clicks in a time $T$ can be derived just as
in the Poisson case,

$\displaystyle p(N)=$ $\displaystyle\frac{1}{N!}\int db\,P_{\rm
cl}(b)[\eta\phi T]^{N}e^{-\eta\phi T}.$ (26)

The click variance is now given by

$\displaystyle\Delta N^{2}$ $\displaystyle=\overline{N}+\eta^{2}T^{2}\int
db\,P_{\rm cl}(b)\left[\phi-\braket{\phi}\right]^{2}$ (27)
$\displaystyle=\overline{N}+\eta^{2}T^{2}\big{\langle}\left[\phi-\braket{\phi}\right]^{2}\big{\rangle},$

where here $\braket{\cdot}$ denotes an integral over the classical
distribution, so that $\braket{\phi}=\int db\,P_{\rm cl}(b)\phi(b)$. In the
case where $b$ is fixed, $P_{\rm cl}(b)=\delta(b-b^{\prime})$, the final term
on the right hand side vanishes, reproducing Eq. (25). More generally, as a
classical probability distribution satisfies $P_{\rm cl}(b)\geq 0$, Eq. (27)
implies that $\Delta N^{2}\geq\overline{N}$: the click statistics must appear
either as Poisson or “super-Poisson.” Note, however, that as
$\overline{N}\propto\eta$, whereas the super-Poisson contribution scales as
$\eta^{2}$: the statistics revert to Poisson for $\eta\ll 1$.

Now consider the model in Eqs. (2) or (5), where the incident field is
quantized into photons or gravitons. Contrary to what we concluded for the
classical field model, there are quantum states that can produce sub-Poisson
statistics, $\Delta N^{2}<\overline{N}$. We will analyze this in direct
analogy with the classical case. First, consider an incoming radiation mode in
a definite coherent state $\ket{\beta}$. The probability of a single click in
the detector is

$p(1)=\braket{\beta}{K_{1}^{\dagger}K_{1}}{\beta}=g^{2}\delta
t^{2}|\beta|^{2}.$ (28)

Moreover, the state of the incident radiation mode after this measurement is
still $\ket{\beta}$.131313After a measurement described by a Krauss operator
$K_{a}$, the state $\ket{\psi}\to K_{a}\ket{\psi}/\|K_{a}\ket{\psi}\|$ where
$a$ is the measurement outcome, which is the generalization of the usual
projection postulate. In this particular case, since $K_{1}\propto b$, the
measurement leaves a coherent state $\ket{\beta}$ unchanged to leading order
in $g\delta t\ll 1$. This is why the coherent state basis is used in the
theory of photodetection. Thus in a time interval $T\gg\delta t$, we can have
$N$ clicks which will again be exactly Poisson distributed,

$p(N|\,\beta)=\frac{1}{N!}[\eta\phi T]^{N}e^{-\eta\phi T},$ (29)

where now $\phi=|\beta|^{2}\ell^{2}/4L^{3}$, so that $\overline{N}=\eta\phi
T$. Here, $\hat{\phi}\sim b^{\dagger}b$ is now an operator, and
$\phi=\phi(\beta)=\braket{\beta}{\hat{\phi}}{\beta}$ is its expectation value
in a specific coherent state $\ket{\beta}$. As claimed, $\phi$ is a direct
measure of the number of incident gravitons per unit time. In the general
case, the state of the mode is described by a density matrix $\rho_{\rm rad}$,
which can be described with the Glauber-Sudarshan representation of Eq. (8).
In analogy to the classical case, one then finds that the variance in the
observed detector counts is

$\Delta N^{2}=\overline{N}+\eta^{2}T^{2}\int
d\beta\,P(\beta)\left[\phi-\braket{\phi}\right]^{2}\\!,$ (30)

where now $\braket{\phi}={\rm Tr}[\rho_{\rm rad}\phi]$. See App. F for the
detailed calculation. In spite of the superficial similarity to Eq. (27),
there is a key difference: in a general radiation field state, $P(\beta)$ can
be negative for certain values of $\beta$, and the integral in Eq. (30) is not
necessarily positive. Therefore, observation of sub-Poisson detector clicks,
$\Delta N^{2}<\overline{N}$, can be explained by the quantum model of the
radiation field but not the classical version.

Detecting sub-Poisson fluctuations is a formidable task. This was only
accomplished in the 1980’s in optics [41, 42]. As discussed in Sec. I.1, there
are two reasons. The first is that we need a source state that can produce a
Glauber distribution which has $P(\beta)<0$ for at least some values of
$\beta$. A coherent state, where $P(\beta)=\delta(\beta-\beta^{\prime})$, does
not satisfy this condition, although it can be achieved for example with a
phase-squeezed state (see App. C). However, in the gravitational case, perhaps
the more important issue is that one needs a high-efficiency detector. This is
clear from Eq. (30): if $\eta\ll 1$, the actual data will be Poisson
distributed. Exactly as for the classical model in Eq. (27), any deviation
from Poisson statistics is suppressed by $\eta$. Heuristically, this is a
statement that the detected sample statistics are a very poor approximation to
the underlying population statistics in the signal, and are dominated by the
intrinsic Poisson statistics of the detector itself.

Modern photocounting devices can achieve $\eta\sim 0.9$ [1, 2]. What about a
gravitational intensity detector? With our toy model, we can make a crude
estimate of the sensitivity of a real detector like CAST. Taking $B_{0}\simeq
10~{}{\rm T}$ and $\ell\simeq\delta t=10~{}{\rm m}$, the efficiency is

$\eta=\frac{B_{0}^{2}}{2M_{\rm Pl}^{2}}\ell\delta t\simeq 8.2\times 10^{-34}.$
(31)

Thus, the data will always look Poissonian. In other words, there is no way
for CAST to make a counting statistics measurement that can differentiate
between a classical gravitational wave and a quantized beam of gravitons.
While this is very suggestive that any realistic detector cannot make such a
measurement, a natural question is whether it is possible in principle. We
analyze this in Sec. II.

Finally, it is important to note that there are observables in intensity
detectors other than sub-Poisson statistics that can distinguish between the
classical and quantized models. Examples include antibunching as observed in
$g^{2}(\tau)$ and higher-order autocorrelation measurements [34]. One can also
consider entangled states (e.g., simultaneous detection of entangled graviton
pairs). We have not systematically studied all of these, but in each example,
the basic issue of detector efficiency should be highly prohibitive
independent of the specific observable considered.

###  I.4 Amplitude measurements

Linear detectors of radiation fields also admit tests of quantization
analogous to sub-Poisson counting statistics. Consider the amplitude detector
described by Eq. (20). We read out the detector by measuring its amplitude
quadrature $X$, which corresponds to a measurement of the gravitational field
amplitude $X_{b}=(b+b^{\dagger})/\sqrt{2}$. As a first guess, we could try to
look at the variance of this observable, in direct analogy with the counting
statistics of the previous section.

One finds a result which is very similar in spirit to the number counting
statistics. The detector data will have variance given by

$\Delta X^{2}=\frac{1}{2}+g^{2}\delta t^{2}\begin{cases}\braket{\Delta
X_{b}^{2}}\\!,&{\rm classical},\\\ \big{[}\\!\braket{\Delta
X_{b}^{2}}\\!-\frac{1}{2}\big{]},&{\rm quantum}.\end{cases}$ (32)

Writing this more explicitly, in the classical radiation case, we have

$\Delta X^{2}=\frac{1}{2}+2g^{2}\delta t^{2}\int db\,P_{\rm cl}(b)\left[{\rm
Re}~{}b-\braket{{\rm Re}~{}b}\right]^{2}\\!,$ (33)

while in the model with quantized radiation, we have instead

$\Delta X^{2}=\frac{1}{2}+2g^{2}\delta t^{2}\int d\beta\,P(\beta)\left[{\rm
Re}~{}\beta-\braket{{\rm Re}~{}\beta}\right]^{2}\\!.$ (34)

A detailed calculation is provided in App. F. In both the classical and
quantum cases, the first $1/2$ term represents the vacuum fluctuations of the
detector mode itself (i.e., noise at the “standard quantum limit”).141414As is
well known, this is not a fundamental limit. The $1/2$ can in principle be
reduced arbitrarily close to $0$ by using various quantum measurement
techniques, such as the preparing the detector mode in a squeezed state
instead of the vacuum. However, this only helps in the case that the phase of
the incoming signal is known (see footnote 11). Moreover, in a realistic
device, squeezing can at most achieve a reduction of this factor $1/2$ by a
few orders of magnitude before becoming limited by optical losses. The second
terms are essentially the variance of the gravitational mode $X_{b}$, which is
transduced onto the detector mode $X$ with efficiency $g\delta t\ll 1$.

Much like the case of number counting statistics, the key difference is that
the integral over $P_{\rm cl}(b)$ is positive-definite, while the integral
over the quantum $P(\beta)$ can be negative. In particular, this means that
with particular incoming quantum states—for example, a highly squeezed
gravitational mode—it is possible that the detector’s data will actually
_decrease_ in variance compared to its native vacuum fluctuations. This
behavior is impossible with an incoming classical signal.

However, observing this effect is probably impossible in practice, because
$g\delta t$ is astronomically small in a real detector. In particular, the
maximum quantum effect, even with an infinitely squeezed state, would be to
render $\Delta X^{2}=1/2-g^{2}\delta t^{2}$. In this case, the gravitational
amplitude is fluctuating less than in its vacuum state. This is directly
analogous to the counting statistics: quantum mechanics allows for radiation
states which will produce “less noisy” data than any classical state (because
they can act to reduce the noise in the detector itself), but observing this
requires overcoming a detector efficiency $\eta\sim g\delta t\ll 1$.

In Refs. [27, 28], it was argued that observing shot noise induced from a
highly squeezed, large-amplitude gravitational state would constitute evidence
for the quantization of gravity. We agree that a highly squeezed state can
lead to an exponential enhancement of the measured strain noise (see App. C),
where $\braket{\Delta X_{b}^{2}}\gg 1$. However, observing this enhanced noise
would not provide any information about the quantization of the field, since
either the quantum or classical radiation models can produce the same data.
This does not require an exotic or highly tuned classical model; a simple
Gaussian distribution of classical plane waves with the same width
$\braket{\Delta X^{2}_{b}}$ would suffice to reproduce the same signal. The
only “uniquely quantum” signal would be if the minute fluctuations of the
noise have $\braket{\Delta X_{b}^{2}}<1/2$, i.e., if the gravitational
amplitude is fluctuating less than its vacuum state. As we have argued above,
observing this would be extremely challenging with a realistic detector.

For completeness, we also note that there are other quadrature observables
that can distinguish between the classical and quantized cases. In principle,
for example, one could fully reconstruct the underlying density matrix using a
tomographic set of measurements [43, 29], which in particular would allow one
to check if $P(\beta)<0$. However, we expect any of these observables to be
similarly limited by the detector efficiency issues discussed here.

##  II Can the quantization signatures be detected in principle?

In the previous section, we described experimental measurements which could
clearly demonstrate that the gravitational radiation field is quantized. Even
if we allowed for sources of the appropriate non-classical graviton states,
all such measurements appear prohibitively difficult with realistic
terrestrial detectors. The fundamental obstruction was encoded in Eq. (31):
signatures of non-classicality are highly suppressed by the miniscule
detection efficiency. In this section, we turn to the question of whether it
is possible even in principle, given the known laws of physics, to construct a
detector capable of observing these subtle quantization effects in the
gravitational radiation field.

To do so, let us return to the example of an amplitude detector based on
graviton to photon conversion, i.e. a futuristic version CAST. Naively, the
basic question is whether the laws of physics will allow us to increase
$B_{0}$ and/or the detector size $\ell$ such that the detector efficiency
$\eta\approx 1$. In this case, the explicit challenge is the small graviton-
photon conversion probability, $p(g\to\gamma)$, which is controlled by the
inverse Gertsenshtein effect, and well studied as a classical wave mixing
problem, see Refs. [9, 44, 45]. The probability that a graviton converts to a
photon over the length of the detector, however, is exactly the detector
efficiency in the language of Eq. (23), with $\delta t=\ell$ set by how long
the graviton takes to traverse the detector. Accordingly, we can also compute
the efficiency from the inverse Gertsenshtein effect,

$\eta=\frac{B_{0}^{2}}{2M_{\rm Pl}^{2}}\ell_{\rm
osc}^{2}\sin^{2}\left(\frac{\ell}{\ell_{\rm osc}}\right)\\!.$ (35)

Here $\ell$ denotes the length over which the oscillation can occur within a
homogeneous magnetic field $B_{0}$ orthogonal to the propagation direction,
exactly as in our discussion in Sec. I.

However, this general result includes an additional effect: self-interactions
of the electromagnetic field, as encoded in the Euler-Heisenberg Lagrangian.
These imply that increasing $B_{0}$ will not arbitrarily enhance the detection
efficiency. Instead, a larger $B_{0}$ eventually reduces the coherence of the
photon and graviton, thereby reducing the mixing. This is captured by the
oscillation length $\ell_{\rm osc}$ that is determined from $4/\ell_{\rm
osc}^{2}=m_{\rm EH}^{4}/\omega^{2}+2B_{0}^{2}/M_{\rm Pl}^{2}$, where for a
single polarization $m_{\rm EH}^{2}=7\omega^{2}(e^{2}/180\pi^{2})(B_{0}/B_{\rm
crit})^{2}$ and $B_{\rm crit}=m_{e}^{2}/e$ with $m_{e}$ and $e$ the electron
mass and charge, respectively. The result as stated neglects any possible
plasma mass for the photon, and further assumes $\omega\gg m_{\rm EH}$, which
is justified for all parameters we will consider.

For CAST, $B_{0}=10~{}{\rm T}$, $f=10^{18}~{}{\rm Hz}$, and $\ell=10~{}{\rm
m}$, so that in natural Heaviside-Lorentz units

$\frac{2B_{0}^{2}}{M_{\rm Pl}^{2}}\simeq\left(5.6~{}{\rm
pc}\right)^{-2}\\!,\hskip 14.22636pt\frac{m_{\rm
EH}^{4}}{\omega^{2}}\simeq\left(0.17~{}{\rm AU}\right)^{-2}\\!.$ (36)

Accordingly $\ell_{\rm osc}\simeq 0.35~{}{\rm AU}\gg\ell$, implying that the
Euler-Heisenberg corrections are irrelevant. Indeed, when $\ell\ll\ell_{\rm
osc}$ the efficiency in Eq. (35) reduces to Eq. (31), as it must,
demonstrating the equivalence of the classical and quantum derivation of the
Gertsenshtein process. If we seek to maximize $\eta$, $\ell_{\rm osc}$ must be
accounted for. Doing so for $\ell=1~{}{\rm AU}$, the maximum efficiency occurs
for $B_{0}\simeq 6~{}{\rm T}$, but remains small, $\eta\simeq 5\times
10^{-14}$. Over cosmological distances sizeable values can be achieved: for
$\ell=1~{}{\rm Mpc}$ ($1~{}{\rm Gpc}$) and $B_{0}=0.1~{}{\rm G}$ ($1~{}{\rm
mG}$) we can achieve $\eta\simeq 1\%$ ($70\%$). While such large devices are
certainly unrealistic, we do not see a physical law that would prohibit such
an instrument existing in the Universe, thereby allowing
$\eta\sim{\mathcal{O}}(1)$. The observation that CAST has already reached the
dilute graviton regime, while even a Gpc version barely meets the requirements
to test quantization underlines how much more difficult the latter task is.

##  III Outlook and relation to other tests

Gravity can be consistently quantized at low energy densities as an effective
quantum field theory [46, 47]. An important question is to determine what
kinds of measurements could be used to verify that this is the way gravity
really operates in nature. We have analyzed Dyson’s suggestion that one could
try to make a detector that is sensitive to single gravitons. The conclusion
is straightforward: it is possible to make such a detector, but simply seeing
the detector click due to incoming gravitational radiation could be explained
either with quantized gravitons or simple classical gravitational waves. Vice
versa, from the point of view of high-frequency gravitational wave searches
venturing into the parameter space shown in Fig. 1, the dilute graviton regime
is a red herring: crossing this line in the parameter space does not come with
any particular change in the difficulty of performing a gravitational wave
measurement.

In order to make a measurement that cleanly distinguishes between the quantum
and classical cases, two very substantial hurdles have to be overcome: a
source of non-classical (e.g., squeezed or entangled) gravitational radiation
has to be produced, and then detected with a gravitational detector with high
efficiency. We do not see a realistic path to satisfying either of these
conditions, although we also have no argument that nature is fundamentally
incapable of satisfying them.

There are two known alternatives to testing low-energy quantum gravity. One is
to look at fluctuations in the cosmic microwave background, assuming these are
produced by a phase of cosmic inflation. As the arguments in this paper should
make clear, observation of a Gaussian power spectrum (in either the scalar or
tensor fluctuations) can easily be explained by either classical or quantized
gravitational radiation. Similarly, the measurement of CMB tensor modes with
superhorizon correlations (as proposed, e.g., in Ref. [48]), without
independent evidence pointing to vacuum fluctuations of inflation, can also be
explained by either hypothesis. It may however be possible to produce good
non-classicality witnesses with higher-point observables [49, 50].

The other is to use tabletop experiments to observe quantum effects generated
by the two-body Newton interaction
$\hat{V}=G_{N}m_{1}m_{2}/|\hat{{\mathbf{x}}}_{1}-\hat{{\mathbf{x}}}_{2}|$ (for
a review, see Ref. [51]). It is known that a coherent, entangling Newton
interaction is only consistent with unitarity and Lorentz invariance if
gravitons also exist in the spectrum of scattering states [52, 53]. Thus
experimental verification that the low-energy interaction is coherent and
unitary would provide an indirect proof that the radiation field is also
quantized. We view the clear difficulties with measurements directly on the
gravitational radiation presented in this paper as substantial motivation to
pursue both of these programs.

## Acknowledgements

We thank Markus Aspelmeyer, Diego Blas, Camilo Garcia-Cely, Daniel Green,
Thiago Guerriero, Giacomo Marocco, Maulik Parikh, Grant Remmen, Jess Riedel,
Mariana Safronova, John Selby, Yotam Soreq, Raman Sundrum, Jacob Taylor, Jure
Zupan, and Kathryn Zurek for discussions. DC is supported in part by the U.S.
Department of Energy QuantISED program, grant KA2401032, and by the Heising-
Simons foundation, grant 2023-4467.

## References

  * [1] R. H. Hadfield, “Single-photon detectors for optical quantum information applications,” Nature Photonics 3 no. 12, (Dec., 2009) 696–705.
  * [2] J.-C. Besse, S. Gasparinetti, M. C. Collodo, T. Walter, P. Kurpiers, M. Pechal, C. Eichler, and A. Wallraff, “Single-shot quantum nondemolition detection of individual itinerant microwave photons,” Phys. Rev. X 8 (Apr, 2018) 021003.
  * [3] F. Dyson, “Is a graviton detectable?,” Int. J. Mod. Phys. A 28 (2013) 1330041.
  * [4] T. Rothman and S. Boughn, “Can gravitons be detected?,” Found. Phys. 36 (2006) 1801–1825, arXiv:gr-qc/0601043.
  * [5] S. Boughn and T. Rothman, “Aspects of graviton detection: Graviton emission and absorption by atomic hydrogen,” Class. Quant. Grav. 23 (2006) 5839–5852, arXiv:gr-qc/0605052.
  * [6] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford Master Series in Physics. Oxford University Press, 2007.
  * [7] CAST Collaboration, V. Anastassopoulos et al., “New CAST Limit on the Axion-Photon Interaction,” Nature Phys. 13 (2017) 584–590, arXiv:1705.02290 [hep-ex].
  * [8] M. Gertsenshtein, “Wave resonance of light and gravitational waves.” Sov.Phys.JETP 14 (1962) 84.
  * [9] G. Raffelt and L. Stodolsky, “Mixing of the Photon with Low Mass Particles,” Phys. Rev. D 37 (1988) 1237.
  * [10] A. Ejlli, D. Ejlli, A. M. Cruise, G. Pisano, and H. Grote, “Upper limits on the amplitude of ultra-high-frequency gravitational waves from graviton to photon conversion,” Eur. Phys. J. C 79 no. 12, (2019) 1032, arXiv:1908.00232 [gr-qc].
  * [11] A. Berlin, D. Blas, R. Tito D’Agnolo, S. A. R. Ellis, R. Harnik, Y. Kahn, and J. Schütte-Engel, “Detecting high-frequency gravitational waves with microwave cavities,” Phys. Rev. D 105 no. 11, (2022) 116011, arXiv:2112.11465 [hep-ph].
  * [12] V. Domcke, C. Garcia-Cely, and N. L. Rodd, “Novel Search for High-Frequency Gravitational Waves with Low-Mass Axion Haloscopes,” Phys. Rev. Lett. 129 no. 4, (2022) 041101, arXiv:2202.00695 [hep-ph].
  * [13] L. Hui, “Wave Dark Matter,” Ann. Rev. Astron. Astrophys. 59 (2021) 247–289, arXiv:2101.11735 [astro-ph.CO].
  * [14] E. Armengaud et al., “Conceptual Design of the International Axion Observatory (IAXO),” JINST 9 (2014) T05002, arXiv:1401.3233 [physics.ins-det].
  * [15] N. Aggarwal et al., “Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies,” Living Rev. Rel. 24 no. 1, (2021) 4, arXiv:2011.12414 [gr-qc].
  * [16] A. Ringwald, J. Schütte-Engel, and C. Tamarit, “Gravitational Waves as a Big Bang Thermometer,” JCAP 03 (2021) 054, arXiv:2011.04731 [hep-ph].
  * [17] M. Goryachev, W. M. Campbell, I. S. Heng, S. Galliou, E. N. Ivanov, and M. E. Tobar, “Rare Events Detected with a Bulk Acoustic Wave High Frequency Gravitational Wave Antenna,” Phys. Rev. Lett. 127 no. 7, (2021) 071102, arXiv:2102.05859 [gr-qc].
  * [18] N. Herman, A. Füzfa, L. Lehoucq, and S. Clesse, “Detecting planetary-mass primordial black holes with resonant electromagnetic gravitational-wave detectors,” Phys. Rev. D 104 no. 2, (2021) 023524, arXiv:2012.12189 [gr-qc].
  * [19] A. Berlin, D. Blas, R. Tito D’Agnolo, S. A. R. Ellis, R. Harnik, Y. Kahn, J. Schütte-Engel, and M. Wentzel, “MAGO 2.0: Electromagnetic Cavities as Mechanical Bars for Gravitational Waves,” arXiv:2303.01518 [hep-ph].
  * [20] A. Ito and J. Soda, “Exploring High Frequency Gravitational Waves with Magnons,” arXiv:2212.04094 [gr-qc].
  * [21] K. Schmieden and M. Schott, “Searching for Gravitational Waves with CMS,” arXiv:2209.12024 [gr-qc].
  * [22] M. E. Tobar, C. A. Thomson, W. M. Campbell, A. Quiskamp, J. F. Bourhill, B. T. McAllister, E. N. Ivanov, and M. Goryachev, “Comparing Instrument Spectral Sensitivity of Dissimilar Electromagnetic Haloscopes to Axion Dark Matter and High Frequency Gravitational Waves,” Symmetry 14 no. 10, (2022) 2165, arXiv:2209.03004 [physics.ins-det].
  * [23] V. Domcke, C. Garcia-Cely, S. M. Lee, and N. L. Rodd, “Symmetries and Selection Rules: Optimising Axion Haloscopes for Gravitational Wave Searches,” arXiv:2306.03125 [hep-ph].
  * [24] T. Bringmann, V. Domcke, E. Fuchs, and J. Kopp, “High-Frequency Gravitational Wave Detection via Optical Frequency Modulation,” arXiv:2304.10579 [hep-ph].
  * [25] N. Bao, A. Chatwin-Davies, J. Pollack, and G. N. Remmen, “Cosmological Decoherence from Thermal Gravitons,” JHEP 08 (2020) 065, arXiv:1911.10207 [hep-th].
  * [26] T. Guerreiro, “Quantum Effects in Gravity Waves,” Class. Quant. Grav. 37 no. 15, (2020) 155001, arXiv:1911.11593 [quant-ph].
  * [27] M. Parikh, F. Wilczek, and G. Zahariade, “Quantum Mechanics of Gravitational Waves,” Phys. Rev. Lett. 127 no. 8, (2021) 081602, arXiv:2010.08205 [hep-th].
  * [28] M. Parikh, F. Wilczek, and G. Zahariade, “Signatures of the quantization of gravity at gravitational wave detectors,” Phys. Rev. D 104 no. 4, (2021) 046021, arXiv:2010.08208 [hep-th].
  * [29] T. Guerreiro, F. Coradeschi, A. M. Frassino, J. R. West, and E. Schioppa, Junior., “Quantum signatures in nonlinear gravitational waves,” Quantum 6 (2022) 879, arXiv:2111.01779 [gr-qc].
  * [30] T. Guerreiro, “Nonlinearities in Black Hole Ringdowns and the Quantization of Gravity,” arXiv:2306.09974 [gr-qc].
  * [31] L. Mandel and E. Wolf, Optical Coherence and Quantum Optics.  Cambridge University Press, 1995.
  * [32] L. Davidovich, “Sub-poissonian processes in quantum optics,” Rev. Mod. Phys. 68 (Jan, 1996) 127–173.
  * [33] E. C. G. Sudarshan, “Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams,” Phys. Rev. Lett. 10 (1963) 277–279.
  * [34] R. J. Glauber, “Coherent and incoherent states of the radiation field,” Phys. Rev. 131 (1963) 2766–2788.
  * [35] B. R. Mollow, “Pure-state analysis of resonant light scattering: Radiative damping, saturation, and multiphoton effects,” Phys. Rev. A 12 (Nov, 1975) 1919–1943.
  * [36] 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 (Feb, 2002) 032323.
  * [37] N. Crisosto, P. Sikivie, N. S. Sullivan, D. B. Tanner, J. Yang, and G. Rybka, “ADMX SLIC: Results from a Superconducting $LC$ Circuit Investigating Cold Axions,” Phys. Rev. Lett. 124 no. 24, (2020) 241101, arXiv:1911.05772 [astro-ph.CO].
  * [38] HAYSTAC Collaboration, L. Zhong et al., “Results from phase 1 of the HAYSTAC microwave cavity axion experiment,” Phys. Rev. D 97 no. 9, (2018) 092001, arXiv:1803.03690 [hep-ex].
  * [39] A. Melissinos and A. Das, “The response of laser interferometers to a gravitational wave,” American Journal of Physics 78 no. 11, (2010) 1160–1164.
  * [40] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information.  Cambridge university press, 2010.
  * [41] R. Short and L. Mandel, “Observation of sub-poissonian photon statistics,” Phys. Rev. Lett. 51 (Aug, 1983) 384–387.
  * [42] G. Rempe, F. Schmidt-Kaler, and H. Walther, “Observation of sub-poissonian photon statistics in a micromaser,” Phys. Rev. Lett. 64 (Jun, 1990) 2783–2786.
  * [43] A. I. Lvovsky and M. G. Raymer, “Continuous-variable optical quantum-state tomography,” Reviews of modern physics 81 no. 1, (2009) 299\.
  * [44] D. Ejlli and V. R. Thandlam, “Graviton-photon mixing,” Phys. Rev. D 99 no. 4, (2019) 044022, arXiv:1807.00171 [gr-qc].
  * [45] V. Domcke and C. Garcia-Cely, “Potential of radio telescopes as high-frequency gravitational wave detectors,” Phys. Rev. Lett. 126 no. 2, (2021) 021104, arXiv:2006.01161 [astro-ph.CO].
  * [46] J. F. Donoghue, “General relativity as an effective field theory: The leading quantum corrections,” Phys. Rev. D 50 (1994) 3874–3888, arXiv:gr-qc/9405057.
  * [47] J. F. Donoghue, “Quantum General Relativity and Effective Field Theory,” arXiv:2211.09902 [hep-th].
  * [48] L. M. Krauss and F. Wilczek, “Using Cosmology to Establish the Quantization of Gravity,” Phys. Rev. D 89 no. 4, (2014) 047501, arXiv:1309.5343 [hep-th].
  * [49] J. Maldacena, “A model with cosmological Bell inequalities,” Fortsch. Phys. 64 (2016) 10–23, arXiv:1508.01082 [hep-th].
  * [50] D. Green and R. A. Porto, “Signals of a Quantum Universe,” Phys. Rev. Lett. 124 no. 25, (2020) 251302, arXiv:2001.09149 [hep-th].
  * [51] D. Carney, P. C. E. Stamp, and J. M. Taylor, “Tabletop experiments for quantum gravity: a user’s manual,” Class. Quant. Grav. 36 no. 3, (2019) 034001, arXiv:1807.11494 [quant-ph].
  * [52] A. Belenchia, R. M. Wald, F. Giacomini, E. Castro-Ruiz, v. Brukner, and M. Aspelmeyer, “Quantum Superposition of Massive Objects and the Quantization of Gravity,” Phys. Rev. D 98 no. 12, (2018) 126009, arXiv:1807.07015 [quant-ph].
  * [53] D. Carney, “Newton, entanglement, and the graviton,” Phys. Rev. D 105 no. 2, (2022) 024029, arXiv:2108.06320 [quant-ph].
  * [54] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [55] O. Pisanti, G. Mangano, G. Miele, and P. Mazzella, “Primordial Deuterium after LUNA: concordances and error budget,” JCAP 04 (2021) 020, arXiv:2011.11537 [astro-ph.CO].
  * [56] T.-H. Yeh, K. A. Olive, and B. D. Fields, “The impact of new $d(p,\gamma)$3 rates on Big Bang Nucleosynthesis,” JCAP 03 (2021) 046, arXiv:2011.13874 [astro-ph.CO].
  * [57] B. Carr and F. Kuhnel, “Primordial black holes as dark matter candidates,” SciPost Phys. Lect. Notes 48 (2022) 1, arXiv:2110.02821 [astro-ph.CO].
  * [58] G. Franciolini, A. Maharana, and F. Muia, “Hunt for light primordial black hole dark matter with ultrahigh-frequency gravitational waves,” Phys. Rev. D 106 no. 10, (2022) 103520, arXiv:2205.02153 [astro-ph.CO].
  * [59] C. L. Mehta, “Diagonal coherent-state representation of quantum operators,” Phys. Rev. Lett. 18 (May, 1967) 752–754.
  * [60] W. Rudin, Functional Analysis.  McGraw-Hill, New York, 1980.

##  Appendix A CAST as a single graviton detector

Figure 1 is suggestive that CAST has the required sensitivity to detect
individual gravitons. Here we demonstrate that this is indeed the case. We
will also see that it is extremely difficult to conceive of a suitable source
close to the current sensitivity of CAST. Of course, we emphasize, as
demonstrated in the main body, even a detection would a priori not prove the
source was quantized.

To begin with, the CAST analysis in Ref. [10] reports a characteristic strain
sensitivity of $h_{c}\simeq 5\times 10^{-28}$ to a stochastic gravitational
wave background at $10^{18}\,{\rm Hz}$. This can be translated to an energy
density, $\rho=\rho_{c}\int d\ln\omega\,\Omega_{\rm gw}$, with $\Omega_{\rm
gw}=\omega^{2}h_{c}^{2}M_{\rm Pl}^{2}/2\rho_{c}$ and $\rho_{c}$ the critical
energy density. The CAST sensitivity thus corresponds to the following energy
and number density,

$\displaystyle\rho$ $\displaystyle\simeq 2\times 10^{18}~{}{\rm keV/cm}^{3},$
(37) $\displaystyle n$ $\displaystyle\simeq 4\times 10^{17}~{}{\rm
gravitons/cm}^{3}.$

Although these densities are large, the de Broglie volume is not, and we have
$n\lambda_{\rm dB}^{3}\sim 10^{-5}\ll 1$. Hence, a (weak) signal in CAST due
to GWs would imply GW detection in the single graviton limit.

However, a stochastic gravitational wave background at this amplitude is
firmly excluded by the constraints on extra radiation from BBN and CMB [54,
55, 56]. We thus proceed to estimate the sensitivity of CAST to transient
signals, which are not subject to these bounds. Let us consider a transient
signal with frequency $f=10^{18}\,{\rm Hz}$, amplitude $h$, and duration
$T_{\rm gw}$. Then, comparing the energy density of a plane GW with the CAST
limit yields a sensitivity of

$h\simeq 7\times 10^{-28}(T_{m}/T_{\rm gw})^{1/4},$ (38)

where $T_{m}\simeq 1$ year denotes the measurement time of CAST and the factor
$(T_{\rm gw}/T_{m})^{1/4}\leq 1$ accounts for finite duration of the GW signal
assuming a statistics limited experiment. As we are now considering a
transient rather than stochastic signal we have swapped to $h=\sqrt{2}h_{c}$,
which matches between the conventions of the present work and Ref. [10].
Consequently, the CAST sensitivity implies a detection threshold of
$n\lambda_{\rm dB}^{3}\simeq 10^{-5}(T_{m}/T_{\rm gw})^{1/2}$, which means
that CAST can reach the single graviton threshold for signals lasting $T_{\rm
gw}\sim 4\,{\rm ms}$ at $h\simeq 2\times 10^{-25}$.

A source with these properties is, however, difficult to conceive. The
prototypical example of an ultra high-frequency GW source is the inspiral of
light primordial black holes (PBH). Matching the CAST frequency and the
required $T_{\rm gw}$ fixes $m_{\scriptscriptstyle\textrm{PBH}}\simeq 4\times
10^{-24}~{}M_{\odot}$ (for equal masses in the binary). At such small masses,
the PBHs would dissipate due to Hawking radiation in less than a day (and for
instance are strongly excluded as being the dark matter of the Universe [57]).
Even putting this aside, in order to reach the required amplitude of $h\sim
10^{-25}$ requires a distance between the binary and the detector of
${\mathcal{O}}(\textrm{cm})$, which clearly rules this out as a viable source.
There are other possible sources – for a review see Refs. [15, 58] – although
these appear similarly challenging at best.

In summary, while interesting that CAST does reach the single graviton regime
and, in principle, one could imagine a GW with the required properties,
conceiving of a viable source for such a GW is challenging. Moreover, as
discussed in Sec. II a measurement of sub-poissonian statistics, seems
completely out of reach with any realistic version of this technology: it
would require an improved detection efficiency – with the associated
challenges discussed in Sec. II – and the source itself having a non-classical
distribution. In this sense, our conclusions are very similar to those of
Dyson [3], though we stress that the difficulty lies in finding a viable GW
source and in measuring sub-poissonian statistics, not in the single graviton
detection itself.

##  Appendix B Semi-classical photoelectric effect

Here we briefly expand on the justification that the semi-classical
Hamiltonian in Eq. (4) reproduces the famous observations of the photoelectric
effect encoded in Eq. (3). This is a well known result, and our discussion
closely follows Ref. [31]. As mentioned in the main text, even when the
electromagnetic field is classical, we continue to infer its presence with
quantized electrons, and therefore our measurements remain discrete clicks.
Further, the appearance of the intensity will remain, as the transition is
dictated by quantum mechanics, it must depend on $|{\mathbf{A}}(t)|^{2}\sim
I(t)$. More interesting is to consider the behavior when $\omega\sim\Delta$.
Working in the interaction picture, the transition probability in a small time
interval $\delta t$ is given by

$\displaystyle{\rm TP}=$ $\displaystyle\left|\int_{0}^{\delta
t}dt\,\bra{E}\hat{H}_{{\rm SC},I}(t)\ket{-\Delta}\right|^{2}$ (39)
$\displaystyle=$ $\displaystyle\left|\frac{e}{m}\int_{0}^{\delta
t}dt\,\bra{E}\hat{\bf p}\ket{-\Delta}\cdot{\bf
A}(t)\,e^{i(E+\Delta)t}\right|^{2}\\!.$

If we treat ${\mathbf{A}}$ as describing a nearly monochromatic plane wave
with frequency $\omega$, the relevant contribution will be proportional
$e^{-i\omega t}$, so that the transition probability is controlled by

$\displaystyle\left|\int_{0}^{\delta
t}dt\,e^{i(E+\Delta-\omega)t}\right|^{2}=\left[\frac{\sin[\tfrac{1}{2}(E+\Delta-\omega)\delta
t]}{\tfrac{1}{2}(E+\Delta-\omega)}\right]^{2}\\!.$ (40)

Even for $\omega<\Delta$, this quantity is non-zero. Nevertheless, for
realistic values of $\Delta$ and $\omega$, it is highly suppressed. For
$\Delta\sim\omega\sim{\rm eV}\sim 10^{15}~{}{\rm Hz}$, we expect
$(\omega-\Delta)\delta t\gg 1$, and the final line of Eq. (40) is very well
approximated by $\delta(E+\Delta-\omega)$. To determine the total probability
for excitation into the conduction band, we integrate the transition
probability over all $E$, finding no support when $\omega<\Delta$, consistent
with Eq. (3). The claim that the semi-classical model is fully consistent with
Eq. (3) therefore follows.

##  Appendix C An example non-classical $P(\beta)$ function

In the main text, we showed how a quantum state whose Glauber representation
$P(\beta)$ is negative for some values of $\beta$ can produce output in a
detector which cannot be reproduced by any classical radiation signal, even
one with classical randomness. In this appendix, we explain how $P(\beta)$
works for specific quantum state of interest: the squeezed coherent states. We
first define these states, calculate their number and amplitude variances
$\braket{\Delta n^{2}}$ and $\braket{\Delta X^{2}}$, then show how these
observables can be used to deduce the negativity of $P(\beta)$.

The squeezed coherent states are Gaussian states of a harmonic oscillator
whose variance in one quadrature is below the vacuum value, say $\Delta
X<1/2$, which by Heisenberg uncertainty requires that $\Delta Y>1/2$. They can
be efficiently defined in terms of the squeezing and displacement unitary
operators (see Sec. 21 of Ref. [31]):

$\displaystyle S(z)$
$\displaystyle=\exp\left\\{\frac{1}{2}\left[z^{*}b^{2}-zb^{{\dagger}2}\right]\right\\}\\!,$
(41) $\displaystyle D(\beta)$ $\displaystyle=\exp\left\\{\beta
b^{\dagger}-\beta^{*}b\right\\}\\!.$

Both $\beta$ and $z$ are arbitrary complex numbers. The displacement operator
rotates a given coherent state $\ket{\alpha}$ into $\ket{\beta+\alpha}$. The
squeezing operator takes a coherent state into a complicated superposition of
coherent states which we will describe shortly. To calculate statistics in
these states, it is most useful to use their transformation properties on the
creation and annihilation operators:

$\displaystyle S^{\dagger}bS=\mu b-\nu b^{\dagger},\ \
S^{\dagger}b^{\dagger}S=\mu b^{\dagger}-\nu^{*}b,$ (42) $\displaystyle
D^{\dagger}bD=b+\beta,\ \ D^{\dagger}b^{\dagger}D=b^{\dagger}+\beta^{*},$

where the complex squeezing number $z$ is parametrized $z=re^{i\theta}$, or
alternatively

$\mu=\cosh r,\ \ \nu=e^{i\theta}\sinh r.$ (43)

Often $r$ is called the squeezing amplitude and $\theta$ is the squeezing
angle. In terms of these operators, the squeezed vacuum state is
$\ket{z}:=S(z)\ket{0}$. More generally one can consider a squeezed coherent
state,151515Note that $[S,D]\neq 0$. Our definition here is what has
historically been called a “two-photon squeezed state,” while the alternative
$DS\ket{0}$ states were called “ideal squeezed states.” This choice is
somewhat arbitrary, as one can always express a two-photon squeezed state as
an ideal one, or vice versa. The terminology “two-photon” does not mean that
the states have two photons (indeed they do not have a definite number), but
rather reflects the fact that to generate $S(z)$ from a Hamiltonian
$S=e^{-iHt}$ requires $H$ to be bilinear in the photon operators. For further
discussion see Ref. [31].

$\ket{z,\beta}:=S(z)D(\beta)\ket{0}\\!.$ (44)

Generally speaking, to create a squeezed state from the vacuum or coherent
state, one requires a non-linear interaction to enact the squeeze operator.
These states have been produced in optical fields since the 1980’s, using, for
example, non-linear crystals. In the gravitational case, there is no
definitively known source of squeezed states, although inflation and the non-
linearities of gravitational mergers have been suggested as possible sources
[26, 27, 28, 29, 30].

Let us now study the nature of the “non-classicality” of the squeezed states.
The most efficient way to understand this is by computing the observables
studied in the main text, such as the signal variance $\braket{\Delta X^{2}}$,
and then comparing the specific results to the general results we showed must
hold in the classical case. Recall the definition of the amplitude quadrature
$X=(b+b^{\dagger})/\sqrt{2}$. In a general squeezed coherent state, simple
algebra using Eq. (LABEL:eq:transforms) gives the average

$\braket{X}_{z,\beta}=\frac{1}{\sqrt{2}}\left[(\mu-\nu^{*})\beta+(\mu-\nu)\beta^{*}\right]\\!,$
(45)

with $\braket{X}_{z,\beta}=\braket{z,\beta}{X}{z,\beta}$. The expectation of
$X^{2}$ is

$\braket{X^{2}}_{z,\beta}=\frac{1}{2}|\mu-\nu|^{2}\braket{bb^{\dagger}}_{0}+\braket{X}_{z,\beta}^{2}\\!,$
(46)

where we left $\braket{bb^{\dagger}}_{0}=\braket{0}{bb^{\dagger}}{0}=1$
explicit to highlight the “vacuum fluctuations.” We then obtain the variance

$\braket{\Delta X^{2}}_{z,\beta}=\frac{1}{2}|\mu-\nu|^{2}.$ (47)

For example, the true vacuum is $\mu=1,\nu=0$ (i.e. $r=0$) so that
$\braket{\Delta X^{2}}_{0}=1/2$, the usual vacuum fluctuations. What squeezing
does is to modify the variance. In particular, we can choose values of the
squeezing parameters such that $\braket{\Delta X^{2}}<1/2$ is below the vacuum
value; for example a large amount of squeezing $r\to\infty$ and $\theta=0$
sends $\braket{\Delta X^{2}}\to 0$. This does not violate Heisenberg
uncertainty: one can perform the analogous calculation in the conjugate
variable $Y=-i(b-b^{\dagger})/\sqrt{2}$ and find that for the same state,
$\braket{\Delta Y^{2}}\to\infty$. In this manner, squeezing distributes all
the uncertainty into a specific quadrature, whereas the standard vacuum
distributes it evenly between variables.

The result in Eq. (47) is enough to show that the Glauber function $P(\beta)$
for such states is negative in some region. Compare this result to the general
result in any Glauber state, our expression $\braket{\Delta X_{b}^{2}}_{\rm
qu}$ in Eq. (34). In a squeezed state where $\braket{\Delta X^{2}}<1/2$, Eq.
(34) immediately implies that the integral over $\beta$ must be negative, from
which it follows that $P(\beta)<0$ for some values of $\beta$. It is worth
emphasizing that the non-classicality of the state is independent of the
coherence $\beta$; it depends only on the squeezing parameters.

In principle, one can look to show $P(\beta)<0$ directly by computing the
Glauber function. This turns out to be difficult. For instance, we can compute
the Glauber function using a formula due to Mehta [59], which gives $P$ in
terms of the density matrix:

$P(\beta)=\frac{e^{|\beta|^{2}}}{\pi^{2}}\int
d\alpha\braket{-\alpha}{\rho}{\alpha}e^{|\alpha|^{2}}e^{\beta\alpha^{*}-\beta^{*}\alpha}.$
(48)

For example, with a coherent state $\rho=\ket{\beta_{0}}\bra{\beta_{0}}$, this
formula produces the expected $P(\beta)=\delta(\beta-\beta_{0})$. However,
inserting a squeezed state, one finds a highly singular distribution in which
$P(\beta)$ is expressed as an infinite series of derivatives acting on Dirac
delta functions $\sim(\partial/\partial\beta)^{n}\delta(\beta)$. A similar
results holds for a Fock state $\ket{n}$ with definite graviton number. This
highlights the fact that in general $P(\beta)$ is a distribution in the sense
of a Dirac delta function, not an actual function. A general theorem due to
Schwartz, however, shows that any distribution of order greater than zero
(which roughly means that it has points which behave like derivatives of the
Dirac delta) must always be negative in some region. See, for instance,
Chapter 6 of Ref. [60].

Finally, we mention that one can also produce sub-Poisson counting statistics
with squeezed signals. Let $n=b^{\dagger}b$ be the number operator of the
signal, i.e., the graviton mode. We use $n$ to distinguish this from $N$, the
number of clicks observed in the detector. Using the same algebraic tools as
above, one finds the average

$\displaystyle\braket{n}_{z,\beta}$
$\displaystyle=|\nu|^{2}+|\mu\beta^{*}-\nu^{*}\beta|^{2}=\sinh^{2}r$ (49)
$\displaystyle+|\beta|^{2}\big{[}\cosh^{2}r+\sinh^{2}r-\sinh
2r\cos(2\phi-\theta)\big{]},$

where we invoked Eq. (43) and wrote the coherence parameter as
$\beta=|\beta|e^{i\phi}$. The variance takes more work, but one obtains

$\displaystyle\braket{\Delta
n^{2}}_{z,\beta}=\braket{n}_{z,\beta}+\sinh^{2}r\cosh 2r$ (50) $\displaystyle\
\ \ +2|\beta|^{2}\sinh r\big{[}\sinh 3r-\cosh 3r\cos(\theta-2\phi)\big{]}.$

Taking $r=0$, we recover the Poissonian result, $\braket{\Delta
n^{2}}=\braket{n}$, as expected for a coherent state. The statistics of the
squeezed state are more general. To exhibit sub-Poisson variance, the final
term in Eq. (LABEL:eq:dn2-app1) must dominate. Therefore, for simplicity, let
us take $\theta=2\phi$, and consider a large-amplitude state $|\beta|\gg 1$.
Then the above results combine as

$\frac{\braket{\Delta n^{2}}-\braket{n}}{\braket{n}}=-1+e^{-2r}.$ (51)

By taking $r\gg 1$ (but $\ll\tfrac{1}{3}\ln|\beta|$) this can be arbitrarily
close to $-1$, which would mean a state with _no_ fluctuations in the number
flux $\braket{\Delta n^{2}}\to 0$, which is the most sub-Poisson distribution
possible.

##  Appendix D Beamsplitters

A general beamsplitter is defined by a scattering matrix

$\begin{pmatrix}b_{1}^{\prime}\\\
b_{2}^{\prime}\end{pmatrix}=\begin{pmatrix}t&r\\\
r^{*}&t\end{pmatrix}\begin{pmatrix}b_{1}\\\ b_{2}\end{pmatrix}\\!,$ (52)

with $b_{i}$ the input modes and $b_{i}^{\prime}$ the output modes. Let $U$ be
the unitary that implements $b_{1}\to
b_{1}^{\prime}=Ub_{1}U^{\dagger}=tb_{1}+rb_{2}$, etc. The action on a single-
excitation state is:

$\displaystyle U\ket{10}$
$\displaystyle=Ub_{1}^{\dagger}U^{\dagger}\ket{00}=(t^{*}b_{1}^{\dagger}+r^{*}b_{2}^{\dagger})\ket{00}$
(53) $\displaystyle=t^{*}\ket{10}+r^{*}\ket{01}\\!,$

where we used $U^{\dagger}\ket{00}=\ket{00}$.161616This in turn follows
because
$b_{1}U^{\dagger}\ket{00}=U^{\dagger}(Ub_{1}U^{\dagger})\ket{00}=U^{\dagger}(tb_{1}+rb_{2})\ket{00}=0$,
and similarly for $b_{2}$. The action on a pair of coherent states can be
computed similarly. The coherent states are generated by the displacement
operator of Eq. (41)

$\displaystyle\ket{\beta_{1}\beta_{2}}$
$\displaystyle=D(\beta_{1},\beta_{2})\ket{00},$ (54) $\displaystyle
D(\beta_{1},\beta_{2})$
$\displaystyle=\exp\left\\{(\beta_{1}b_{1}^{\dagger}\\!-\\!\beta_{1}^{*}b_{1})+(\beta_{2}b_{2}^{\dagger}\\!-\\!\beta_{2}^{*}b_{2})\right\\}\\!.$

Writing $D$ out as a Taylor series and inserting factors of $U^{\dagger}U=1$,
it is clear that this transforms as

$UD(\beta_{1},\beta_{2})U^{\dagger}=D(\beta_{1}^{\prime},\beta_{2}^{\prime}),$
(55)

where

$\beta_{1}^{\prime}=t^{*}\beta_{1}+r\beta_{2},\ \
\beta_{2}^{\prime}=t^{*}\beta_{2}+r^{*}\beta_{1}.$ (56)

Again using $U^{\dagger}\ket{00}=\ket{00}$, we then see that an incoming pair
of coherent states is transformed by the beamsplitter into a new pair of
coherent states,

$U\ket{\beta_{1}\beta_{2}}=UD(\beta_{1},\beta_{2})U^{\dagger}\ket{00}=\ket{\beta_{1}^{\prime},\beta_{2}^{\prime}}\\!.$
(57)

Choosing $t=r=1/\sqrt{2}$ justifies Eqs. (9) and (10) from the main text.

##  Appendix E Detector toy model

In this appendix, we provide a more detailed discussion of the detector
Hamiltonian of Eq. (12) used in the main text. As described there, the basic
idea is to use is the fluctuations in a single mode of the electromagnetic
field with an external, homogeneous magnetic field as the detector.

To determine the interaction between the electromagnetic detector modes and
the gravitational field, we begin with the action for the coupling between
electromagnetism and gravity,

$S=\int
d^{4}x\,\sqrt{-g}\left(-\frac{1}{4}\,g^{\mu\alpha}g^{\nu\beta}F_{\mu\nu}F_{\alpha\beta}\right)\\!.$
(58)

Expanding the metric around flat spacetime
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$, one obtains the interaction to leading
order in the fluctuations $h$:

$\displaystyle S_{\rm int}$ $\displaystyle=\frac{1}{2}\int
d^{4}x\,h_{\mu\nu}T^{\mu\nu},$ (59) $\displaystyle T^{\mu\nu}$
$\displaystyle=F^{\mu\alpha}F^{\nu}_{\hphantom{\nu}\alpha}-\tfrac{1}{4}\eta^{\mu\nu}F_{\alpha\beta}F^{\alpha\beta},$

where here all indices are raised and lowered by $\eta_{\mu\nu}$. In the
Hamiltonian framework, we can write this as an interaction Hamiltonian

$V=-\frac{1}{2}\int_{\rm\ell^{3}}d^{3}{\mathbf{x}}\,h_{\mu\nu}T^{\mu\nu}.$
(60)

The integration is performed over the detector volume, $\ell^{3}$, which is
the region we assume our magnetic field is localized to.

Turning to the magnetic field, we assume it is homogeneous, has magnitude
$B_{0}$ and is oriented along the $z$-axis. We expand
$F_{\mu\nu}=F_{0,\mu\nu}+\delta F_{\mu\nu}$ into fluctuations around this
background,171717More precisely, we are taking the external magnetic field to
be static in the transverse traceless frame, whose coordinates are given by
freely falling observers. We make this assumption mainly for simplicity,
although it will not impact the qualitative conclusions we reach. In general,
the difference between frames for high-frequency gravitational waves must be
treated carefully, see for instance the discussion in Refs. [11, 12, 23].
Nevertheless, if the frequency exceeds the mechanical resonances of the
system, as is the case for the setups we consider, we expect this to be a good
approximation to a more complex detector model [24, 19]. where

$F_{0,\mu\nu}=\begin{pmatrix}0&0&0&0\\\ 0&0&B_{0}&0\\\ 0&-B_{0}&0&0\\\
0&0&0&0\end{pmatrix}\\!.$ (61)

In transverse traceless gauge, this reduces the interaction Hamiltonian to the
simple form

$\\!\\!V=B_{0}\int_{\ell^{3}}d^{3}{\mathbf{x}}\,\Big{[}h_{xz}\delta
B_{x}+h_{yz}\delta B_{y}-(h_{xx}+h_{yy})\delta B_{z}\Big{]}.$ (62)

This expression is completely general and only relies on the orientation of
the external magnetic field.

Now we want to determine the part of the coupling in Eq. (62) that drives the
electromagnetic mode that we will detect. We can decompose the electromagnetic
gauge field (in the Schrödinger picture) as

$\delta
A_{\mu}({\mathbf{x}})=\sum_{{\mathbf{p}},r}\epsilon_{r,\mu}({\mathbf{p}})u_{{\mathbf{p}}}({\mathbf{x}})a^{\dagger}_{{\mathbf{p}},r}+{\rm
h.c.},$ (63)

where $\epsilon_{r}({\mathbf{p}})$ are a complete set of polarization vectors,
and the $u_{{\mathbf{p}}}({\mathbf{x}})$ form a complete, discrete set of
modes:

$\int_{\ell^{3}}d^{3}{\mathbf{x}}\,\
u_{{\mathbf{p}}}({\mathbf{x}})u_{{\mathbf{p}}^{\prime}}^{*}({\mathbf{x}})=\frac{\delta_{{\mathbf{p}},{\mathbf{p}}^{\prime}}}{2\omega_{{\mathbf{p}}}},\
\ \
\epsilon_{r}({\mathbf{p}})\cdot\epsilon^{*}_{r^{\prime}}({\mathbf{p}})=\delta_{r,r^{\prime}}.$
(64)

In general, the nature of these modes depends on the details of the detector.
For a device like CAST, which detects X-ray wavelength photons
($\lambda\lesssim 1~{}{\rm nm}$) in a detector volume of order $m^{3}$, the
modes of interest can be modelled as travelling plane waves. In an experiment
like a microwave axion cavity, one instead detects a single mode whose
wavelength is the size of the cavity, and one should really use stationary
standing waves. For simplicity here, we will use plane waves,

$u_{{\mathbf{p}}}({\mathbf{x}})=\frac{e^{-i{\mathbf{p}}\cdot{\mathbf{x}}}}{\sqrt{2\omega_{{\mathbf{p}}}\ell^{3}}},\
\ \ {\mathbf{p}}=\frac{2\pi{\mathbf{n}}}{\ell},$ (65)

with ${\mathbf{n}}$ an arbitrary vector of integers, though this result can be
easily generalized to cavity modes with appropriate boundary conditions. The
normalizations are chosen so that the energy density in the field modes is

$H_{\rm
EM}=\frac{1}{2}\int_{\ell^{3}}d^{3}{\mathbf{x}}\left[{\mathbf{E}}^{2}+{\mathbf{B}}^{2}\right]=\sum_{{\mathbf{p}},r}\omega_{{\mathbf{p}}}a^{\dagger}_{{\mathbf{p}},r}a_{{\mathbf{p}},r},$
(66)

plus the usual infinite zero-point energy.

Meanwhile, the gravitational radiation field $h_{\mu\nu}$ is freely
propagating no matter what the detector boundary conditions are, so we expand
it as in Eq. (7), viz.,

$h_{\mu\nu}({\mathbf{x}})=M_{\rm
Pl}^{-1}\sum_{{\mathbf{k}},s}\epsilon_{s,\mu\nu}({\mathbf{k}})v_{{\mathbf{k}}}({\mathbf{x}})b_{{\mathbf{k}},s}^{\dagger}+{\rm
h.c.},$ (67)

where we take the polarization tensors $\epsilon_{s,\mu\nu}({\mathbf{k}})$ to
be real, and normalized according to,

$\epsilon_{s,\mu\nu}({\mathbf{k}})\,\epsilon_{s^{\prime}}^{\mu\nu}({\mathbf{k}})=\delta_{s,s^{\prime}}.$
(68)

Conventionally the two polarizations are labeled as
$s=+,\times$.181818Concretely, we work in the transverse traceless gauge, and
define the polarization tensors as
$\epsilon_{+}^{ij}=(u^{i}u^{j}-v^{i}v^{j})/\sqrt{2}$,
$\epsilon_{\times}^{ij}=(u^{i}v^{j}+v^{i}u^{j})/\sqrt{2}$ with
$\bf{v}=(\hat{\bf{e}}_{z}\times\hat{\bf{k}})/|\hat{\bf{e}}_{z}\times\hat{\bf{k}}|=\hat{\bf
e}_{\phi}$ and $\bf{u}=\bf{v}\times\hat{\bf{k}}$. The mode functions are just
the usual plane waves, normalized again to fix the correct expression for the
kinetic energy, and so that $h$ is a dimensionless strain:

$v_{{\mathbf{k}}}({\mathbf{x}})=\frac{e^{-i{\mathbf{k}}\cdot{\mathbf{x}}}}{\sqrt{2\omega_{{\mathbf{k}}}L^{3}}},\
\ \ {\mathbf{k}}=\frac{2\pi{\mathbf{n}}}{L}.$ (69)

To connect our notation with a classical gravitational wave of a single mode,
we would take amplitudes $b_{s}=h_{s}M_{\rm Pl}\sqrt{\omega L^{3}}$, so that

$h^{\mu\nu}(x)=\frac{1}{\sqrt{2}}\big{[}h_{+}\epsilon_{+}^{\mu\nu}({\mathbf{k}})+h_{\times}\epsilon_{\times}^{\mu\nu}({\mathbf{k}})\big{]}e^{ik\cdot
x}+{\rm h.c.},$ (70)

and the total strain is $h^{2}\equiv|h_{+}|^{2}+|h_{\times}|^{2}$.

Now, suppose we are monitoring a single mode $a=a_{{\mathbf{p}}_{0},r_{0}}$ of
the electromagnetic field, with wavevector
${\mathbf{p}}_{0}=\omega_{0}\hat{{\mathbf{x}}}$ and polarization
$r_{0}$.191919If the detector is an electromagnetic cavity, this can be done
easily by just isolating a single cavity frequency. For a detector with a
continuum of modes, one can isolate a single mode for example with filter
cavities near the final photodetectors. We use the single-mode language for
simplicity, but the results generalize in a straightforward way to finite-
bandwidth detectors. We take this to be a mode such that its magnetic field at
$t=0$ is aligned along the $y$-axis. Then only the $\delta B_{y}$ term in Eq.
(62) contributes. Using the mode expansion of Eq. (63) with
${\mathbf{B}}=\nabla\times{\mathbf{A}}$, and the gravitational mode expansion
Eq. (67), we find that the interaction Hamiltonian between the mode $a$ and
the gravitational field is given by

$V=\frac{B_{0}}{M_{\rm
Pl}}\int_{\ell^{3}}d^{3}{\mathbf{x}}\,({\mathcal{O}}_{b}+{\mathcal{O}}_{b}^{\dagger})({\mathcal{O}}_{a}+{\mathcal{O}}_{a}^{\dagger}),$
(71)

where we have defined the operators

${\mathcal{O}}_{b}=\sum_{{\mathbf{k}},s}\epsilon_{s,yz}({\mathbf{k}})v^{*}_{{\mathbf{k}}}({\mathbf{x}})b_{{\mathbf{k}},s},\quad{\mathcal{O}}_{a}=i\omega_{0}u^{*}_{{\mathbf{p}}_{0}}({\mathbf{x}})a.$
(72)

Equation (71) represents the precise coupling of the detector mode
$a=a_{{\mathbf{p}}_{0},r_{0}}$ to the entire gravitational field; notice that
it involves a sum over all the gravitational modes.

Finally, we can simplify the interaction potential by performing the volume
integral. We use the integral

$\displaystyle W({\mathbf{p}},{\mathbf{k}})=$
$\displaystyle\int_{-\ell/2}^{\ell/2}d^{3}{\mathbf{x}}\,e^{i({\mathbf{p}}-{\mathbf{k}})\cdot{\mathbf{x}}}$
(73) $\displaystyle=$
$\displaystyle\prod_{i=x,y,z}\frac{2\sin[(p^{i}-k^{i})\ell/2]}{p^{i}-k^{i}}.$

Again the precise nature of this window function depends on the relative size
of the detector and signal mode of interest. For simplicity, consider the case
where $\lambda_{\rm sig}\ll\ell$, and where we the detector integrates over a
reasonable number of signal periods $\omega\delta t\gg 1$ (as in, for example,
CAST). Then,
$W({\mathbf{p}},{\mathbf{k}})\to\ell^{3}\delta_{{\mathbf{p}},{\mathbf{k}}}$,
and we have simply

$V\rightarrow\frac{iB_{0}}{2\sqrt{2}M_{\text{Pl}}}\sqrt{\frac{\ell^{3}}{L^{3}}}\big{(}a^{\dagger}b-ab^{\dagger}\big{)},$
(74)

where $b=b_{{\mathbf{p}}_{0},\times}$. This completes the derivation of Eq.
(12) in the main text.

We end with a note on our final two assumptions. We have assumed the detector
is large compared to the signal wavelength; this assumption can be relaxed in
a straightforward way by just leaving the window functions $W$ under an
integral. Formally, the result in Eq. (74) also contains a contribution
proportional to $ab_{-{\mathbf{p}}_{0},\times}$ and its conjugate, which we
have dropped by assuming $\omega\delta t\gg 1$. In the case of a detector
which averages over many periods of the signal, these naturally average out in
the interaction Hamiltonian Eq. (14) (the “rotating wave approximation”). In a
slower detector, like a resonant amplitude detector, these terms do in fact
contribute by adding an additional quanta of graviton vacuum fluctuations. As
we will see in App. F, this additional quanta will not substantially change
any of our conclusions, and can safely be ignored.

##  Appendix F Variance calculations

In this appendix we provide the detailed computations of the variances in the
click rates – Eqs. (27) and (30) – and amplitude measurements – Eq. (32) –
given in the main text. Our discussion partly follows Ref. [31], and we refer
there for further details.

###  F.1 Number variance

Consider the intensity detector discussed in Sec. I.3. We begin with the
classical calculation. First, assume a classical incoming radiation state with
definite Fourier amplitude $b$. The mean number of clicks observed in a time
$T$ is then fixed by Eq. (22) to $\overline{N}=\eta\phi T$, with
$\phi=\phi(b)=|b|^{2}\ell^{2}/4L^{3}$ the incoming gravitational flux, as in
Eq. (23). Even so, the observed number of detector clicks will be stochastic.
As the events will each be independent, the observed number of clicks $N$ in a
single observation of time $T$ will be Poisson distributed,

$p(N|\,b)=\frac{1}{N!}[\eta\phi(b)T]^{N}e^{-\eta\phi(b)T},$ (75)

where we explicitly note this is the distribution in the semi-classical model.
We can confirm the average is as expected,

$\displaystyle\overline{N}$ $\displaystyle=\sum_{N=0}^{\infty}p(N|\,b)N$ (76)
$\displaystyle=\eta\phi(b)T\sum_{N=1}^{\infty}\frac{1}{(N-1)!}[\eta\phi(b)T]^{N-1}e^{-\eta\phi(b)T}$
$\displaystyle=\eta\phi(b)T.$

Similarly,

$\displaystyle\overline{N(N-1)}$
$\displaystyle=[\eta\phi(b)T]^{2}\sum_{N=2}^{\infty}\frac{1}{(N-2)!}[\eta\phi(b)T]^{N-2}e^{-\eta\phi(b)T}$
$\displaystyle=[\eta\phi(b)T]^{2},$ (77)

from which we conclude

$\displaystyle\Delta
N^{2}=\overline{N(N-1)}+\overline{N}-\overline{N}^{2}=\eta\phi(b)T,$ (78)

as one expects in a Poisson distribution. These are elementary calculations;
we include them as the calculations in the more general scenarios will mirror
these. Note that we use overlines to denote averages taken over many
observations.

Now instead of a definite value for $b$, suppose that it is a classically
random variable with distribution $P_{\rm cl}(b)$. Physically this could
correspond to a stochasticity in the amplitude of the input wave. By the law
of total probability, we now have

$\displaystyle p(N)=$ $\displaystyle\int db\,P_{\rm cl}(b)p(N|\,b)$ (79)
$\displaystyle=$ $\displaystyle\frac{1}{N!}\int db\,P_{\rm
cl}(b)[\eta\phi(b)T]^{N}e^{-\eta\phi(b)T}.$

Despite the superficial similarity, for a general $P_{\rm cl}(b)$ Eq. (79) is
not the Poisson distribution. This can be confirmed by explicit calculation.
We have

$\displaystyle\overline{N}$ $\displaystyle=\int db\,P_{\rm
cl}(b)\sum_{N=0}^{\infty}p(N|\,b)N$ (80) $\displaystyle=\int db\,P_{\rm
cl}(b)\,\eta\phi(b)T$ $\displaystyle=\eta T\\!\braket{\phi}\\!,$

directly from Eq. (76). Here we use the notation $\braket{\cdot}$ to denote an
integral over $P_{\rm cl}(b)$ is performed, for example $\braket{\phi}=\int
dbP_{\rm cl}(b)\phi(b)$. We emphasize that this is distinct from expectation
values over the distribution of the clicks observed in a time $T$, which we
represent with an overline. Continuing, Eq. (77) implies

$\displaystyle\overline{N(N-1)}=\int dbP_{\rm
cl}(b)[\eta\phi(b)T]^{2}=\eta^{2}T^{2}\\!\braket{\phi^{2}}\\!,$ (81)

so that the variance is now

$\displaystyle\Delta
N^{2}=\overline{N}+\eta^{2}T^{2}\big{\langle}(\phi-\braket{\phi})^{2}\big{\rangle}\geq\bar{N},$
(82)

as quoted in the main text in Eq. (27).

These same arguments can now be generalized to the case where the incident
graviton flux is quantized. Consider first an incoming field of definite
coherent state $\ket{\beta}$, such that the flux takes a fixed value of
$\phi=\phi(\beta)=|\beta|^{2}\ell^{2}/4L^{3}$. Clicks in the detector will
remain independent, so that again the observed number will be Poisson
distributed,

$p(N|\,\beta)=\frac{1}{N!}[\eta\phi(\beta)T]^{N}e^{-\eta\phi(\beta)T}.$ (83)

The mean and variance are identical to the fixed classical case with
$|b|^{2}\to|\beta|^{2}$.

A general quantized initial state is described by a Glauber $P(\beta)$, and
although this is not a probability distribution, the general distribution is
given by

$p(N)=\frac{1}{N!}\int
d\beta\,P(\beta)\,[\eta\phi(\beta)T]^{N}e^{-\eta\phi(\beta)T}.$ (84)

A careful derivation justifying this form for a stationary incident field is
given in Ref. [31]. Crucially, however, Eq. (84) is _not_ the quantum
expectation value of Eq. (75) in the state $\ket{\beta}$, with
$|b|^{2}=b^{\dagger}b$ interpreted as an operator. Rather, it is the
expectation value of the _normal ordered_ version of $p(N|b)$. This latter
statement follows from a general result sometimes called the optical
equivalence theorem, which says that for an operator
$\mathcal{O}(b,b^{\dagger})$, we have
$\braket{:\\!\mathcal{O}(b,b^{\dagger})\\!:}=\int
d\beta\,P(\beta)\mathcal{O}(\beta,\beta^{*})$. This normal ordering will be
crucial in what follows. Note that $\braket{\phi}=\braket{:\\!\phi\\!:}$ since
$\phi\sim b^{\dagger}b$ is already normal ordered.

Again, for a general $P(\beta)$, Eq. (84) is not the Poisson distribution. The
deviations can be even more striking than in the classical case. Given Eq.
(84), one can proceed as above, and finds in similar fashion

$\displaystyle\overline{N}$ $\displaystyle=\eta T\\!\braket{\phi}\\!,$ (85)
$\displaystyle\Delta N^{2}$ $\displaystyle=\overline{N}+\eta^{2}T^{2}\int
d\beta\,P(\beta)\left[\phi(\beta)-\braket{\phi}\right]^{2}\\!,$

with the second result confirming Eq. (30). As the integral is now performed
over $P(\beta)$, which can take on negative values, $\Delta
N^{2}-\overline{N}$ can take on either sign, and exhibit sub-Poisson
statistics, which are strictly forbidden for a purely classical incident
state. An explicit example of this was provided in App. C. The result cannot
be arbitrarily negative, however. As $\phi$ is proportional to the number
operator,

$\displaystyle\Delta N^{2}$
$\displaystyle=\overline{N}+\eta^{2}T^{2}\left[\langle:\\!\phi^{2}\\!:\rangle-\langle:\\!\phi\\!:\rangle^{2}\right]$
(86) $\displaystyle=\bar{N}(1-\eta
T\ell^{2}/4L^{3})+\eta^{2}T^{2}\braket{\Delta\phi^{2}}$
$\displaystyle\geq\bar{N}(1-\eta T\ell^{2}/4L^{3}).$

As expected, the size of the deviation is controlled by $\eta$.

###  F.2 Amplitude variance

Consider next the amplitude detector discussed in Sec. I.4. We will compute
its output noise spectrum,

$\Delta X^{2}=\overline{X^{2}}-\overline{X}^{2}$ (87)

in the presence of some incoming gravitational radiation. Here the overline
denotes an expectation value taken over $p(X)$. We will consider the classical
and quantum calculations as well as the fixed and variable incident flux
scenarios simultaneously. We take the interaction as in Eq. (12), namely

$V=ig(a^{\dagger}b-ab^{\dagger}),$ (88)

where for the time being $b,\,b^{\dagger}$ can be either c-numbers or
operators. To compute the variance, we will need the Krauss operators and POVM
to second order in $g\delta t\ll 1$. The Krauss operator to this order is

$\displaystyle K_{X}$ $\displaystyle=\braket{X}{0}+g\delta
t\braket{X}{a^{\dagger}}{0}b$ (89) $\displaystyle+\frac{1}{2}g^{2}\delta
t^{2}\left[\braket{X}{a^{{\dagger}2}}{0}b^{2}-\braket{X}{aa^{\dagger}}{0}b^{\dagger}b\right]\\!.$

Note again that $\ket{0}$ means the vacuum $\ket{n=0}$, not the position
eigenstate $\ket{X=0}$. This gives the POVM elements to the same order:

$\displaystyle E_{X}$ $\displaystyle=K_{X}^{\dagger}K_{X}$
$\displaystyle=\braket{0}{X}\braket{X}{0}$ $\displaystyle+g\delta
t\left[\braket{0}{X}\braket{X}{a^{\dagger}}{0}b+\braket{0}{a}{X}\braket{X}{0}b^{\dagger}\right]$
$\displaystyle+\frac{1}{2}g^{2}\delta
t^{2}\left[\braket{0}{X}\braket{X}{a^{{\dagger}2}}{0}b^{2}-\braket{0}{X}\braket{X}{aa^{\dagger}}{0}b^{\dagger}b\right]$
$\displaystyle+\frac{1}{2}g^{2}\delta
t^{2}\left[\braket{0}{a^{2}}{X}\braket{X}{0}b^{{\dagger}2}-\braket{0}{a^{\dagger}a}{X}\braket{X}{0}b^{\dagger}b\right]$
$\displaystyle+g^{2}\delta
t^{2}\left[\braket{0}{a}{X}\braket{X}{a^{\dagger}}{0}b^{\dagger}b\right]\\!.$
(90)

At this stage, the dependence on the detector is reduced to numerical factors
involving matrix elements of creation and annihilation operators
$a,a^{\dagger}$. The dependence on the gravitational signal enters through the
$b,b^{\dagger}$.

Now we can assume some state (either classical or quantum mechanical) for the
gravitational signal, and use it to compute the average value of the detector
amplitude output $\overline{X}$, by

$\overline{X}=\int dX\braket{E_{X}}X.$ (91)

This follows as in the notation adopted, $p(X)=\braket{E_{X}}$, which is the
appropriate generalization of Eq. (19). The integral over $X$ reduces the
inner products in Eq. (90) to simple numbers, via manipulations of the form

$\int dX\braket{0}{X}\braket{X}{0}X=\braket{0}{X}{0}=0.$ (92)

One finds after straightforward algebra of this type that

$\overline{X}=\frac{1}{\sqrt{2}}g\delta t\braket{b+b^{\dagger}}=g\delta
t\braket{X_{b}}.$ (93)

This final expectation value can be treated either classically, by integrating
over a distribution $P_{\rm cl}(b)$, or quantum mechanically, by integrating
in the Glauber representation $P(\beta)$. A similar calculation results in

$\overline{X^{2}}=\frac{1}{2}+\frac{1}{2}g^{2}\delta
t^{2}\braket{b^{2}+b^{{\dagger}2}+2b^{\dagger}b}.$ (94)

This result is valid both classically and quantum mechanically; to obtain it
we did not need to invoke any commutator between $b,\,b^{\dagger}$. The first
$1/2$ term comes from the commutator $[a,a^{\dagger}]=1$, i.e., the detector’s
vacuum fluctuations in the quadrature variables $X,Y$.

The key difference between quantum and classical radiation variables now
appears when we try to re-express this result in terms of
$X_{b}=(b+b^{\dagger})/\sqrt{2}$. In the quantum case we need to use a
commutator $[b,b^{\dagger}]=1$, and in particular one has

$\displaystyle\overline{X^{2}}_{\rm cl}$
$\displaystyle=\frac{1}{2}+g^{2}\delta t^{2}\braket{X_{b}^{2}}\\!,$ (95)
$\displaystyle\overline{X^{2}}_{\rm qu}$
$\displaystyle=\frac{1}{2}+g^{2}\delta
t^{2}\left[\braket{X_{b}^{2}}-\frac{1}{2}\right]\\!.$

Putting these results together, we obtain the variance in the as

$\Delta X^{2}=\frac{1}{2}+g^{2}\delta t^{2}\begin{cases}\braket{\Delta
X_{b}^{2}}\\!,&{\rm classical},\\\ \big{[}\\!\braket{\Delta
X_{b}^{2}}\\!-\frac{1}{2}\big{]},&{\rm quantum}.\end{cases}$ (96)

In terms of explicit integration over a classical probability or quantum
Glauber representation,

$\Delta X^{2}=\frac{1}{2}+2g^{2}\delta t^{2}\begin{cases}\int db\,P_{\rm
cl}(b)\left[{\rm Re}~{}b-\braket{{\rm Re}~{}b}\right]^{2}\\!,&{\rm cl}.\\\
\int d\beta\,P(\beta)\left[{\rm Re}~{}\beta-\braket{{\rm
Re}~{}\beta}\right]^{2}\\!,&{\rm qu}.\end{cases}$ (97)

In the classical case the expectation value is taken over $P_{\rm cl}(b)$, and
therefore $\Delta X^{2}\geq 1/2$. In the quantum case, the integral is over
$P(\beta)$, which again can allow us to evade the classical bound. However,
from Eq. (96) we see that even for an extremely squeezed gravitational state
that achieves $\braket{\Delta X_{b}^{2}}\to 0$, we only have demonstrably
quantum behavior in the range $1>2\Delta X^{2}\geq 1-g^{2}\delta t^{2}$.
Again, although we see a distinctly quantum state is required, it is
insufficient: we will also need a detector capable of registering the
$g^{2}\delta t^{2}$ effect. Notice that a large-amplitude coherent state does
not make this easier.