QUANTUM OPTIMALITY OF PHOTON COUNTING FOR TEMPERATURE MEASUREMENT OF THERMAL ASTRONOMICAL SOURCES

Regardless of the frequency range, the electromagnetic field input to an antenna or telescope is described by a time-dependent positive-frequency field operator $E(t)$ (in units of $\sqrt{{\rm{photons}}\;{{\rm{s}}}^{-1}}$) given by

$$\begin{eqnarray}&&E(t)={\displaystyle \int }_{0}^{\infty }{a}_{\nu }\;{e}^{-i2\pi \nu t}\;d\nu .\end{eqnarray} \tag{ 4 }$$

Here, the ${\{{a}_{\nu }\}}_{\nu \gt 0}$ are un-normalized single-frequency annihilation operators satisfying $\left[{a}_{\nu },{a}_{{\nu }^{\prime }}^{\dagger }\right]=\delta (\nu -{\nu }^{\prime })$ which implies that $\left[E(t),{E}^{\dagger }(t\prime )\right]=\delta (t-t\prime )$. As in Lieu et al. (2015 and Zmuidzinas (2015), we are assuming the input field to be of a single polarization and in a single spatial mode to focus on the central issue—the additional generality does not substantially alter the result. If the input field is from a thermal source, it is a Gaussian field (Shapiro 2009, Section III) with the quantum expectation values (Mandel & Wolf 1995; Zmuidzinas 2003a):

$$\begin{eqnarray}&&\left\langle{a}_{\nu }\right\rangle=0,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\left\langle{a}_{\nu }\;{a}_{{\nu }^{\prime }}\right\rangle=0\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\left\langle{a}_{{\nu }^{\prime }}^{\dagger }\;{a}_{\nu }\right\rangle={n}_{{\mathsf{th}}}(\nu )\delta (\nu -{\nu }^{\prime }),\end{eqnarray} \tag{ 7 }$$

where ${n}_{{\mathsf{th}}}(\nu )$ is the mean occupation number in a thermal state given by the Planck formula

$$\begin{eqnarray}&&{n}_{{\mathsf{th}}}(\nu )=\displaystyle \frac{1}{{e}^{\frac{h\nu }{{{kT}}_{{\rm{s}}}}}-1},\end{eqnarray} \tag{ 8 }$$

for T_s the source temperature (in K) and k, Boltzmann's constant. In terms of field operators, these relations imply

$$\begin{eqnarray}&&\langle E(t)\rangle =0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{K}^{(p)}(t,t\prime ):= \langle E(t)E(t\prime )\rangle =0,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{K}^{(n)}(t,t\prime ):= \left\langle{E}^{\dagger }(t)E(t\prime )\right\rangle={\displaystyle \int }_{0}^{\infty }{n}_{{\mathsf{th}}}(\nu ){e}^{-i2\pi \nu (t\prime -t)}\;d\nu .\end{eqnarray} \tag{ 11 }$$

Here, the functions ${K}^{(p)}(t,t\prime )$ and ${K}^{(n)}(t,t\prime )$ are the phase-sensitive and (normally ordered) phase-insensitive correlation functions of the field respectively (see Equations (61)–(62) of Shapiro 2009). Note that these functions (as well as the mean) depend only upon the time difference $\tau := t\prime -t$, indicating the statistical stationarity of the field. Using the Gaussian moment-factoring theorem (Mandel & Wolf 1995), all higher moments of the field operators can be expressed in terms of these functions, which therefore constitute a complete description of the field.

The thermal field described above is broadband. In practice, a measurement operates on only a finite band of the input field that is determined either by insertion of filters or the response of the measuring device. We will accordingly assume that the input is passed through a filter centered at frequency ${\nu }_{0}$ and with a flat profile over the band $[{\nu }_{0}-{\rm{\Delta }}\nu /2,{\nu }_{0}+{\rm{\Delta }}\nu /2]$. In a typical radio astronomy measurement, e.g., we may have ${\nu }_{0}=1$ GHz, and ${\rm{\Delta }}\nu$ a few MHz. Further, in radio astronomy, the source temperature T_s is typically such that $h{\nu }_{0}/{{kT}}_{{\rm{s}}}\ll 1$, so that the Rayleigh–Jeans approximation to Equation (8) may be used (Burke & Graham-Smith 2010). As a further simplification, we may assume that the mean occupation number of Equation (8) is approximately flat in the band $[{\nu }_{0}-{\rm{\Delta }}\nu /2,{\nu }_{0}+{\rm{\Delta }}\nu /2]$ at the value ${n}_{0}\simeq \ {{kT}}_{{\rm{s}}}/h{\nu }_{0}$. As a result of the above assumptions, for the field at the output of the filter, we replace ${n}_{{\mathsf{th}}}(\nu )$ in Equations (7) and (11) by

$$\begin{eqnarray}n(\nu )=\left\{\begin{array}{ll}{n}_{0}=\displaystyle \frac{{{kT}}_{{\rm{s}}}}{h{\nu }_{0}} & {\rm{if}}\;{\rm{}}\nu \in [{\nu }_{0}-{\rm{\Delta }}\nu /2,{\nu }_{0}+{\rm{\Delta }}\nu /2]\\ 0 & {\rm{otherwise}}.\end{array}\right.\end{eqnarray} \tag{ 12 }$$

In this quasi-monochromatic regime, $n(\nu )$ is essentially the (dimensionless) power spectral density of the field, in that multiplication by $h{\nu }_{0}$ gives the average power per unit frequency (in ${\rm{W}}\;{{\rm{Hz}}}^{-1}$) of the field.

The field at the output of the bandpass filter continues to satisfy Equations (9) and (10), but Equation (11) is modified to

$$\begin{eqnarray}&&{K}^{(n)}(t,t\prime )\equiv {K}^{(n)}(\tau )={n}_{0}\cdot {\rm{\Delta }}\nu \;\mathrm{sinc}\;[{\rm{\Delta }}\nu \tau ]\;{e}^{-i2\pi {\nu }_{0}\tau },\end{eqnarray} \tag{ 13 }$$

where $\mathrm{sinc}\;(x)=\mathrm{sin}(\pi x)/(\pi \;x)$ is the sinc function. The coherence time ${\tau }_{{\rm{c}}}$ of the output radiation (Mandel 1959), which is the approximate "width" of $\left|{K}^{(n)}(\tau )\right|$, is then

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}:= {\displaystyle \int }_{-\infty }^{\infty }{\left|{g}^{(1)}(\tau )\right|}^{2}\;d\tau ={\displaystyle \int }_{-\infty }^{\infty }{\left|\displaystyle \frac{{K}^{(n)}(\tau )}{{K}^{(n)}(0)}\right|}^{2}\;d\tau =\displaystyle \frac{1}{{\rm{\Delta }}\nu },\end{eqnarray} \tag{ 14 }$$

and is also the separation between the peak of $\left|{K}^{(n)}(\tau )\right|$ at $\tau =0$ and the first zero. For the radio frequency example above, we have ${\tau }_{{\rm{c}}}\sim 1\;\mu {\rm{s}}$. In the optical regime, ${\tau }_{{\rm{c}}}$ is typically much smaller even downstream of an optical filter, on the of order of nanoseconds.

The measurement on the field takes place in a finite time interval ${\mathcal{T}}$ of duration T_s, which we take to be $[-T/2,T/2]$. In the remainder of the paper, we assume that the parameters ${\nu }_{0},{\rm{\Delta }}\nu$ (or ${\tau }_{{\rm{c}}}$), and T are fixed and known, while n₀ (or equivalently, T_s) is the single unknown parameter that we wish to estimate. We further assume that ${\tau }_{{\rm{c}}}\ll T$, which is usually the case for the detection of faint sources, for which T could range from seconds to hours to days.

Quantum information and metrology using electromagnetic fields (Helstrom 1976; Holevo 2011), and Gaussian quantum information in particular (Olivares 2012; Weedbrook et al. 2012), is usually set up and studied on a finite set of modes that are excited in the problem under consideration. In order to obtain a modal description of the field from that in terms of its mean and correlation functions in the previous subsection, we expand the field in the measurement interval ${\mathcal{T}}$ in terms of a complete orthonormal set of positive-frequency "Fourier-series" traveling-wave modes given by

$$\begin{eqnarray}{\phi }_{m}(t)=\left\{\begin{array}{ll}\displaystyle \frac{1}{\sqrt{T}}\mathrm{exp}\left(-i\displaystyle \frac{2\pi {mt}}{T}\right) & {\rm{if}}\;{\rm{}}t\in [-T/2,T/2]\\ 0 & {\rm{otherwise,}}\end{array}\right.\end{eqnarray} \tag{ 15 }$$

for $m=0,1,2,\ldots$. The modal annihilation operators

$$\begin{eqnarray}&&{a}_{m}:= {\displaystyle \int }_{-\infty }^{\infty }E(t){\phi }_{m}^{*}(t){dt},\;\;\;\;m=0,1,2,\ldots \end{eqnarray} \tag{ 16 }$$

then satisfy the (normalized) commutation relations

$$\begin{eqnarray}&&\left[{a}_{m},{a}_{n}^{\dagger }\right]={\delta }_{{mn}}.\end{eqnarray} \tag{ 17 }$$

As in standard quantum mechanics, the field state is described by a density operator ρ (a positive semidefinite operator with $\mathrm{Tr}\;\rho =1$) on a Hilbert space ${\mathcal{H}}={\displaystyle \otimes }_{m=0}^{\infty }{{\mathcal{H}}}_{m},$ where ${{\mathcal{H}}}_{m}$ is the (infinite-dimensional) Hilbert space of the mth mode. Our next task is to obtain ρ.

Since E(t) is a Gaussian field, the density operator ρ is in a so-called Gaussian state (Shapiro 2009; Holevo 2011). In order to define Gaussian states, we need some notation. For each m, we define the quadrature operators

$$\begin{eqnarray}&&{q}_{m}=\displaystyle \frac{{a}_{m}+{a}_{m}^{\dagger }}{\sqrt{2}};\;\;{p}_{m}=\displaystyle \frac{{a}_{m}-{a}_{m}^{\dagger }}{\sqrt{2}i}\end{eqnarray} \tag{ 18 }$$

satisfying the canonical commutation relations $[{q}_{m},{q}_{n}]=[{p}_{m},{p}_{n}]=0$ and $\left[{q}_{m},{p}_{n}\right]=i{\delta }_{{mn}}$. Consider the vector ${\boldsymbol{R}}\equiv ({R}_{1},{R}_{2},\ldots ):= ({q}_{1},{p}_{1},{q}_{2},{p}_{2},\ldots )$ of quadrature operators. The mean vector $\bar{{\boldsymbol{R}}}$ in the state ρ is

$$\begin{eqnarray}&&\bar{{\boldsymbol{R}}}:= {\langle {\boldsymbol{R}}\rangle }_{\rho }=(\mathrm{Tr}\;(\rho {R}_{1}),\mathrm{Tr}\;(\rho {R}_{2}),\ldots )\end{eqnarray} \tag{ 19 }$$

and the covariance matrix ${\boldsymbol{\sigma }}$ has the ($i,j$)th matrix element

$$\begin{eqnarray}{\sigma }_{{ij}} & := & \displaystyle \frac{1}{2}\langle ({R}_{i}-{\bar{R}}_{i})({R}_{j}-{\bar{R}}_{j})+({R}_{j}-{\bar{R}}_{j})({R}_{i}-{\bar{R}}_{i}){\rangle }_{\rho }\\ & = & {\rm{Re}}\left\{\mathrm{Tr}\;\rho ({R}_{i}-{\bar{R}}_{i})({R}_{j}-{\bar{R}}_{j})\right\}.\end{eqnarray} \tag{ 20 }$$

With this notation, a Gaussian state is a state whose Wigner characteristic function

$$\begin{eqnarray}&&{\chi }_{\rho }({\boldsymbol{\xi }}):= \mathrm{Tr}\;\left[\rho \;\mathrm{exp}\left(-{{\boldsymbol{\xi }}}^{T}\;{\rm{\Omega }}\;{\boldsymbol{R}}\right)\right]\end{eqnarray} \tag{ 21 }$$

is of the Gaussian form (Olivares 2012)

$$\begin{eqnarray}&&{\chi }_{\rho }({\boldsymbol{\xi }})=\mathrm{exp}\left(-i{{\boldsymbol{\xi }}}^{T}\;{\rm{\Omega }}\;\bar{{\boldsymbol{R}}}-\displaystyle \frac{1}{2}{{\boldsymbol{\xi }}}^{T}{\rm{\Omega }}{\boldsymbol{\sigma }}{{\rm{\Omega }}}^{T}{\boldsymbol{\xi }}\right),\end{eqnarray} \tag{ 22 }$$

where ${\boldsymbol{\xi }}=({\xi }_{1}^{(1)},{\xi }_{2}^{(1)},{\xi }_{1}^{(2)},{\xi }_{2}^{(2)},\ldots )$ is the vector of Fourier variables corresponding to the phase-space coordinates of each mode. The matrix

$$\begin{eqnarray}&&{\rm{\Omega }}=\underset{m}{\displaystyle \oplus }{\boldsymbol{\omega }}\end{eqnarray} \tag{ 23 }$$

is block-diagonal in the 2 × 2 blocks

$$\begin{eqnarray}{\boldsymbol{\omega }}=\left(\begin{array}{cc}0 & 1\\ -1 & 0\end{array}\right),\end{eqnarray} \tag{ 24 }$$

with one block per mode.

From Equation (22), we see that a Gaussian state is completely described by its mean and covariance matrix. To calculate these for filtered thermal radiation, we must compute the first- and second-order moments of the form $\left\langle{a}_{m}\right\rangle$, $\left\langle{a}_{m}\;{a}_{n}\right\rangle$, and $\left\langle{a}_{m}^{\dagger }\;{a}_{n}\right\rangle$ for each m and n. In the Appendix, we present the detailed calculations for the above quantities in the long observation time limit ${\tau }_{{\rm{c}}}\ll T$. The results are

$$\begin{eqnarray}&&\left\langle{a}_{m}\right\rangle=0\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\left\langle{a}_{m}\;{a}_{n}\right\rangle=0\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\left\langle{a}_{m}^{\dagger }\;{a}_{n}\right\rangle\simeq {n}_{0}\;{\rm{rect}}\;\left[\displaystyle \frac{{\nu }_{0}-\displaystyle \frac{m}{T}}{{\rm{\Delta }}\nu }\right]\cdot {\delta }_{{mn}},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}{\rm{rect}}(x)=\left\{\begin{array}{ll}1 & {\rm{if}}\;{\rm{}}| x| \leqslant 1/2\\ 0 & {\rm{otherwise}}\end{array}\right.\end{eqnarray} \tag{ 28 }$$

is the rectangle function. In other words, in this ${\tau }_{{\rm{c}}}\ll T$ limit, we have $\left\langle{a}_{m}^{\dagger }\;{a}_{n}\right\rangle\simeq 0$ for $m\ne n$ and the average number of photons

$$\begin{eqnarray}&&\left\langle{a}_{m}^{\dagger }\;{a}_{m}\right\rangle\simeq {n}_{0}\;\;{\rm{for}}\;\;m\in {\mathcal{M}}=\left\{\left({\nu }_{0}-\displaystyle \frac{{\rm{\Delta }}\nu }{2}\right)T,\right.\\ &&\quad \left.\left({\nu }_{0}-\displaystyle \frac{{\rm{\Delta }}\nu }{2}\right)T+1,\ldots ,\left({\nu }_{0}+\displaystyle \frac{{\rm{\Delta }}\nu }{2}\right)T\right\}\end{eqnarray} \tag{ 29 }$$

i.e., in $M:= T{\rm{\Delta }}\nu$ "approximately single-frequency" modes that are within the bandwidth of the filter. Together with Equations (25)–(26), Equation (29) implies that the covariance matrix ${\boldsymbol{\sigma }}$ of ρ is

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}=\underset{m\in {\mathcal{M}}}{\displaystyle \oplus }\left(\displaystyle \frac{2{n}_{0}+1}{2}\;{{\mathbb{1}}}_{2}\right)\;\underset{m\notin {\mathcal{M}}}{\displaystyle \oplus }\left({{\mathbb{1}}}_{2}\right),\end{eqnarray} \tag{ 30 }$$

where ${{\mathbb{1}}}_{2}$ is the 2 × 2 identity matrix. This implies that $M={\rm{\Delta }}\nu T$ of the modes (those in the set ${\mathcal{M}}$ defined in Equation (29)) are in independent thermal states each with n₀ photons on average given by (Olivares 2012):

$$\begin{eqnarray}{\rho }_{{\mathsf{th}}}({n}_{0}) & = & \displaystyle \frac{1}{{n}_{0}+1}\underset{k=0}{\overset{\infty }{\displaystyle \sum }}{\left(\displaystyle \frac{{n}_{0}}{{n}_{0}+1}\right)}^{k}| k\rangle \langle k| \\ & = & \displaystyle \frac{1}{\pi \;{n}_{0}}{\displaystyle \int }_{{\mathbb{C}}}\mathrm{exp}\left(-\displaystyle \frac{{| \alpha | }^{2}}{{n}_{0}}\right)| \alpha \rangle \langle \alpha | \;{d}^{2}\alpha .\end{eqnarray} \tag{ 31 }$$

In the first representation, ${\rho }_{{\mathsf{th}}}({n}_{0})$ is a mixture of number states with a Bose–Einstein distribution, while in the second, it is a zero-mean circularly symmetric Gaussian distribution of coherent states $\{| \alpha \rangle \}$, i.e., of the eigenstates of the annihilation operator—$a| \alpha \rangle =\alpha \;| \alpha \rangle$—of the relevant mode. The remaining modes are all in the vacuum state $| 0\rangle \langle 0|$. Thus, the overall state is

$$\begin{eqnarray}&&{\rho }_{{n}_{0}}=\left(\underset{m\in {\mathcal{M}}}{\displaystyle \otimes }{\rho }_{{\mathsf{th}}}^{(m)}({n}_{0})\right)\left(\underset{m\notin {\mathcal{M}}}{\displaystyle \otimes }{| 0\rangle }^{(m)(m)}\langle 0| \right).\end{eqnarray} \tag{ 32 }$$

Since the modes $m\notin {\mathcal{M}}$ carry no information about n₀, it is sufficient to make measurements on just the modes in the set ${\mathcal{M}}$, effectively reducing the problem to one involving a finite number of modes.

The classical Cramér–Rao bound (Rao 1945; Cramér 1946; Van Trees 2001) provides a lower bound on the variance ${\rm{Var}}\;\hat{\theta }$ of any unbiased estimator $\hat{\theta }$ of an unknown parameter θ indexing a family of probability distributions $\{{P}_{\theta }\}$ on a given sample space. An unbiased estimator of θ is one that satisfies ${\mathsf{E}}[\hat{\theta }]=\theta$, where the statistical expectation value is taken with respect to ${P}_{\theta }$. The Cramér–Rao bound is widely used to provide limits on and benchmarks for the performance of communication and measurement systems (Van Trees 2001). It has also found application in the design of astronomical instruments (Zmuidzinas 2003b) and in measurements of the cosmic microwave background (Yadav et al. 2007).

In a quantum estimation problem, instead of probability distributions $\{{P}_{\theta }\}$, we are provided with a family $\{{\rho }_{\theta }\}$ of density operators of a given quantum-mechanical system depending on the unknown parameter. The additional feature of the quantum estimation problem over its classical counterpart is the freedom of choosing the quantum measurement that generates a probability distribution from ${\rho }_{\theta }$. All possible quantum measurements can be mathematically described by an object called a positive-operator-valued measure (POVM) (Helstrom 1976; Holevo 2011), which subsumes the well-known observables of standard quantum mechanics. The quantum Cramér–Rao bound (henceforth "q-CR bound") (Helstrom 1967, 1968, 1973, 1976; Holevo 2011) provides a lower bound on the variance of any unbiased estimator $\hat{\theta }$ of a parameter θ indexing a family of density operators $\{{\rho }_{\theta }\}$ optimized over all possible POVMs subject to the unbiasedness condition.

We simply state the result of the q-CR bound here, referring to Helstrom (1967, 1968, 1973, 1976) and Holevo (2011) for details. We are given a family $\{{\rho }_{\theta }\}$ of density operators depending on the parameter of interest θ. The operator equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{\theta }}{\partial \theta }=\displaystyle \frac{1}{2}\left({L}_{\theta }{\rho }_{\theta }+{\rho }_{\theta }{L}_{\theta }\right),\end{eqnarray} \tag{ 33 }$$

