Characterizing photon number statistics using conjugate optical homodyne detection

TL;DR

This paper investigates how to determine photon number statistics of unknown quantum states using conjugate optical homodyne detection, addressing quantum noise limitations and potential quantum communication applications.

Contribution

It introduces a method to recover photon number statistics from conjugate quadrature measurements without phase scanning or randomization, enhancing quantum state characterization.

Findings

Photon number statistics can be reconstructed from conjugate quadrature data.

The method works without phase scanning or phase randomization.

Potential applications in quantum key distribution are discussed.

Abstract

We study the problem of determining the photon number statistics of an unknown quantum state by simultaneously measuring conjugate quadratures with double homodyne detectors. Classically, the sum of the squared outputs of the two homodyne detectors is proportional to the intensity (thus the photon number) of the input light. Quantum mechanically, due to vacuum noise, the above photon number measurement is intrinsically noisy. We quantify the information gain in a single-shot measurement and discuss potential applications of this technology in quantum key distribution. We also show that the photon number statistics can be recovered in repeated measurements on an ensemble of identical input states without scanning the phase of the input state or randomizing the phase of the local oscillator used in homodyne detection.

Equations

$$Z=X_3^2+P_4^2$$

$${\hat{X}}_3=\frac{1}{\sqrt{2}}\left[{\hat{a}}_3^{†}\exp (i\theta )+{\hat{a}}_3\exp (-i\theta )\right],$$

$${\hat{P}}_4=\frac{i}{\sqrt{2}}\left[{\hat{a}}_4^{†}\exp (i\theta )-{\hat{a}}_4\exp (-i\theta )\right].$$

$$\hat{Z}={\hat{X}}_3^2+{\hat{P}}_4^2.$$

$${\hat{a}}_3=\frac{1}{\sqrt{2}}({\hat{a}}_1+{\hat{a}}_2),$$

$${\hat{a}}_4=\frac{1}{\sqrt{2}}({\hat{a}}_1-{\hat{a}}_2),$$

$$[{\hat{a}}_j,{\hat{a}}_j^{†}]=1,j\in \left\{1,2,3,4\right\}$$

$$\hat{Z}={\hat{n}}_1+{\hat{n}}_2+{\hat{a}}_1^{†}{\hat{a}}_2^{†}{e}^{i2\theta }+{\hat{a}}_1{\hat{a}}_2{e}^{-i2\theta }+1$$

$$⟨\hat{Z}⟩=⟨{\psi }_1{0}_2|\hat{Z}|{\psi }_1{0}_2⟩=⟨n_1⟩+1,$$

$$⟨\Delta Z^2⟩=⟨(\hat{Z}-⟨\hat{Z}⟩{)}^2⟩=⟨{\hat{Z}}^2⟩-⟨\hat{Z}{⟩}^2.$$

$$⟨{\hat{Z}}^2⟩=⟨{\hat{n}}_1^2⟩+3⟨{\hat{n}}_1⟩+2.$$

$$⟨\Delta Z^2⟩=⟨\Delta n_1^2⟩+⟨n_1⟩+1,$$

$$g^{(2)}(0)=\frac{⟨{\hat{a}}_1^{†}{\hat{a}}_1^{†}{\hat{a}}_1{\hat{a}}_1⟩}{⟨{\hat{a}}_1^{†}{\hat{a}}_1{⟩}^2}.$$

$$⟨{\hat{a}}_1^{†}{\hat{a}}_1^{†}{\hat{a}}_1{\hat{a}}_1⟩=⟨{\hat{Z}}^2⟩-4⟨\hat{Z}⟩+2.$$

$$g^{(2)}(0)=\frac{⟨{\hat{Z}}^2⟩-4⟨\hat{Z}⟩+2}{(⟨\hat{Z}⟩-1{)}^2}.$$

$$\begin{gathered} P_{X_3,P_4}(x_3,p_4)=\frac{1}{\pi }\sum _{m,n=0}^{\infty }{\rho }_{mn}\frac{\exp \left[i(n-m)\theta \right]}{(m!n!{)}^{1/2}} \\ \times (x_3-i{p}_4{)}^m(x_3+i{p}_4{)}^n\exp \left[-(x_3^2+p_4^2)\right]. \end{gathered}$$

$$P_R(r)={\int }_0^{2\pi }{P}_{X_3,P_4}(r\cos \varphi ,r\sin \varphi )d\varphi .$$

