Study of higher-order correlation functions and photon statistics using multiphoton-subtracted states and quadrature measurements

TL;DR

This paper investigates high-order correlation functions and photon statistics in quantum optics by using multiphoton subtraction and quadrature measurements to improve state preparation and measurement accuracy.

Contribution

It introduces a method to analyze quantum states after photon subtraction, enhancing the precision of quantum state generation and characterization.

Findings

Successful preparation of photon-subtracted states with 1 to 10 photons

Application of maximum likelihood estimation for parameter accuracy

Exploration of photon bunching regimes in gamma-compounded Poisson states

Abstract

The estimation of high order correlation function values is an important problem in the field of quantum computation. We show that the problem can be reduced to preparation and measurement of optical quantum states resulting after annihilation of a set number of quanta from the original beam. We apply this approach to explore various photon bunching regimes in optical states with gamma-compounded Poisson photon number statistics. We prepare and perform measurement of the thermal quantum state as well as states produced by subtracting one to ten photons from it. Maximum likelihood estimation is employed for parameter estimation. The goal of this research is the development of highly accurate procedures for generation and quality control of optical quantum states.

Tables

Table 1. Table 1: Results of statistical reconstruction from experimental data

State	$\mu \pm {\sigma }_{\mu }$	$a\pm {\sigma }_a$	Sample size	Fidelity	${\chi }^2$ significance level,
Thermal $m=0$	3.034±0.022	0.999±0.016	50000	99.999	0.370
$m=1$	5.983±0.051	1.605±0.036	25000	99.659	0.112
$m=2$	9.063±0.093	2.515±0.088	12500	99.777	0.550
$m=3$	12.261±0.150	3.149±0.147	7500	99.572	0.014
$m=4$	15.538±0.223	4.331±0.281	4500	99.799	0.565
$m=5$	17.957±0.244	5.198±0.350	4500	99.831	0.111
$m=6$	21.050±0.362	6.378±0.600	2500	99.928	0.214
$m=7$	24.732±0.410	7.045±0.665	2500	99.832	0.021
$m=8$	27.795±0.433	8.847±0.882	2500	99.945	0.207
$m=9$	30.536±0.999	11.259±2.680	500	99.896	0.021
$m=10$	33.114±1.267	11.337±3.156	358	99.980	0.123

Equations

$$P\left(k\right)=\frac{G^{\left(k\right)}\left(0\right)}{k!};$$

$$E\left[k\left(k-1\right)...\left(k-m+1\right)\right]=G^{\left(m\right)}\left(1\right),$$

$$g^{\left(m\right)}=\frac{G^{\left(m\right)}\left(1\right)}{{\mu }^m},m=1,2,...$$

$${\rho }^{out}=Na{\rho }^{in}{a}^{†}.$$

$$G_1\left(z\right)=\frac{G^{\left(1\right)}\left(z\right)}{G^{\left(1\right)}\left(1\right)}=\frac{G^{\left(1\right)}\left(z\right)}{\mu }.$$

$$G_1\left(z\right)=\frac{G^{\left(1\right)}\left(z\left(1-p\right)\right)}{G^{\left(1\right)}\left(1-p\right)}.$$

$$G_m\left(z\right)=\frac{G^{\left(m\right)}\left(z\right)}{\mu {\mu }_1\cdot \cdot \cdot {\mu }_{m-1}},m=1,2,...$$

$$g^{\left(m\right)}=\frac{{\mu }_1{\mu }_2\cdot \cdot \cdot {\mu }_{m-1}}{{\mu }^{m-1}},m=2,3,...$$

$$g^{\left(2\right)}=\frac{{\mu }_1}{\mu }.$$

$$g^{\left(m+1\right)}=g^{\left(m\right)}\frac{{\mu }_m}{\mu },m=1,2,...$$

$$g_m^{\left(n\right)}=\frac{G_m^{\left(n\right)}\left(1\right)}{{\mu }_m^n}=\frac{g^{\left(m+n\right)}{\left(g^{\left(m\right)}\right)}^{n-1}}{{\left(g^{\left(m+1\right)}\right)}^n},n=1,2,...;m=1,2,...$$

$$G_0\left(z|\lambda \right)=\exp \left(-\lambda \left(1-z\right)\right).$$

$$P\left(\lambda \right)=\frac{b^a{\lambda }^{a-1}{e}^{-b\lambda }}{\Gamma \left(a\right)},$$