has a unique Hermitian solution ${L}_{\theta }={L}_{\theta }^{\dagger }$ when ${\rho }_{\theta }$ has no zero eigenvalues (Bhatia 2007). The operator ${L}_{\theta }$ is called the symmetric logarithmic derivative (SLD) in analogy with the classical case. The quantity

$$\begin{eqnarray}&&{I}_{Q}(\theta )=\mathrm{Tr}\;\left({\rho }_{\theta }\;{L}_{\theta }^{2}\right)={\left\langle{L}_{\theta }^{2}\right\rangle}_{{\rho }_{\theta }}\end{eqnarray} \tag{ 34 }$$

is known as the quantum Fisher information and the q-CR bound reads

$$\begin{eqnarray}&&{\rm{Var}}\;\hat{\theta }\geqslant \displaystyle \frac{1}{{I}_{Q}(\theta )},\end{eqnarray} \tag{ 35 }$$

and is valid for any unbiased estimator satisfying ${\mathsf{E}}[\hat{\theta }]=\theta$, where the expectation is over the probability distribution induced by the POVM on the state ${\rho }_{\theta }$.

The q-CR bound was originally developed in the context of quantum optical communication in the years following the invention of the laser, so it was natural for the early work to focus on the experimentally important Gaussian states of light, particularly on the estimation of the mean vector (Equation (19)) of Gaussian states. The problem of estimating the average photon number in a thermal state was also considered by Helstrom (1968). Very recently, the estimation of a general parameter indexing single-mode and multi-mode Gaussian states has been considered by Pinel et al. (2013) and Monras (2013), respectively.