$$P_R(r)=2\exp (-r^2)\sum _{n=0}^{\infty }\frac{{\rho }_{nn}}{n!}{r}^{2n+1}.$$

$$P_Z(z)=\exp (-z)\sum _{n=0}^{\infty }\frac{{\rho }_{nn}}{n!}{z}^n.$$

$$P(Z=z|n)=\exp (-z)\frac{z^n}{n!}.$$

$$P(N=n|z)=\frac{P(Z=z|n)P_N(n)}{P_Z(z)}=\exp (-z)\frac{z^n}{n!}.$$

$$\sigma =⟨\Delta n^2⟩=z.$$

$$P(N\ge n_{th}|z)=\sum _{n=n_{th}}^{\infty }\exp (-z)\frac{z^n}{n!}.$$

$$\begin{array}{r}P(N\ge n_{th}|z)\le \sum _{n=n_{th}}^{\infty }\exp (-z)\frac{z^n}{\sqrt{2\pi }{n}^{n+\frac{1}{2}}\exp (-n)}\\ <\sum _{n=n_{th}}^{\infty }\exp (-z)\frac{z^n}{\sqrt{2\pi }{n}_{th}^n\exp (-n)}\\ =\frac{\exp (-z)}{\sqrt{2\pi }}\sum _{n=n_{th}}^{\infty }{\left(\frac{ze}{n_{th}}\right)}^n=\frac{\exp (-z)}{\sqrt{2\pi }}\frac{a^{n_{th}}}{1-a}\end{array}$$

$$\eta ={\int }_T^{\infty }{P}_Z(z|1)dz$$

$$D={\int }_T^{\infty }{P}_Z(z|0)dz$$

$$p(Z_k=z_k|\{p(n)\},n_k)=p(n_k)\frac{\exp (-z_k)z_k^{n_k}}{n_k!}.$$

$${\mathcal{L}}_c=\prod _{k=1}^{M}p(Z_k=z_k|\{p(n)\},n_k).$$

$$\mathcal{L}=\prod _{k=1}^M\sum _{n_k=0}^{N_{max}}p(Z_k=z_k|\{p(n)\},n_k),$$

$$\prod _{k=1}^M\sum _{n_k=0}^{N_{max}}=\sum _{n_1=0}^{N_{max}}\cdots \sum _{n_M=0}^{N_{max}}.$$

Full text

I Introduction

A single photon detector (SPD) is a workhorse of modern quantum optics experiments. It is commonly used to determine the number of photons in a given light pulse in a single-shot measurement scenario Hadfield09 ; Eisaman11 . While photon-number-resolving detection schemes have been demonstrated using either a singls SPD or SPD arrays Rosenberg05 ; Waks03 ; Cahall17 ; Banaszek03 ; Fitch03 ; Zambra05 , most commercial SPDs are threshold detectors which can only distinguish between vacuum and non-vacuum states. Alternatively, given a large ensemble of identical input states, optical homodyne detection (OHD) can be used to completely reconstruct the quantum state of the ensemble, including the photon number statistics. This method is known as the optical homodyne tomography Vogel89 ; Smithey93 ; Lvovsky09 .

In OHD, a strong local oscillator (LO) is mixed with a weak input signal at a beam splitter. The interference signals can be strong enough to be detected by low cost, highly efficient photo-diodes working at room temperature. This makes OHD an appealing solution in practice, for example, in chip-size implementation Raffaelli18 . Note that the LO in OHD also functions as a mode selector: only the signals in the same spatiotemporal mode as the LO will be detected. This can be advantageous in certain applications such as quantum key distribution (QKD) Gisin02 ; Scarani09 ; Lo14 ; Diamanti16 , where the mode selecting function can effectively suppress broadband background noise originating from the communication channel Qi10 ; Heim14 ; Kumar15 ; Tobias19 .

Given the LO is sufficiently strong, a DC-balanced homodyne detector measures quadrature $X_{\theta}$ of the input signal, where $\theta$ is the phase of the LO Yuen83 ; Abbas83 . To reconstruct the full quantum state, repeated measurements are required for all values of $\theta\in[0,2\pi]$. Obviously, a crucial requirement for such a reconstruction scheme is that the phase of the input state is well controlled. However, in many applications, this requirement can not be easily satisfied. For example, in QKD, the quantum signals detected by the receiver come from a channel controlled by an adversary. In this case, we cannot make any assumptions about the phase of incoming signals. Fortunately, it has been shown that the photon number statistics of an input state with an unknown phase can still be fully recovered by using either one or two optical homodyne detectors, given the phase of the LO is uniformly randomized Munroe95 ; Schiller96 ; Leonhardt96 ; Richter98 ; Chrzanowski13 . In the case of QKD, phase randomization can be implemented by using a phase modulator driven by a random pattern, as demonstrated in Zhao07 . Nevertheless, the requirements of truly random numbers and high-speed modulator introduce additional complexities.