$$G\left(z|a,b\right)=\underset{0}{\overset{\infty }{\int }}{G}_0\left(z|\lambda \right)P\left(\lambda \right)d\lambda =\frac{1}{{\left(1+\frac{\left(1-z\right)}{b}\right)}^a}.$$

$$G\left(z|\mu ,a\right)=\frac{1}{{\left(1+\frac{\mu \left(1-z\right)}{a}\right)}^a}.$$

$$g^{\left(m\right)}=\frac{{\left(a\right)}_m}{a^m}=\frac{a\left(a+1\right)...\left(a+m-1\right)}{a^m}.$$

$${\left(a\right)}_k=a\left(a+1\right)...\left(a+k-1\right).$$

$$a^{\left(1\right)}=a+1.$$

$${\mu }^{\left(1\right)}=\frac{\mu \left(a+1\right)}{a}.$$

$$a^{\left(k\right)}=a+k;$$

$${\mu }^{\left(k\right)}=\frac{\mu \left(a+k\right)}{a}.$$

$$a^{\left(k\right)}=1+k,{\mu }^{\left(k\right)}=\mu \left(1+k\right).$$

$$\rho \left(x,x^{\prime }\right)=\sum _{k=0}^{\infty }P\left(k|\mu ,a\right){\phi }_k\left(x\right){\phi }_k^{\ast }\left(x^{\prime }\right)$$

$${\phi }_k\left(x\right)=\frac{1}{{\left(2^{k}k!\sqrt{\pi }\right)}^{1/2}}{H}_k\left(x\right)\exp \left(-\frac{x^2}{2}\right),k=0,1,2,...$$

$$P\left(x|\mu ,a\right)=\sum _{k=0}^{\infty }P\left(k|\mu ,a\right){\left|{\phi }_k\left(x\right)\right|}^2.$$

$${\sigma }^2=\mu +\frac{1}{2}.$$

$${\beta }_2=-6{\left(\frac{\mu }{2\mu +1}\right)}^2\frac{\left(a-1\right)}{a}.$$

$${\beta }_1=\frac{E\left[{\left(x-E\left(x\right)\right)}^3\right]}{{\sigma }^3},$$

$${\beta }_2=\frac{E\left[{\left(x-E\left(x\right)\right)}^4\right]}{{\sigma }^4}-3.$$

$$\mu ={\sigma }^2-\frac{1}{2},$$

Full text

\supit

aInstitute of Physics and Technology, Russian Academy of Sciences, 117218, Moscow, Russia;

\supitbNational Research Nuclear University “MEPhI”, 115409, Moscow, Russia;

\supitcNational Research University of Electronic Technology MIET, 124498, Moscow, Russia;

\supitdM.V. Lomonosov Moscow State University, 119991, Moscow, Russia

1 Introduction

Preparation and measurement of optical quantum states are key problems in applied quantum information technologies. Among the currently used states of light, the thermal state plays a special role. It serves as a testbed for various effects based on quantum and classical correlations, while being easy to prepare.

The pioneering work of Brown and Twiss [1], which is considered to be the first quantum optics experiment, explored correlation in thermal light using a beam splitter and a coincidence circuit with two detectors. Since then thermal states have been used in many applications including ghost imaging [2, 3, 4], quantum illumination [5], and “thermal laser” [6]. A recent demonstration of classical teleportation [7] also used thermal states. In this paper we consider the properties of photon number statistics and autocorrelation functions in photon-subtracted thermal states.

Photon addition and subtraction is of the great interest in quantum optics, because it provides a tool for direct tests of basic commutation relations [8] and enables the preparation of Schrodinger cat and other exotic quantum states [9]. It can also be used for probabilistic linear no-noise amplification [10]. One- and two-photon subtracted thermal states were demonstrated for the first time in [11]. The measurement of photon statistics with photon number resolving detectors was demonstrated in [12]. In the present work we provide a comprehensive description of multiphoton subtracted thermal states, based on a general approach, suitable for any photon number distribution. We demonstrate the technique of high-fidelity preparation and reconstruction of up to 10-photon subtracted thermal states, using detectors not capable of resolving the photon numbers. The technique discussed in the paper can also be used in some metrological applications [13, 14].

2 Probability generating functions and autocorrelation