With the density operator for filtered thermal radiation now at hand from our modal decomposition in Section 2.2, we may invoke the required q-CR bound from Helstrom (1968), Monras (2013), and Pinel et al. (2013). For completeness, however, we re-derive the q-CR bound for estimating the average photon number n₀, to begin with, in a single-mode thermal state ${\rho }_{{\mathsf{th}}}({n}_{0})$. From the first representation in Equation (31), we see that

$$\begin{eqnarray}&&{\rho }_{{\mathsf{th}}}({n}_{0})=\displaystyle \frac{1}{{n}_{0}+1}\;{\left(\displaystyle \frac{{n}_{0}}{{n}_{0}+1}\right)}^{N},\end{eqnarray} \tag{ 36 }$$

where $N={a}^{\dagger }\;a={\displaystyle \sum }_{k=0}^{\infty }k| k\rangle \langle k|$ is the number operator. Since ${n}_{0}\gt 0$, ${\rho }_{{\mathsf{th}}}({n}_{0})$ has no zero eigenvalues and a unique SLD exists. To find it, we compute the derivative

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\mathsf{th}}}({n}_{0})}{\partial {n}_{0}}=\displaystyle \frac{1}{{({n}_{0}+1)}^{2}}\left[\displaystyle \frac{N}{{n}_{0}}{\left(\displaystyle \frac{{n}_{0}}{{n}_{0}+1}\right)}^{N}-{\left(\displaystyle \frac{{n}_{0}}{{n}_{0}+1}\right)}^{N}\right]\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&={\rho }_{{\mathsf{th}}}({n}_{0})\left[\displaystyle \frac{N}{{n}_{0}({n}_{0}+1)}-\displaystyle \frac{1}{{n}_{0}+1}\right].\end{eqnarray} \tag{ 38 }$$