In this paper, we show that the conjugate optical homodyne detection scheme Richter98 ; Walker86 ; Noh91 can be used to determine the photon number statistics without controlling the phase of the input quantum state or randomizing the phase of the LO. In this scheme, two homodyne detectors are employed to measure conjugate quadratures of the input state simultaneously. We define a measurement observable $Z$ as the sum of the squared outputs of the two homodyne detectors. In classical electrodynamics, the outputs of two conjugate homodyne detectors correspond to the in-phase and out-of-phase components of an electromagnetic wave, so the observable $Z$ defined above is proportional to the intensity (or the photon number) of the input signal. In quantum mechanics, canonically conjugate quadrature components of quantum optical fields do not commute with each other and thus cannot be determined simultaneously and noiselessly due to Heisenberg’s uncertainty principle. So, a single-shot measurement of $Z$ is intrinsically noisy. Nevertheless, by repeating the $Z$ measurement on a large ensemble of identical states, the photon number statistics can still be determined, as we will show in this paper. We will also study the problem of quantifying the photon number statistics of an unknown input state based on a single-shot measurement. This problem is seldom discussed in previous studies. Nevertheless, the solutions developed here may find important applications in QKD and other areas.

This paper is organized as follows: in Sec.II, we present the theory of conjugate homodyne detection. In Sec.III, we study the case of single-shot measurement and discuss its potential applications in QKD. In Sec.IV, we study the case of repeated measurements on an ensemble of identical input states. Finally, we conclude this paper with a brief summary in Sec.V.

II Conjugate homodyne detection

The basic setup of a conjugate homodyne detection system is shown in Fig.1. The input quantum state $|\psi\rangle$ is split into two parts by a symmetric beam splitter ($\textrm{BS}_{1}$ in Fig.1). One output of the beam splitter (mode 3 in Fig.1) is measured by an optical homodyne detector with an LO phase $\theta$; the other output state (mode 4 in Fig.1) is measured by another optical homodyne detector with an LO phase $\theta+\pi/2$. The two LOs can be generated from a common laser using a beam splitter and a $\pi/2$ phase shifter. The common phase $\theta$ is defined using the phase of the input state as a reference. Provided that the input state has an unknown phase which may change from pulse to pulse, we have no control of $\theta$ (i.e. $\theta$ is a random variable with an unknown distribution). We remark that the setup shown in Fig.1 (plus the beam splitter and the $\pi/2$ phase shifter for generating two LOs from a single laser) can be conveniently implemented with a compact commercial $90^{0}$ optical hybrid Optoplex . In this Section, we assume noiseless homodyne detectors with unity efficiency. We will discuss the case of non-unity detection efficiency in Sec.IV.

Given the LOs are sufficiently strong, the outputs of the two homodyne detectors are quadrature components of mode 3 ($X_{3,\theta}$) and mode 4 ($X_{4,\theta+\pi/2}$). For simplicity, we use $X_{3}$ and $X_{4}$ to represent $X_{3,\theta}$ and $X_{4,\theta+\pi/2}$ in the rest of the paper.

We define a new parameter

[TABLE]

as an estimation of the photon number of the input state $|\psi\rangle$. Intuitively, $Z$ is proportional to the intensity of the input light in classical electrodynamics.

In quantum optics, the two homodyne outputs are represented by operators $\hat{X}_{3}$ and $\hat{P}_{4}$, which are defined in terms of photon annihilation operator $\hat{a}$ and photon creation operator $\hat{a}^{\dagger}$ as

[TABLE]

[TABLE]

We define an operator $\hat{Z}$ as

[TABLE]

Using the transformation relations of a symmetric beam splitter Loudon00

[TABLE]

[TABLE]

and the commutation relation

[TABLE]

it can be shown that

[TABLE]

where $\hat{n}_{1}=\hat{a}_{1}^{\dagger}\hat{a}_{1}$ and $\hat{n}_{2}=\hat{a}_{2}^{\dagger}\hat{a}_{2}$ are photon number operators of mode 1 and 2. Obviously, $\hat{Z}$ is a Hermitian operator. Thus, it is a valid observable.