The approach is based on exploiting the properties of generating functions. Probability generating functions $G\left(z\right)$ contain all information about the random distribution of photons. In particular, the probability of detecting $k$ photons is the derivative of order $k$ evaluated at zero ${G^{\left(k\right)}}\left(0\right)$. The factorial moment of order $m$ $E\left[{k\left({k-1}\right)...\left({k-m+1}\right)}\right]$ is equal to the $m$-th order derivative evaluated at $z=1:$ ${G^{\left(m\right)}}\left(1\right).$

[TABLE]

[TABLE]

where $E$ stands for the expected value.

Autocorrelation function of order $m,$ which can be measured in an experiment with $m$ photon detectors, is determined by factorial moments, and thus can be expressed in terms of the generating function derivative at $z=1.$

[TABLE]

Here $\mu=G^{\left(1\right)}\left(1\right)$ is the mean number of photons in the initial state. This implies the equality $g^{\left(1\right)}=1.$

We use a beam splitter with low reflection probability $p$ to separate individual photons from the beam. Using the events at the subtracted photon detector, we select the data pertaining only to the states from which we have annihilated a photon; everything else is discarded. Let the initial distribution of the number of photons have the probability mass function (pmf) $P\left(k\right).$ Taking into account the probability of exactly one photon being split off at the beam splitter $kp{\left(1-p\right)}^{k-1},$ we have the probability of obtaining the state with $k-1$ photons equal to $N\cdot P\left(k\right)kp{\left(1-p\right)}^{k-1},$ where $N$ is the normalisation constant. Thus the pmf of the initial state’s photon number distribution is modified by a factor of $kp{\left(1-p\right)}^{k-1}.$ As $p\rightarrow 0$ the factor goes to $k,$ since $p{\left(1-p\right)}^{k-1}\rightarrow 1.$ We can express photon subtraction using creation and annihilation operators:

[TABLE]

Based on this analysis we can show that, as the reflection probability goes to zero $p\rightarrow 0,$ the generating function $G_{1}\left(z\right)$ of the photon number distribution for the photon-subtracted state is determined by the derivative of the generating function of the initial state’s distribution:

[TABLE]

In cases where $p$ is not small enough to ignore the above expression becomes:

[TABLE]

Repeated application of (5) allows us to derive an expression for the probability generating function $G_{m}\left(z\right)$ of a state from which $m$ photons have been subtracted:

[TABLE]

Here $\mu_{i}$ is the mean number of photons for the state with $i$ subtracted quanta. Iterated measurement of these quantities allows us to compute autocorrelation functions of arbitrary order:

[TABLE]

In particular, the above reply implies that the second order correlation function is equal to the ratio of mean photon number in the state with one photon subtracted $\mu_{1}$ to the mean photon number in the original state $\mu:$

[TABLE]

In general, the following recurrence allows us to calculate the autocorrelation function of order $m+1$ from the function of order $m:$

[TABLE]

It is important to note that the autocorrelation characteristics of the conditional states obtained by subtracting quanta can be expressed in terms of correlations of the original state. The correlation function $g_{m}^{\left(n\right)}$ of order $n$ for the distribution of the $m$-subtracted state can be expressed through the correlation functions of the original distribution:

[TABLE]

3 Compound Poisson distribution

Poisson distribution with parameter $\lambda$ has the probability generating function:

[TABLE]

Let the parameter $\lambda$ be a random variable with Gamma probability distribution with the probability density function (pdf)

[TABLE]

where parameters $a>0,\,b>0.$ And $\Gamma\left(a\right)$ is the gamma function.

The resulting compound distribution has the following probability generating function:

[TABLE]

This compound distribution for positive values of $a$ is known as the negative binomial distribution. It has the mean of $\mu=\nicefrac{{a}}{{b}}$ and can be reparameterised as:

[TABLE]

Thus the gamma-compouded Poisson distribution has two parameters: $\mu$ is the mean number of photons and $a$ is the photon clusterisation factor. This parameter can also be interpreted as the degree of coherency: as $a$ rises, the photon distribution converges to a Poisson distribution with the same mean, which gives photon number distribution of the coherent state. It is easy to see that (15) is equal to the generating function of the Poisson distribution in the $a\rightarrow\infty$ limit. The case of $a=1$ corresponds to the thermal state.