Comparing with Equation (33), we obtain the SLD operator

$$\begin{eqnarray}&&{L}_{{n}_{0}}=\displaystyle \frac{N}{{n}_{0}({n}_{0}+1)}-\displaystyle \frac{1}{{n}_{0}+1}.\end{eqnarray} \tag{ 39 }$$

Note that ${L}_{{n}_{0}}$ commutes with ${\rho }_{{\mathsf{th}}}({n}_{0})$, reflecting the fact that the $\{{\rho }_{{\mathsf{th}}}({n}_{0})\}$ commute with each other. The quantum Fisher information is

$$\begin{eqnarray}&&{I}_{Q}({n}_{0})={\left\langle{L}_{{n}_{0}}^{2}\right\rangle}_{{\rho }_{{\mathsf{th}}}({n}_{0})}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{1}{{({n}_{0}+1)}^{2}}\left[{\left\langle\displaystyle \frac{{N}^{2}}{{n}_{0}^{2}}-\displaystyle \frac{2N}{{n}_{0}}+1\right\rangle}_{{\rho }_{{\mathsf{th}}}({n}_{0})}\right]\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{1}{{({n}_{0}+1)}^{2}}\left[\displaystyle \frac{2{\langle N\rangle }_{{\rho }_{{\mathsf{th}}}({n}_{0})}^{2}+{\langle N\rangle }_{{\rho }_{{\mathsf{th}}}({n}_{0})}}{{n}_{0}^{2}}-\displaystyle \frac{2{\langle N\rangle }_{{\rho }_{{\mathsf{th}}}({n}_{0})}}{{n}_{0}}+1\right]\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{n}_{0}({n}_{0}+1)},\end{eqnarray} \tag{ 43 }$$