II.1 Expectation value of $\hat{Z}$

As shown in Fig.1, the joint input state of mode 1 and 2 is given by $|\psi_{1}0_{2}\rangle$. From (8), the expectation value of $\hat{Z}$ can be determined to be

[TABLE]

where $\langle n_{1}\rangle$ is the average photon number of the input quantum state. The constant $1$ on the RHS of (9) can be interpreted as vacuum noise contribution Note1 . Eq.(9) shows that the average photon number of the input state can be estimated by subtracting $1$ from the expectation value of $\hat{Z}$.

II.2 Variance of $\hat{Z}$

The variance of $\hat{Z}$ is given by

[TABLE]

From (8), it is straightforward to show

[TABLE]

Using (9)-(11), we obtain

[TABLE]

where $\langle\Delta n_{1}^{2}\rangle=\langle\hat{n}_{1}^{2}\rangle-\langle n_{1}\rangle^{2}$ is the variance of photon number distribution of the input state. Clearly, the additional noise of the proposed scheme is $\langle n_{1}\rangle+1$.

II.3 Second-order correlation function $g^{(2)}(0)$

The single-time second-order correlation function $g^{(2)}(0)$ is an important parameter for characterizing a photon source Glauber63 . In the case of a single homodyne detector with a phase randomized LO, it has been shown that $g^{(2)}(0)$ can be determined from the measurement statistics Roumpos13 ; Luders18 . Here, we show that $g^{(2)}(0)$ can also be determined from the statistics of $Z$ measurement.

As defined in Glauber63

[TABLE]

Using (7), the numerator on the RHS of (13) can be written as $\langle\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}\hat{a}_{1}\hat{a}_{1}\rangle=\langle\hat{a}_{1}^{\dagger}\hat{a}_{1}\hat{a}_{1}^{\dagger}\hat{a}_{1}-\hat{a}_{1}^{\dagger}\hat{a}_{1}\rangle=\langle\hat{n}_{1}^{2}\rangle-\langle\hat{n}_{1}\rangle$. Using Eq.(9) and Eq.(11), we have

[TABLE]

From (9), the denominator on the RHS of (13) is simply $(\langle\hat{Z}\rangle-1)^{2}$. Finally, we have

[TABLE]

II.4 Probability density function of $\hat{Z}$

Unlike photon number $n$, the parameter $Z$ measured in our scheme is a continuous variable. Given an arbitrary input state described by the density matrix $\rho$, we would like to determine the probability density function (PDF) of $Z$.

The joint PDF of quadrature components $X_{3}$ and $P_{4}$ of conjugate homodyne detection has been derived in Richter98 as

[TABLE]

The data pair $\left(x_{3},p_{4}\right)$ can be interpreted as the Cartesian coordinates of a point, which relates to the Polar Coordinates $\left(r,\phi\right)$ by $x_{3}=r\cos\phi$ and $p_{4}=r\sin\phi$. The marginal distribution of $r$ is given by

[TABLE]

Note the term $(x_{3}-ip_{4})^{m}(x_{3}+ip_{4})^{n}$ on the RHS of (16) is transformed into $r^{2}\exp[i(n-m)\phi]$ in the Polar Coordinates. Obviously, only terms with $n=m$ have non-zero contribution in (17). It is straightforward to show

[TABLE]

Since $Z=X_{3}^{2}+P_{4}^{2}=R^{2}$, the PDF $P_{Z}(z)$ can be determined from (18) as

[TABLE]

Eq.(19) shows the relation between the PDF of $Z$ and the photon number distribution $\rho_{nn}$ of the input state. Note that $P_{Z}(z)$ is only dependent on the diagonal terms of the input state, a feature we would expect from a “phase-insensitive” photon detector. Experimentally, Eq.(19) implies that the photon number statistics can be determined without scanning (or randomizing) the phase of the LO.

III Single-shot measurement