This distribution describes a multimode thermal state, where $a$ is the number of modes [15]. It can be shown that the same distribution also applies to the single-mode multiphoton-subtracted thermal state [16].

We note that taking optical losses into account corresponds to a simple scaling of the mean photon number: $\mu\rightarrow\mu t,$ where $t$ is the transmission coefficient: $t=1-\gamma,$ with $\gamma$ equal to the absorbed proportion of energy.

Another interesting fact is that the distribution generated by (15) is well defined not only for the positive values of $a,$ but also for negative integers: $a=-n,$ where $n=1,2,\ldots.$ Such distributions can be used to model states with $n$ photons as long as the condition $0<\mu\leq n$ is satisfied. The case $\mu=n$ corresponds to the case, where the prepared state contains precisely $n$ photons. In case of $\mu<n,$ we can interpret the distribution as describing the situation after each of $n$ photos was subject to annihilation with probability of survival $\theta=\frac{\mu}{n}.$ The resulting binomial distribution with probability mass function $f\left(k;n,\theta\right)$ then gives the probability of $k$ photons surviving, out of initial set of $n.$

Using the results from section 2, it is easy to show that autocorrelation function of order $m$ for the gamma-compounded Poisson is given by:

[TABLE]

Here $\left(a\right)_{k}$ is the rising factorial:

[TABLE]

In particular: ${g^{\left(2\right)}}=\frac{{a+1}}{a},$ ${g^{\left(3\right)}}=\frac{{\left({a+1}\right)\left({a+2}\right)}}{{a^{2}}},$ ${g^{\left(4\right)}}=\frac{{\left({a+1}\right)\left({a+2}\right)\left({a+3}\right)}}{{a^{3}}}.$

For example a thermal state has $a=1$ and $g^{\left(2\right)}=2;$ one-photon Fock state has $a=-1$ and $g^{\left(2\right)}=0$ and a coherent state has $g^{\left(2\right)}=1$ with $a\rightarrow\infty.$

Using (5) it is easy to show that the conditional state obtained after photon subtraction has the following important property: the photon number distribution keeps the compound Poisson type, the parameter $b$ is unchanged and the parameter $a$ is simply incremented by one:

[TABLE]

From the above we derive the expression for the mean number of photons after subtraction:

[TABLE]

Here $a^{\left(1\right)}$ and $\mu^{\left(1\right)}$ are the parameters of the photon-subtracted state’s distribution.

After subtracting $m$ photons we obtain:

[TABLE]

[TABLE]

In case of a thermal source $a=1,$ so

[TABLE]

In particular $\mu^{\left(1\right)}=2\mu$ and after subtracting one photon we have raised the expected number of photons in the state by a factor of two. This increase does not represent a paradox, since the increased photon count is expected in the conditional state’s distribution. Subtracting a photon lowers the photon count by one, but at the same time the subtraction event is more likely for the states with higher photon counts in the original state. In case of thermal states, the second effect is powerful enough to double the expected number of photons.

4 Quadrature measurements of states

In case of the mixture with the compound-Poisson probability distribution $P\left({k\left|{\mu,a}\right.}\right)$ over Fock components we have the following density matrix in the in-phase quadrature basis:

[TABLE]

Here $\varphi_{k}\left(x\right),\;k=0,1,2,\ldots$ are the Hermite functions, which form the eigenbasis of the harmonic oscillator. These functions can be written in closed form as

[TABLE]

where $H_{k}$ are the Hermite polynomials.

The corresponding probability distribution over the quadratures is then given by

[TABLE]

This is an even function, so all the odd moments are equal to zero.

The variance of the quadrature distribution (25) does not depend on parameter $a$ and is completely determined by the mean photon count:

[TABLE]

It is also possible to obtain closed form expressions for higher order moments. In particular the skewness $\beta_{1}$ is equal to zero and the excess kurtosis $\beta_{2}$ is given by

[TABLE]

The above quantities for a random variable $x$ are defined as

[TABLE]

[TABLE]

Here $E$ stands for the expected value and $\sigma$ is the standard deviation.

It is possible to estimate the mean photon number $\mu$ and the clusterisation factor $a$ of the mixture from the variance and excess kurtosis of the quadrature data:

[TABLE]

[TABLE]

These formulae give the estimates of the parameters $\mu$ and $a,$ according to the method of moments. We use these estimates as starting points for a more powerful maximum likelihood estimation procedure.