where we have applied Gaussian moment factoring (Mandel & Wolf 1995) to obtain Equation (42)—see also Equation (25) of (Helstrom 1968) and Equation (19) of (Pinel et al. 2013) for the final result.

According to Equation (32), the input field is the tensor product of $M={\rm{\Delta }}\nu T$ modes each in the thermal state ${\rho }_{{\mathsf{th}}}({n}_{0})$. From the additivity of the quantum Fisher information for tensor-product states (which follows directly from Equations (33)–(34)), we get the total Fisher information

$$\begin{eqnarray}&&{I}_{Q}^{{\rm{total}}}({n}_{0})=\displaystyle \frac{{\rm{\Delta }}\nu T}{{n}_{0}({n}_{0}+1)},\end{eqnarray} \tag{ 44 }$$

leading to the sought q-CR bound

$$\begin{eqnarray}&&{\rm{Var}}\;{\hat{n}}_{0}\geqslant \displaystyle \frac{{n}_{0}({n}_{0}+1)}{{\rm{\Delta }}\nu T}\end{eqnarray} \tag{ 45 }$$

valid for any unbiased estimator ${\hat{n}}_{0}$ for n₀. The relative sensitivity of any unbiased estimator ${\hat{n}}_{0}$ therefore satisfies

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{Var}}\;{\hat{n}}_{0}}{{n}_{0}^{2}}\geqslant \displaystyle \frac{{n}_{0}+1}{{n}_{0}{\rm{\Delta }}\nu T}.\end{eqnarray} \tag{ 46 }$$

This lower limit on the relative sensitivity of any unbiased estimator of n₀ is our main result. From Equation (12), we see that the relative sensitivity of an unbiased estimator ${\hat{T}}_{{\rm{s}}}$ of the source temperature similarly obeys the limit

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{Var}}\;{\hat{T}}_{{\rm{s}}}}{{T}_{{\rm{s}}}^{2}}\geqslant \displaystyle \frac{1+\displaystyle \frac{h{\nu }_{0}}{{{kT}}_{{\rm{s}}}}}{{\rm{\Delta }}\nu T}.\end{eqnarray} \tag{ 47 }$$