In most of previous studies on OHD, one assumption is that an ensemble of identical input states are available. This allows precise quantum state characterization based on repeated measurements. However in some applications, we may need to make an estimation of photon numbers based on a single-shot measurement. One example is continuous-variable (CV) QKD Ralph99 ; Hillery00 ; GMCSQKD , where the transmitter (Alice) encodes random bits on quadratures of either squeezed states or weak coherent states and the receiver (Bob) measures either one or two quadratures using optical homodyne detection. In these protocols, it is crucial to upper bound the photon number of the quantum state received by Bob for various reasons. First, as shown in Qin16 ; Qin18 , an adversary may attack Bob’s QKD system by sending strong laser pulses and forcing the homodyne detectors to work in nonlinear or saturation region. It thus important to experimentally verify that the photon number of the incoming signal is within the normal range. Second, in the recent security proofs of CV-QKD with discrete modulations Ghorai19 ; Lin19 ; Kaur19 , numerical approaches are employed to determine lower bounds of secure key rates. To reduce the dimension of the relevant quantum space from infinite to finite, one important assumption is that the photon number received by Bob is finite. The scheme proposed in this paper provides a practical way to upper bound the incoming photon number. More specifically, as we will show below, given a single measurement of $Z$, it is possible to determine a threshold value of photon number $n_{th}$, such that the probability that the received signal containing $n_{th}$ or more photons is negligible.

Given the input state is a Fock state $\rho=|n\rangle\langle n|$, the likelihood of a measurement output of $z$ can be determined from (19) as

[TABLE]

Using Bayes’ rule the likelihood of $n$ photons coming in given the measurement output $z$ is

[TABLE]

In the lase step of (21), we have assumed that the prior $P_{N}(n)$ is a uniform distribution. This leads to a uniform PDF of $P_{Z}(z)$ from (19).

The distribution in (21) is the Poisson distribution, which quantifies the uncertainty of the photon number $n$ given a single measurement of $Z$. The uncertainty of the photon number can be determined from (21) as

[TABLE]

Using (21), given a single $Z$ measurement, the probability that the received signal containing $n_{th}$ or more photons is

[TABLE]

where $n_{th}$ is assumed to be much larger than $z$.

Using Stirling’s formula $n!\geq\sqrt{2\pi}\times n^{n+\frac{1}{2}}\times\exp(-n)$, we can derive the following inequality from (23)

[TABLE]

where $a=\dfrac{ze}{n_{th}}$ is assumed to be much less than 1.

From (24), by choosing an appropriate $n_{th}$, $P(N\geq n_{th}|z)$ can be negligible. In practical QKD, the average photon number of the quantum state received by Bob is typically much less than one, so $z$ is most likely a small number. In this case, a moderate $n_{th}$ can serve as the upper bound of the incoming photon number. For example, if $z=3$, the probability that the received signal containing $n_{th}=30$ or more photons is less than $3\times 10^{-19}$.

We remark that the proposed scheme can be easily implemented in CV QKD based on conjugate homodyne detection (also called “heterodyne detection” in literature) Weedbrook04 ; Qi15 . In this case, Bob measures both X and P quadratures simultaneously for key generation. So the $Z$ information is automatically available to Bob without any changes to the QKD system or the measurement procedures.

In principle, the proposed system may also act as a threshold single-photon detector (SPD) which can discriminate vacuum state from non-vacuum states but cannot resolve photon number. This kind of detector is commonly using in discrete-variable (DV) QKD. The basic strategy is to choose an optimal threshold value $T\in\left[0,\infty\right)$ for the $Z$ measurement: if the measurement result $z$ is smaller (larger) than $T$, the detected quantum state is assumed to be vacuum (non-vacuum) state.

Two important parameters of a threshold SPD are detection efficiency and dark count probability. The detection efficiency $\eta$ is defined as the conditional probability that the detector reports a non-vacuum state given the input is a single photon Fock state. The dark count probability $D$ is defined as the conditional probability that the detector reports a non-vacuum state given the input is vacuum. From (20), these two parameters can be determined by

[TABLE]

[TABLE]

Fig.2 shows the simulation results of $\eta$ and $D$ as a function of the threshold value $T$. By choosing an appropriate threshold value, we can either achieve a high detection efficiency or a low dark count. Unfortunately, we cannot have both at the same time. Similar conclusions had been drawn in previous studies based on the single homodyne detection scheme Grice96 . In Fig.2, we also present the ratio $R=\eta/D$, which is an important figure of merit in applications like QKD. A state-of-the-art SPD can provide a R-value as high as $10^{8}$ Hadfield09 . In comparison, the R-vale of the proposed scheme is less than 10 in the region with a reasonable detection efficiency.

In next section, we discuss the case of determining photon number distribution from repeated measurements.

IV Repeated measurements