We have prepared and measured eleven states: the thermal state and the photon-subtracted states with the number of subtracted photons $m$ varying from 1 to 10. The results are presented in 1. Here we denote the number of subtracted photons as $m$ with $m=0$ referring to the original thermal state. The error bounds are expressed using standard deviation of parameters. Large bounds for the cases with nine and ten subtracted photons are a consequence of a small sample size. Fidelity was calculated between the reconstructed state and the ideal theoretical state. $\chi^{2}$ significance test was used to check for consistency between the obtained tomographic models and the data. In our case all of the significance levels are higher than one percent and the models are consistent with the data.

According to (LABEL:autocorr_order_m), the product of all the mean photon numbers gives the autocorrelation function of the eleventh order: ${g^{\left({11}\right)}}=\frac{{\prod\limits_{m=1}^{10}{\mu_{m}}}}{{\mu^{10}}}.$ Using the data from 1 we obtain $\ln{g^{\left({11}\right)}}=17.53\pm 0.10$. The theoretical value is equal to $\ln 11!=17.50.$

5 The hierarchy of compound Poisson distributions

In section 3 of the article we let the parameter $\lambda$ of the Poisson distribution be a random variable with a gamma distribution, which gives the gamma-compound Poisson distribution. However, this model can be viewed as the first order approximation, with the original Poisson corresponding to the zeroth order. If we let the mean photon count $\mu$ of the compound Poisson be a random variable, we can continue the process of obtaining successively more sophisticated models. The physical justification for this operation is as follows. The stabilisation of the mean photon count inside the exposition time of $\tau$ is not absolute. In general the different frames are formed in different conditions and the mean time may differ from one to the other. The second level model allows us to take this variation into account. Various groups of frames may also be inhomogeneous, which would lead to a third level model and so on.

Following [17] we derive the expressions for the probability generating functions of the distributions in this hierarchy. We start from the first level gamma compound Poisson probability generating function in the following form:

[TABLE]

Now let the mean photon count $\mu$ be a random variable described by a gamma probability distribution with parameters $a_{2}$ and $b_{2},$ while the $b_{1}$ is constant. In this case we obtain the second level model. Note that we can transform a zeroth level model (12) into a first level one by means of a simple formal substitution:

[TABLE]

Averaging over the mean $\mu$ in (32) we obtain the generating function of the second level compound Poisson model with parameters $a_{2}$ and $b_{2}$:

[TABLE]

Here $\mu$ is the general expectation, $b_{1}$ and $b_{2}$ are parameters of the first and second levels of hierarchy, respectively. We can also define the clusterisation parameters for the first and second levels as $a_{1}=\mu b_{1},$ $a_{2}=\mu b_{2}.$ We note that as $a_{2}\rightarrow\infty,$ $b_{2}\rightarrow\infty,$ while $\nicefrac{{a_{2}}}{{b_{2}}}=\mu=const,$ the second level model converges to the first order one.

Iteratively repeating the above procedure we obtain the following recurrence relations between probability generating functions of various levels of the hierarchy:

[TABLE]

For each of the levels along with parameters $b_{r}$ we can also define the clusterisation parameters $a_{r}=\mu b_{r},\;r=1,2,\ldots$

The experiments demonstrate a relatively small inconsistency of the data with the first level model and the necessity of employing the second level model.

Figure 1 depicts the comparison between models for the state with one subtracted photon.

Here the experimental sample is used to form the histogram, the dashed red line is the first level model with $a^{ideal}=2$ and the solid green line corresponds to a second level model with $a_{1}=2,$ $a_{2}=8.46.$ $\mu$ is equal to $5.98$ in both cases. The second level model clearly fits the data better than the ideal theoretic model with no corrections (in particular, the corrected model has no minimum at [math]). Note that the second level corrections tend to zero as the parameter $a_{2}\rightarrow\infty.$

6 Conclusion

Photon statistics of a family of photon-subtracted thermal states have been described using the probability generating functions of the photon number distribution. The correlation functions of various orders were expressed in terms of quantities measurable using only detectors incapable of resolving the number of photons. Up to ten-photon subtracted states have been experimentally realised and measured with $>99\%$ fidelity. Our results showcase the flexibility of this model in analysis of quantum optical experimental data.

Acknowledgements.

The work was supported by the Russian Science Foundation, grant no 14-12-01338.