Given a sequence of M independent measurements of $Z$, denoted as $\{z_{i}\}=z_{1},\dots,z_{M}$, we would like to infer the underlying photon number statistics given by the distribution $p(n)=\rho_{nn}$ for $n=0,1,\dots,N_{max}$. The observed data statistics can be described in terms of a mixture model of $N_{max}+1$ gamma distributions ${\rm Gamma}(n+1,1)$ with mixing coefficients given by the photon number distribution $p(n)$. In this model the conditional probability of observing an outcome $Z_{k}=z_{k}$ reads,

[TABLE]

The complete data likelihood function $\mathcal{L}_{c}$ for the entire measurement sequence can be written as,

[TABLE]

Here, $p(n),n=0,\cdots,N_{max}$ are the unknown parameters (Fock state probabilities) we wish to infer and $n_{1},\dots,n_{M}$ are random variables, each in range $[0,N_{max}]$, that determine which mixture component $n_{i}\in[0,N_{max}]$ has generated an outcome $z_{i}$. Have we known the values of $n_{1},\dots,n_{M}$ then we could use the maximum-likelihood estimation (MLE) based on $\mathcal{L}_{c}$ to infer the most likely values of the parameters $p(n)$. Unfortunately the variables $n_{i},i={1,M}$ are unobserved (latent). Therefore, in order to use MLE we first need to marginalize the complete data likelihood $\mathcal{L}_{c}$ over the latent variables. The resulting marginal likelihood of the observed data reads,

[TABLE]

where we used the following convention,

[TABLE]

In principle one now can try and maximize the function $\mathcal{L}$ with respect to parameters $p(n)$, but in practice the complexity of this task will grow exponentially with the number of measurements $M$. Fortunately, there is a way to determine the maximum of $\log\mathcal{L}$ that avoids explicit maximization of $\mathcal{L}$. This method is widely used for inference in mixture models and is called expectation-maximization (EM) algorithm EM . Applicability of EM in the context of homodyne measurement of photon statistics has been advocated previously Banaszek98 . EM is an iterative procedure that uses the expected value of the complete data likelihood $\mathcal{L}_{c}$ with respect to the latent random variables $n_{1},\dots,n_{M}$ as a maximization objective function for determining $p^{t}(n)$ – the $t$-th iteration estimate of the parameters $p(n)$. Here is how it works:

 

First, set initial estimates of the parameters $p^{0}(n)$ to some values. We chose a uniform prior

[TABLE]

as it is uninformative and easy to implement. Due to the multimodality of Eq. (29) the choice of prior will bias the EM reconstruction given below. It is an open problem to identify the best prior for this application111In principal, one can choose any set of positive real numbers such that $\sum\limits_{n=0}^{N_{max}}p^{0}(n)=1$.. Next, repeat the following steps until convergence criteria are satisfied.

• •

At the $t$-th iteration, update probabilities $p(n_{k}|z_{k},\{p^{t}(n)\})$ for all latent variables $n_{1},\dots,n_{M}$ using Bayes rule with $p^{t}(n)$ as a prior,

    $p(n_{k}|z_{k},\{p^{t}(n)\})=\frac{p(Z_{k}=z_{k}|\{p^{t}(n)\},n_{k})}{\sum\limits_{n_{k}=0}^{N_{max}}p(Z_{k}=z_{k}|\{p^{t}(n)\},n_{k})}$

• •

Calculate the expected value of the complete data $\log$ likelihood with respect to the updated distribution $p(n_{k}|z_{k},\{p^{t}(n)\})$,

    $Q(\{p(n)\}|\{p^{t}(n)\})={\rm E}_{\{n_{k}\}|\{z_{k}\},\{p^{t}(n)\}}[\log\mathcal{L}_{c}]=$

$\sum\limits_{k=1}^{M}\sum\limits_{j=0}^{N_{max}}p(j|z_{k},\{p^{t}(n)\})\log[p(Z_{k}=z_{k}|\{p(n)\},j)]$

• •

Find parameter values $\{p_{max}(n)\}$ that maximize $Q(\{p(n)\}|\{p^{t}(n)\})$ and set the next iteration estimates of photon number probabilities $\{p^{t+1}(n)\}=\{p_{max}(n)\}$. Note that the values of $\{p_{max}(n)\}$ can be calculated analytically,

          $p^{t+1}(n)=\frac{1}{M}\sum\limits_{j=1}^{M}p(n|z_{j},\{p^{t}(n)\})$.

It can be shown EM that iterative maximization of $Q(\{p(n)\}|\{p^{t}(n)\})$ also results in the maximization of the marginal likelihood $\mathcal{L}$ in Eq.(29). Therefore, after a sufficient number of iterations $t_{s}$ our MLE estimator of the photon number statistics is given by the distribution $\{p^{t_{s}}(n)\}$.

IV.1 Simulation results

To illustrate our EM-based photon number statistics inference scheme we applied it to a sequence of $32768$ simulated homodyne measurement outcomes for a coherent state $\rho=|\alpha\rangle\langle\alpha|$ with the mean photon number $|\alpha|^{2}=5$. To reconstruct the photon number statistics from the simulated measurement data, we have selected a mixture model with $N_{max}=20$ ($21$-component model). The results of EM reconstruction (after 9 iterations) are plotted on Fig.3 and demonstrate a good quantitative agreement with the true state.

IV.2 Experimental results

We apply the theory described above in experiments to reconstruct photon number distributions of a weak coherent state and a thermal state.

In Sec.II, we have assumed perfect homodyne detectors with unity efficiency. Here we consider practical detector with non-unity detection efficiency. It is well known that a realistic photo-detector with efficiency $\eta$ can be modeled by placing a virtual beam splitter (with a transmittance of $\eta$) in front of an ideal detector Yuen78 . By assuming the four photo-detectors have identical efficiency, we can model the conjugate homodyne detection using the setup shown in Fig.4(a). In Appendix A we will show that given the LOs are strong enough, the setup in Fig.4(a) is equivalent to that in Fig.4(b), where the four virtual beam splitters in front of the photo-detectors are replaced by a common virtual beam splitter (with the same transmittance) at the input of the first beam splitter. This is convenient in practice since we can apply the theory in Sect. II directly to the experimental results by assuming the photo-detectors are ideal. The photon number distribution reconstructed this way is related to that of the input state by the Bernoulli transformation Campos89 . By further applying the inverse Bernoulli transformation, the photon number distribution of the input state can be determined.

In our experiments, either a weak coherent source or a thermal source is employed as the input. The state after the virtual beam splitter in Fig.4(b) is still a coherent state (or a thermal state) with a reduced average photon number. This is convenient since we can simply redefine the state after the virtual beam splitter as the input state and assume the detection efficiency is one.

The experimental setup is shown in Fig.5. A 1550nm continuous wave (CW) laser is employed as the LO. The conjugate optical homodyne detection system is constructed by a commercial $90^{o}$ optical hybrid (Optoplex) and two $350$ MHz balanced amplified photodetectors (Thorlabs). Variable optical attenuators are used to balance the detection efficiency of different channels and control the average photon number of the input state. The outputs of the two balanced photodetectors are sampled by a two-channel data acquisition board (Texas Instruments). The overall efficiency of the detection system is 0.5, with electrical noises variances (in shot-noise unit) of $\sigma^{2}_{X_{3}}=0.21,\sigma^{2}_{P_{4}}=0.16$ for the quadratures $X_{3}$ and $P_{4}$ respectively.

In the first experiment, a heavily attenuated laser source is used to provide a weak coherent state input. By adjusting the variable optical attenuator, the average photon number within one sampling window has been set to 5 (after correction of detection efficiency). In the second experiment, an amplified spontaneous emission (ASE) source is used to provide a thermal state input (with an average photon number of 15.3). Note, while the output of the ASE source contains multiple modes, the homodyne detector selectively measures the one matched with the LO mode. Limited by the memory size of the data acquisition board, 32728 quadrature pairs are sampled in each measurement.

Using (15), the $g^{(2)}(0)$ factors of the two sources have been determined to be $1.11\pm 0.02$ (weak coherent source) and $1.94\pm 0.02$ (thermal source) correspondingly, where the uncertainty quantifies statistic fluctuation. These results match with the theoretical values of 1 (ideal coherent state) and 2 (ideal thermal state) reasonably well.

We apply the EM algorithm described above to the experimental data in order to reconstruct the photon number statistics of the two light sources. In Fig.6 we plot reconstruction results for a weak coherent state. The blue bars represent the photon number statistic of a coherent state $|\alpha\rangle$ with the mean number of photons $|\alpha|^{2}=5$, assuming noiseless detectors. The red bars correspond to the numerically synthesized data from the coherent state $|\alpha\rangle$ with added Gaussian noise which mimics the effect of noisy photo detectors. The green bars depict the photon number statistics obtained from raw experimental data by using the EM algorithm. We note that the red and green histograms look remarkably similar which implies that experimental data come from a coherent state affected by the detector noise. In Fig.7 we depict EM reconstruction results for a thermal state with a mean photon number of 15.3. The black line represents the theoretical distribution of an ideal thermal state. The red bar are EM reconstruction results using numerically simulated quadrature measurement data with added detector noise. Finally, the blue bars correspond to EM reconstruction of actual experimental data. We notice that for the thermal state the three distributions are very close visually. This is because the detector noise is much smaller than the mean number of photons in this case. Therefore, the noise effects are not as pronounced as in the case on a weak coherent state.

V Summary

Classically, the intensity of a single-mode light pulse can be determined by measuring two conjugate quadratures simultaneously. Quantum mechanically, the above measurement process is intrinsically noisy. In this paper, we develop theoretical tools to reconstruct photon number statistics of a single-mode quantum state by performing conjugate homodyne detection. Comparing with previous studies based on single homodyne detection, no LO phase randomization is required in our scheme. We also show how to determine an upper bound of incoming photon number based on a single-shot measurement. Such a technology could be useful in CV-QKD to prevent bright pulse attack or validate crucial assumptions in security proofs.

We acknowledge helpful comments from Ryan S. Bennink, Warren Grice, Charles C. W. Lim and Nicholas A. Peters. This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). The authors acknowledge support from ORNL laboratory directed research and development program.

Appendix A Detector efficiency in the conjugate homodyne scheme

In this Appendix, we will show that given the LOs are strong enough, the setup shown in Fig.4(a) is equivalent to that in Fig.4(b), where the four virtual beam splitters in front of the photo-detectors are replaced by a common virtual beam splitter (with the same transmittance) at the input path of the first beam splitter.

Given the LO is strong enough, a single DC-balanced homodyne detector using two realistic (non-unity efficiency) and identical photo-detectors can be modeled by the one with ideal photo-detectors by placing a virtual beam splitter at the signal input Leonhardt93 . This allows us to replace the four virtual beam splitters in Fig.4(a) by two virtual beam splitters with the same transmittance (one at each of the output port of the first beam splitter). We will show that these two virtual beam splitters can be further replaced by one as shown in Fig.4(b). More specifically, we will show the two models in Fig.8 are equivalent to each other.

For simplicity, we define the transmittance of the virtual beam splitter as $\eta=\cos\gamma$.

From Fig.8(a) and using the transmission relations of lossless beam splitter, we have

[TABLE]

where $\hat{X}_{V1}=\dfrac{\cos\gamma}{\sqrt{1+\sin^{2}\gamma}}\hat{X}_{2}+\dfrac{\sqrt{2}\sin\gamma}{\sqrt{1+\sin^{2}\gamma}}\hat{X}_{5}$. Since the inputs of mode 2 and mode 5 in Fig.8(a) are vacuum, the unitary transformation described above yields another vacuum state. So $\hat{X}_{V1}$ at the RHS of (A1) can be interpreted as the X-quadrature of vacuum state.

Similarly, the P-quadrature of mode 4 in Fig.8(a) is given by

[TABLE]

where $\hat{P}_{V2}=-\dfrac{\cos\gamma}{\sqrt{1+\sin^{2}\gamma}}\hat{P}_{2}+\dfrac{\sqrt{2}\sin\gamma}{\sqrt{1+\sin^{2}\gamma}}\hat{P}_{6}$.

We can apply the same process in the model shown in Fig.8(b) and have the following relations:

[TABLE]

where $\hat{X}_{V3}=\dfrac{\sin\gamma}{\sqrt{1+\sin^{2}\gamma}}\hat{X}_{5^{\prime}}+\dfrac{1}{\sqrt{1+\sin^{2}\gamma}}\hat{X}_{2^{\prime}}$

[TABLE]

where $\hat{P}_{V4}=\dfrac{\sin\gamma}{\sqrt{1+\sin^{2}\gamma}}\hat{P}_{5^{\prime}}-\dfrac{1}{\sqrt{1+\sin^{2}\gamma}}\hat{P}_{2^{\prime}}$.

It is easy to show $X_{V1}$, $P_{V2}$, $X_{V3}$, $P_{V4}$ are independent and identically distributed random variables. From (A1)-(A4), the joint probability of $X_{3}$ and $P_{4}$ is the same as that of $X_{3^{\prime}}$ and $P_{4^{\prime}}$. So the two models given in Fig.8 are equivalent.