Full analysis of multi-photon pair effects in spontaneous parametric down conversion based photonic quantum information processing

Let us consider an n-mode bosonic system associated with an infinite-dimensional Hilbert space ${{\mathcal{H}}^{\otimes n}}$ and n pairs of annihilation and creation operators, ${{\{{{\hat{a}}_{i}},\hat{a}_{i}^{\dagger }\}}_{i=1,\cdots ,n}}$, respectively, which satisfy the commutation relations

$$\begin{eqnarray}&&[{{\hat{a}}_{i}},\hat{a}_{j}^{\dagger }]={{\delta }_{ij}}.\end{eqnarray} \tag{ 2 }$$

From these, one may construct the quadrature field operators:

$$\begin{eqnarray}&&{{\hat{x}}_{i}}=\frac{1}{\sqrt{2}}(\hat{a}_{i}^{\dagger }+{{\hat{a}}_{i}}),\quad {{\hat{p}}_{i}}=\frac{i}{\sqrt{2}}(\hat{a}_{i}^{\dagger }-{{\hat{a}}_{i}}).\end{eqnarray} \tag{ 3 }$$

It is easy to verify that the commutation relations now translate to $[{{\hat{x}}_{i}},{{\hat{p}}_{j}}]=i{{\delta }_{ij}}$. In the n-mode bosonic system, a quantum state with density operator $\hat{\rho }$ is described by its characteristic function

$$\begin{eqnarray}&&\chi (x)={\rm Tr}\left[ \hat{\rho }\hat{\mathcal{W}}(x) \right],\end{eqnarray} \tag{ 4 }$$

where,

$$\begin{eqnarray}&&\hat{\mathcal{W}}(x)={\rm exp} \left[ -{\rm i}{{x}^{T}}\hat{R} \right],\end{eqnarray} \tag{ 5 }$$

is a Weyl operator, $\hat{R}={{[{{\hat{x}}_{1}},...,{{\hat{x}}_{n}},{{\hat{p}}_{1}},...,{{\hat{p}}_{n}}]}^{T}}$ is a $2n$ vector consisting of quadrature operators, and $x=[{{x}_{1}},\cdots ,{{x}_{2n}}]$ is a $2n$ real vector.

Definition. The Gaussian state is defined as the quantum state whose characteristic function is given by a Gaussian function:

$$\begin{eqnarray}&&\chi (x)={\rm exp} \left[ -\frac{1}{4}{{x}^{T}}\gamma x-{\rm i}{{d}^{T}}x \right],\end{eqnarray} \tag{ 6 }$$

where γ is a $2n\times 2n$ matrix and d is a $2n$ vector called the covariance matrix and the displacement vector, respectively. For example, coherent state $|\alpha \rangle$ is a Gaussian state. Its characteristic function and displacement vector are described by the Gaussian form in (6) with

$$\begin{eqnarray}&&{{\gamma }^{{\rm coh}}}=I,\quad d=\sqrt{2}\left[ \begin{array}{ccccccccccccccc} \operatorname{Re}\alpha \\ {\rm Im}\alpha \\ \end{array} \right]\end{eqnarray} \tag{ 7 }$$

where I is a 2-by-2 identity matrix (thus the vacuum state is simply given by $\{\gamma =I,d=0\}$). The most important Gaussian state in this paper is the TMSV. Its covariance matrix is given by

$$\begin{eqnarray}{{\gamma }^{{\rm TMSV}}}(\mu )=\left[ \begin{array}{ccccccccccccccc} {{\gamma }^{+}}(\mu ) & {\bf 0} \\ {\bf 0} & {{\gamma }^{-}}(\mu ) \\ \end{array} \right],\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}{{\gamma }^{\pm }}(\mu )=\left[ \begin{array}{ccccccccccccccc} 2\mu +1 & \pm 2\sqrt{\mu (\mu +1)} \\ \pm 2\sqrt{\mu (\mu +1)} & 2\mu +1 \\ \end{array} \right],\end{eqnarray} \tag{ 9 }$$

and $\mu ={{\lambda }^{2}}/(1-{{\lambda }^{2}})$ while d = 0. It is worth noting that μ corresponds to the average photon number per mode, i.e. $\mu =\langle \psi |\hat{n}\otimes \hat{I}|\psi \rangle =\langle \psi |\hat{I}\otimes \hat{n}|\psi \rangle$, where $\hat{I}$ and $\hat{n}$ are the identity and photon number operators, respectively.

Partial trace of Gaussian states. The covariance matrix of the reduced state after partial trace is simply given by the submatrix corresponding to the remaining system. For example, the covariance matrix of the reduced state of the TMSV ${{\rho }_{{\rm S}}}={\rm T}{{{\rm r}}_{I}}[|\psi \rangle \langle \psi {{|}_{SI}}]$ is

$$\begin{eqnarray}{{\gamma }^{{\rm S}}}=\left[ \begin{array}{ccccccccccccccc} 2\mu +1 & 0 \\ 0 & 2\mu +1 \\ \end{array} \right],\end{eqnarray} \tag{ 10 }$$

which corresponds to the covariance matrix of a thermal state with average photon number μ:

$$\begin{eqnarray}&&{{\hat{\rho }}_{{\rm th}}}=\mathop{\sum }\limits_{n=0}^{\infty }\frac{{{\mu }^{n}}}{{{(\mu +1)}^{n+1}}}|n\rangle \langle n|.\end{eqnarray} \tag{ 11 }$$

Gaussian unitary operation is defined as the unitary operation that transforms Gaussian states to other Gaussian states. Any Gaussian unitary operation acting on a Gaussian state can be described by symplectic transformations of the covariance matrix and the displacement vector of the state:

$$\begin{eqnarray}&&\gamma \to {{S}^{T}}\gamma S,\quad d\to {{S}^{T}}d,\end{eqnarray} \tag{ 12 }$$

where S is a symplectic matrix corresponding to the Gaussian unitary operation and T is a transpose operation. For any covariance matrix, there exists a symplectic transformation that diagonalizes the covariance matrix (symplectic diagonalization). If the unitary operation includes only a linear optical process (beam splitters and phase shifts), then ${{S}^{T}}={{S}^{-1}}$ and such a matrix S is called an orthogonal symplectic matrix. The explicit expression of the symplectic matrix for phase shifting and beam splitting are given below.

Phase shift:

$$\begin{eqnarray}R(\phi )=\left[ \begin{array}{ccccccccccccccc} {\rm cos} \phi & {\rm sin} \phi \\ -{\rm sin} \phi & {\rm cos} \phi \\ \end{array} \right].\end{eqnarray} \tag{ 13 }$$

Beam splitter on mode A and B:

$$\begin{eqnarray}S_{AB}^{t}=\left[ \begin{array}{ccccccccccccccc} \sqrt{t} & \sqrt{1-t} & 0 & 0 \\ -\sqrt{1-t} & \sqrt{t} & 0 & 0 \\ 0 & 0 & \sqrt{t} & \sqrt{1-t} \\ 0 & 0 & -\sqrt{1-t} & \sqrt{t} \\ \end{array} \right],\end{eqnarray} \tag{ 14 }$$

where t is the transmittance of the beam splitter. Throughout the paper, we often simplify the description of a block diagonalized matrix like (14) as

$$\begin{eqnarray*}S_{AB}^{t}={{\left[ \begin{array}{ccccccccccccccc} \sqrt{t} & \sqrt{1-t} \\ -\sqrt{1-t} & \sqrt{t} \\ \end{array} \right]}^{\,\oplus 2}}.\end{eqnarray*}$$

A major detection device in experimental photonic QIP is a photon detector, which discriminates only zero or non-zero photons. Such a device is mathematically described by a set of measurement operators

$$\begin{eqnarray}&&{{\hat{\Pi }}_{{\rm off}}}=|0\rangle \langle 0|,\quad {{\hat{\Pi }}_{{\rm on}}}=\hat{I}-|0\rangle \langle 0|,\end{eqnarray} \tag{ 15 }$$

where $\hat{I}$ is an identity operator. Similar to the state, we can define the characteristic function of the measurement operator $\hat{\Pi }$ as ${{\chi }_{\Pi }}(x)={\rm Tr}[\hat{\Pi }\hat{\mathcal{W}}(x)]$. In general, the probability of detecting the state $\hat{\rho }$ with the measurement operator $\hat{\Pi }$ is given by

$$\begin{eqnarray}&&{\rm Tr}\left[ \hat{\rho }\hat{\Pi } \right]={{\left( \frac{1}{2\pi } \right)}^{n}}\int {\rm d}x\;{{\chi }_{\rho }}(x){{\chi }_{\Pi }}(-x).\end{eqnarray} \tag{ 16 }$$

Now suppose $\hat{\rho }$ is a single-mode Gaussian state with ${{\chi }_{\rho }}(x)={\rm exp} [-\frac{1}{4}{{x}^{T}}\gamma x]$ and is measured by an on-off detector. The probability of obtaining the 'on' outcome (i.e. detecting non-zero photons) is calculated to be

$$\begin{eqnarray}\begin{array}{rcl} {{P}_{{\rm on}}} & = & {\rm Tr}\left[ \hat{\rho }{{{\hat{\Pi }}}_{{\rm on}}} \right]=1-{\rm Tr}\left[ \hat{\rho }|0\rangle \langle 0| \right] \\ {} & = & 1-\frac{1}{2\pi }\int {\rm d}x\;{{\chi }_{\rho }}(x){{\chi }_{{{\Pi }_{{\rm off}}}}}(-x) \\ {} & = & 1-\frac{1}{2\pi }\int {\rm d}x\;{\rm exp} \left[ -\frac{1}{4}{{x}^{T}}(\gamma +I)x \right] \\ {} & = & 1-\frac{2}{\sqrt{{\rm det} (\gamma +I)}}, \\ \end{array}\end{eqnarray} \tag{ 17 }$$

where the last line performs the Gaussian integration. The above measurement is an ideal one, i.e. unit efficiency and no dark counts. We will discuss the detector loss and dark counts in the next subsection.

In the SPDC QIP experiments, we often consider the coincidence photon counts of a multi-mode quantum state. Let ${{\hat{\rho }}^{\gamma }}$ be a density matrix of an m-mode Gaussian state with covariance matrix γ. Then the following formula is useful to calculate the coincidence counts:

$$\begin{eqnarray}\begin{array}{rcl} {\rm Tr}\left[ {{{\hat{\rho }}}^{\gamma }}|0\rangle \langle 0{{|}^{\otimes m}} \right] & = & {{\left( \frac{1}{2\pi } \right)}^{m}}\int {\rm d}x\;{\rm exp} \left[ -\frac{1}{4}{{x}^{T}}(\gamma +I)x \right] \\ {} & = & {{\left( \frac{1}{2\pi } \right)}^{m}}\sqrt{\frac{{{(4\pi )}^{2\,m}}}{{\rm det} (\gamma +I)}} \\ {} & = & \frac{{{2}^{m}}}{\sqrt{{\rm det} (\gamma +I)}}. \\ \end{array}\end{eqnarray} \tag{ 18 }$$

Linear loss. The optical channel with transmittance t (i.e. $1-t$ loss) is modeled by combining the channel with a vacuum environment via a beam splitter of transmittance t and then tracing out the environment mode. The lossy channel is known as one of the Gaussian channels, i.e. it consists of Gaussian operations (but not necessarily unitary) [24]. Suppose a single-mode Gaussian state with covariance matrix γ is transmitted through a lossy channel ${{\mathcal{L}}^{t}}$ where t is the channel transmittance. Then ${{\mathcal{L}}^{t}}$ transforms the covariance matrix of the state as

$$\begin{eqnarray}&&{{\mathcal{L}}^{t}}\gamma ={{K}^{T}}\gamma K+\alpha ,\end{eqnarray} \tag{ 19 }$$

where $K=\sqrt{t}I$ and $\alpha =(1-t)I$. Note that the photon detector with efficiency η can be modeled by a lossy channel with transmittance η followed by a lossless photon detector.

Detector dark counts. The dark counts are modeled by a Poissonian process [25] or a phase insensitive amplification process [14] that provide the same expression for the on-off detector operators:

$$\begin{eqnarray}&&{{\hat{\Pi }}_{{\rm off}}}(\nu )=(1-\nu )|0\rangle \langle 0|,\quad {{\hat{\Pi }}_{{\rm on}}}(\nu )=\hat{I}-{{\hat{\Pi }}_{{\rm off}}}(\nu ),\end{eqnarray} \tag{ 20 }$$

where ν is the dark count probability.

Mode mismatch. Imperfection of mode matching at a beam splitter degrades the interference visibility. Let ξ be a phenomenological mode match factor representing the overlap between two pulses, i.e. $0\leqslant \xi \leqslant 1$ where $\xi =1$ corresponds to a perfect mode matching while $\xi =0$ means the two pulses are in completely different modes and there is no interference between them. Then the effect of mode mismatch is modeled by adding two (effective) beam splitters with transmittance ξ. Figures 1(a) and (b) depict an example of the mode mismatch between two single-mode pulses in the temporal mode and its translation into the beam splitting model (spatial mode), respectively, where the two pulses are partially overlapped. In the spatial mode A of figure 1(a), one can choose a temporal mode expansion in which one of the orthonormal modes is fully occupied by the pulse (i.e. includes A₁ and A₃). At the beam splitter, the mode expansion should be transformed such that the temporal part A₁ and A₃ are separated where these two modes should be a mixture with vacua to fulfil the normalization condition. This mode transformation is mathematically equivalent to splitting the pulse via a (virtual) beam splitter with transmittance ξ into two spatial modes A₁ and A₃ (and combining it with a vacuum from the other port of the beam splitter). The same thing happens in spatial mode B where the pulse is split into modes B₁ and B₂. Then A_i and B_i are combined at the real beam splitter where the pulses are interfered only between A₁ and B₁ (figure 1(b)).

Zoom In Zoom Out Reset image size

Here we give an example. Suppose we would like to measure the coincidence counts after interfering with the two Gaussian state pulses via a beam splitter with transmittance t. Let ${{\gamma }_{{{A}_{1}}{{B}_{1}}}}$ be the covariance matrix of the two-mode state (i.e. two pulses) before the beam splitter. Then after the beam splitting, a whole state is described by a six-mode covariance matrix

$$\begin{eqnarray}\begin{array}{ll} {} & {} & {{\tilde{\gamma }}_{{{A}_{1}}{{B}_{1}}{{A}_{2}}{{B}_{2}}{{A}_{3}}{{B}_{3}}}} \\ {} & {} & =S_{{\rm BS}}^{T}S_{{\rm MM}}^{T}({{\gamma }_{{{A}_{1}}{{B}_{1}}}}\;\otimes \;{{I}_{{{A}_{2}}{{B}_{2}}{{A}_{3}}{{B}_{3}}}}){{S}_{{\rm MM}}}{{S}_{{\rm BS}}}, \\ \end{array}\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&{{S}_{{\rm MM}}}=S_{{{B}_{1}}{{B}_{3}}}^{\xi }S_{{{A}_{1}}{{A}_{2}}}^{\xi },\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{S}_{{\rm BS}}}=S_{{{A}_{3}}{{B}_{3}}}^{t}S_{{{A}_{2}}{{B}_{2}}}^{t}S_{{{A}_{1}}{{B}_{1}}}^{t},\end{eqnarray} \tag{ 23 }$$

and the terms in the rhs are the beam splitter matrices defined in (14).

The coincidence count probability is then given by

$$\begin{eqnarray*}&&{\rm Tr}[{{\rho }^{\tilde{\gamma }}}(\hat{I}-|0\rangle {{\langle 0{{|}^{\,\otimes 3}})}_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}}(\hat{I}-|0\rangle {{\langle 0{{|}^{\,\otimes 3}})}_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}}].\end{eqnarray*}$$

Note that in our model, the mode mismatch is included by a phenomenological factor and thus we can incorporate any kinds of mode mismatch, such as temporal, spectral, spatial, etc.

Before closing the section, we briefly mention the (classical) feedforward operations. Feedforward is an important resource in QIP. For example, it allows one to implement an on-demand single-photon source from heralding of the SPDC photons [26], and even constructing a universal quantum computer [27, 40]. In the feedforward scenario, each operation in the system can be adaptively chosen according to the prior partial measurement outcomes. Therefore, to fully simulate such a system, one has to calculate all possible branches of the measurement outcomes and feedforwarded operations. As will be discussed in section 3.3, since the computational complexity of our method exponentially grows as the system size increases, with our method it is not possible to fully simulate a large scale linear optics quantum computer proposed in [27] with an efficient computation time. Note however, for relatively small size systems, it is possible to trace each feedforward branches and apply our method to calculate the probability of observing each event separately. That would be useful to simulate the currently (or near future) feasible experimental setups such as [26, 40].

The HOM interference is observed when two indistinguishable single-photons are interfered via a 50/50 beam splitter [28]. The visibility of the HOM interference can be unit only when the two photons are fully indistinguishable in any degree of freedom (temporal mode, frequency mode, etc) which is thus a useful measure to evaluate the indistinguishability of the signal and idler photons from the same SPDC source. The HOM interference test for an SPDC source is described as follows. The signal and idler pulses are interfered by a 50/50 beam splitter and the coincidence count of the two detectors are measured by scanning the time delay of the idler pulse. The HOM dip is observed when there is no time delay for the idler pulse, i.e, the overlap of the signal and the idler in time is maximum. Let $P_{{\rm min} }^{CC}$ be the coincidence count (CC) probability without the delay and $P_{{\rm mean}}^{CC}$ be the CC probability with the delay which is enough larger than the pulse width such that no interference occurs between the signal and idler pulses. Following the previous works e.g. [6, 14], we define the visibility of the HOM test by

$$\begin{eqnarray}&&{{V}_{{\rm HOM}}}=\frac{P_{{\rm mean}}^{CC}-P_{{\rm min} }^{CC}}{P_{{\rm mean}}^{CC}}.\end{eqnarray} \tag{ 24 }$$

To see the validity of our method, we simulate the HOM experiment in [6] where the pump power dependence of ${{V}_{{\rm HOM}}}$ is experimentally observed with a standard mode-locked laser (Ti:Sapphire laser with a repetition rate of 76 MHz). The experimental setup and the corresponding theoretical model are shown in figure 2(a) and (b), respectively.

Zoom In Zoom Out Reset image size

The output from the SPDC source is described by a TMSV which fully includes the multi-photon components. The transmission losses in the signal (A) and idler modes (B), that are for example caused by coupling the spatial modes into fibers, are described by beam splitters with the transmittance ${{t}_{A,B}}$, respectively. Also t_B includes the losses in the controllable delay line (the collimator in figure 2(a)) which makes the setup asymmetric (in the sense that ${{t}_{A}}\ne {{t}_{B}}$). ${{\eta }_{A}}$ and ${{\eta }_{B}}$ are the quantum efficiencies of the two photon detectors. The mode mismatch between the signal and idler pulses at the 50/50 beam splitter is characterized by the mode match factor ξ which is in fact very sensitive to the HOM interference visibility as shown below.

The covariance matrix of the state generated from the SPDC is given by $\gamma _{AB}^{{\rm TMSV}}(\mu )$ which is defined in (8). Applying losses in mode A and B with transmittance t_A and t_B, respectively, we get

$$\begin{eqnarray}\begin{array}{rcl} \gamma _{AB}^{\mathcal{L}} & = & \mathcal{L}_{B}^{{{t}_{B}}}\mathcal{L}_{A}^{{{t}_{A}}}\gamma _{AB}^{{\rm TMSV}}(\mu ) \\ {} & = & K_{AB}^{{{t}_{A}}{{t}_{B}}\,T}\gamma _{AB}^{{\rm TMSV}}(\mu )K_{AB}^{{{t}_{A}}{{t}_{B}}}+\alpha _{AB}^{{{t}_{A}}{{t}_{B}}} \\ {} & = & {{\left[ \begin{array}{ccccccccccccccc} 2{{t}_{A}}\mu +1 & \pm 2\sqrt{{{t}_{A}}{{t}_{B}}\mu (\mu +1)} \\ \pm 2\sqrt{{{t}_{A}}{{t}_{B}}\mu (\mu +1)} & 2{{t}_{B}}\mu +1 \\ \end{array} \right]}^{\,\oplus 2}}, \\ \end{array}\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}K_{AB}^{{{t}_{A}}{{t}_{B}}}={{\left[ \begin{array}{ccccccccccccccc} \sqrt{{{t}_{A}}} & 0 \\ 0 & \sqrt{{{t}_{B}}} \\ \end{array} \right]}^{\,\oplus 2}},\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}\alpha _{AB}^{{{t}_{A}}{{t}_{B}}}={{\left[ \begin{array}{ccccccccccccccc} 1-{{t}_{A}} & 0 \\ 0 & 1-{{t}_{B}} \\ \end{array} \right]}^{\,\oplus 2}}.\end{eqnarray} \tag{ 27 }$$

Application of a 50/50 beam splitter with mode matching factor ξ is calculated to be

$$\begin{eqnarray}&&\begin{array}{l} \gamma _{{{A}_{1}}{{A}_{2}}{{A}_{3}}{{B}_{1}}{{B}_{2}}{{B}_{3}}}^{{\rm BS}} \\ \quad =S_{{\rm BS}}^{T}S_{{\rm MM}}^{T}(\gamma _{{{A}_{1}}{{B}_{1}}}^{\mathcal{L}}\oplus {{I}_{{{A}_{2}}{{B}_{2}}{{A}_{3}}{{B}_{3}}}}){{S}_{{\rm MM}}}{{S}_{{\rm BS}}}, \\ \end{array}\end{eqnarray} \tag{ 28 }$$

where ${{S}_{{\rm MM}}}$ and ${{S}_{{\rm BS}}}$ are given by (22) and (23) with $t=1/2$, respectively. ${{I}_{{{A}_{2}}{{B}_{2}}{{A}_{3}}{{B}_{3}}}}$ is a 8-by-8 identity matrix representing the vacua in modes A₂, B₂, A₃ and B₃. Finally photon detection consists of lossy channels corresponding to the detector loss and perfect detectors. The covariance matrix of the system after the lossy channels with ${{\eta }_{A}}$ and ${{\eta }_{B}}$, corresponding to the quantum efficiencies of two detectors, are given by

$$\begin{eqnarray}&&{{\gamma }_{{{A}_{1}}...{{B}_{3}}}}=\mathcal{L}_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}^{{{\eta }_{B}}}\mathcal{L}_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}^{{{\eta }_{A}}}\gamma _{{{A}_{1}}...{{B}_{3}}}^{{\rm BS}}.\end{eqnarray} \tag{ 29 }$$

We are now in a position to derive the coincidence count probability:

$$\begin{eqnarray}\begin{array}{rcl} P_{{\rm min} }^{CC} & = & {\rm Tr}\left[ {{{\hat{\rho }}}^{{{\gamma }_{{{A}_{1}}...{{B}_{3}}}}}}{{\left( \hat{I}-|0\rangle \langle 0{{|}^{\,\otimes 3}} \right)}_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}}{{\left( \hat{I}-|0\rangle \langle 0{{|}^{\,\otimes 3}} \right)}_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}} \right] \\ {} & = & 1-{\rm Tr}\left[ {{{\hat{\rho }}}^{{{\gamma }_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}}}}|0\rangle \langle 0|_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}^{\,\otimes 3} \right]-{\rm Tr}\left[ {{{\hat{\rho }}}^{{{\gamma }_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}}}}|0\rangle \langle 0|_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}^{\,\otimes 3} \right]+{\rm Tr}\left[ {{{\hat{\rho }}}^{{{\gamma }_{{{A}_{1}}...{{B}_{3}}}}}}|0\rangle \langle 0|_{{{A}_{1}}...{{B}_{3}}}^{\otimes 6} \right] \\ {} & = & 1-\frac{8}{\sqrt{{\rm det} ({{\gamma }_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}}+I)}}-\frac{8}{\sqrt{{\rm det} ({{\gamma }_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}}+I)}}+\frac{64}{\sqrt{{\rm det} ({{\gamma }_{{{A}_{1}}...{{B}_{3}}}}+I)}}, \\ \end{array}\end{eqnarray} \tag{ 30 }$$

where ${{\gamma }_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}}$ and ${{\gamma }_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}}$ are the submatrices of ${{\gamma }_{{{A}_{1}}...{{B}_{3}}}}$. The determinants in (30) are explicitly given by

$$\begin{eqnarray}&&\begin{array}{l} \sqrt{{\rm det} ({{\gamma }_{{{A}_{1}}{{A}_{2}}{{A}_{3}}}}+I)} \\ =8{{\left\{ {{\left[ 1+\frac{{{\eta }_{A}}\mu }{4}\left\{ 2({{t}_{A}}+{{t}_{B}})-{{t}_{A}}{{t}_{B}}{{\eta }_{A}}(1-{{\xi }^{2}}) \right\} \right]}^{2}}-{{t}_{A}}{{t}_{B}}\eta _{A}^{2}{{\xi }^{2}}\mu (\mu +1) \right\}}^{1/2}}, \\ \end{array}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&\begin{array}{l} \sqrt{{\rm det} ({{\gamma }_{{{B}_{1}}{{B}_{2}}{{B}_{3}}}}+I)} \\ =8{{\left\{ {{\left[ 1+\frac{{{\eta }_{B}}\mu }{4}\left\{ 2({{t}_{A}}+{{t}_{B}})-{{t}_{A}}{{t}_{B}}{{\eta }_{B}}(1-{{\xi }^{2}}) \right\} \right]}^{2}}-{{t}_{A}}{{t}_{B}}\eta _{B}^{2}{{\xi }^{2}}\mu (\mu +1) \right\}}^{1/2}}, \\ \end{array}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\begin{array}{l} \sqrt{{\rm det} ({{\gamma }_{{{A}_{1}}...{{B}_{3}}}}+I)} \\ =64\left[ {{\left\{ 1+\frac{\mu }{4}\left[ 2({{t}_{A}}+{{t}_{B}})({{\eta }_{A}}+{{\eta }_{B}})-{{t}_{A}}{{t}_{B}}\left\{ {{({{\eta }_{A}}+{{\eta }_{B}})}^{2}}-{{\xi }^{2}}{{({{\eta }_{A}}-{{\eta }_{B}})}^{2}} \right\} \right] \right\}}^{2}} \right. \\ {{\left. \qquad -{{t}_{A}}{{t}_{B}}{{({{\eta }_{A}}-{{\eta }_{B}})}^{2}}{{\xi }^{2}}\mu (\mu +1) \right]}^{1/2}}. \\ \end{array}\end{eqnarray} \tag{ 33 }$$

The coincidence count with delay, $P_{{\rm mean}}^{CC}$ is simply obtained from $P_{{\rm min} }^{CC}$ by setting $\xi =0$, i.e. no overlap between the two modes.

The theoretical visibilities derived above and the experimental plots in [6] are compared in figure 3 where the two channel transmittances were experimentally measured as t_A = 0.42, t_B = 0.29. The detectors used in the experiment were superconducting nanowire single photon detectors (SNSPDs) that gave ${{\eta }_{A}}=0.68$, and ${{\eta }_{B}}=0.70$ and negligible dark counts (less than 1k counts per second). The squeezing parameter ${{\lambda }^{2}}$ corresponds to the photon-pair generation rate which is directly measurable in the experiment (recall that ${{\lambda }^{2}}=\mu /(1+\mu )$). The mode matching factor is not measured and thus given as a fitting parameter. The figure shows a quantitative agreement between the theory and the experiment. Also the theory curves reveal that the visibility is highly sensitive to the mode matching at the beam splitter. Note that [6] also provides the theoretical curves, that are obtained by directly calculating the state vector of a whole system including environments and numerically computed the coincidence counts. One can numerically check that the numerical result in [6] and the theoretical curve in figure 3 are exactly identical. However, our method is much simpler and more systematic than the approach in [6] and allows us to obtain a closed formula. From that closed formula, we can also derive a simple analytical picture of the visibility. Suppose experimental imperfection is only the low efficiency of the detectors ${{\eta }_{A}}={{\eta }_{B}}=\eta \ll 1$ and ${{t}_{A}}={{t}_{B}}=\xi =1$. Assuming $\mu \ll 1$, which is reasonable for the photon-pair experiments, the visibility is simply given by

$$\begin{eqnarray}&&{{V}_{{\rm HOM}}}\approx \frac{1+2\mu }{1+4\mu }\approx 1-2\mu .\end{eqnarray} \tag{ 34 }$$

Zoom In Zoom Out Reset image size

It should be noted that (34) agrees with the one obtained in [14]. Another interesting limit is the ideal case (${{t}_{A}}={{t}_{B}}={{\eta }_{A}}={{\eta }_{B}}=\xi =1$ and $\mu \ll 1$ ):

$$\begin{eqnarray}&&{{V}_{{\rm HOM}}}\approx \frac{2+2\mu }{2+3\mu }\approx 1-\frac{1}{2}\mu ,\end{eqnarray} \tag{ 35 }$$

which shows that the visibility degradation due to the multi-photon emission is accelerated by the low efficiency of the detectors.

Another example is the EPR interference test (also called correlation measurement) for polarization entangled photon pairs. We consider the Sagnac loop entanglement source, which was proposed in [29] and is now widely used for photonic QIP experiments [18, 30–33]. Here we compare our model with the EPR interference experiment performed by a Sagnac loop entanglement source in [18]. The experimental set up and the corresponding theoretical model are illustrated in figures 4(a) and (b), respectively.

Zoom In Zoom Out Reset image size

The Sagnac loop source generates a TMSV into the horizontal and vertical polarization modes in the clockwise direction (labeled as ${{A}_{H}},{{A}_{V}}$ in the theoretical model) and the counter-clockwise direction (${{B}_{H}},{{B}_{V}}$). The entangled state is generated by swapping the vertical polarization modes A_V and B_V via a dichroic polarization beam splitter (DPBS). This operation transforms the covariance matrix of the two TMSV states $\gamma _{{{A}_{H}}{{A}_{V}}{{B}_{H}}{{B}_{V}}}^{{\rm TMSV}2}$ into $\gamma _{{{A}_{H}}{{A}_{V}}{{B}_{H}}{{B}_{V}}}^{{\rm SL}}$ where

$$\begin{eqnarray}\gamma _{{{A}_{H}}{{A}_{V}}{{B}_{H}}{{B}_{V}}}^{{\rm TMSV}2}={{\left[ \begin{array}{ccccccccccccccc} 2\mu +1 & \pm 2\sqrt{\mu (\mu +1)} & 0 & 0 \\ \pm 2\sqrt{\mu (\mu +1)} & 2\mu +1 & 0 & 0 \\ 0 & 0 & 2\mu +1 & \pm 2\sqrt{\mu (\mu +1)} \\ 0 & 0 & \pm 2\sqrt{\mu (\mu +1)} & 2\mu +1 \\ \end{array} \right]}^{\,\oplus 2}},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}\gamma _{{{A}_{H}}{{A}_{V}}{{B}_{H}}{{B}_{V}}}^{{\rm SL}}={{\left[ \begin{array}{ccccccccccccccc} 2\mu +1 & 0 & 0 & \pm 2\sqrt{\mu (\mu +1)} \\ 0 & 2\mu +1 & \pm 2\sqrt{\mu (\mu +1)} & 0 \\ 0 & \pm 2\sqrt{\mu (\mu +1)} & 2\mu +1 & 0 \\ \pm 2\sqrt{\mu (\mu +1)} & 0 & 0 & 2\mu +1 \\ \end{array} \right]}^{\,\oplus 2}}.\end{eqnarray} \tag{ 37 }$$

The EPR interference test is a common and relatively easy way to experimentally estimate the quality of entanglement. Each spatial mode is projected onto a particular polarization basis by a half-wave plate (HWP) and a polarization beam splitter (PBS) and then the photons in a chosen polarization are detected by a single photon detector. The coincidence photon count rate depends on the angle of each polarization basis and the interference fringe is obtained by fixing one of the polarization angles constant and rotating the other one. The visibility is given by

$$\begin{eqnarray}&&{{V}_{{\rm ent}}}=\frac{P_{{\rm max} }^{CC}-P_{{\rm min} }^{CC}}{P_{{\rm max} }^{CC}+P_{{\rm min} }^{CC}},\end{eqnarray} \tag{ 38 }$$

where $P_{{\rm max} }^{CC}$ and $P_{{\rm min} }^{CC}$ are the maximum and minimum count rates in the fringe.

The polarizer (HWP and PBS) with angle θ effectively works as a beam splitter between the horizontal and vertical polarization modes with transmittance ${{{\rm cos} }^{2}}\theta$. The main imperfection in this component is a finite extinction ratio in the PBS which is modeled by a perfect beam splitter followed by losses in each polarization mode with ${{t}_{H}}\lt 1$ and ${{t}_{V}}\gt 0$ where t_H and t_V are the transmittances of the horizontal and vertical polarization modes, respectively (see figure 5(a) and (b)). In the following, to simplify the calculation of the setup in figure 4(b), we use ${{\tilde{\eta }}_{H}}={{t}_{H}}{{\eta }_{H}}$ and ${{\tilde{\eta }}_{V}}={{t}_{V}}{{\eta }_{V}}$ where ${{\eta }_{H}}$ and ${{\eta }_{V}}$ are the quantum efficiencies of the detector for horizontally and vertically polarized photons, respectively. The covariance matrix of the entangled source in (37) is first transformed by perfect beam splitter operations $S_{{{B}_{H}}{{B}_{V}}}^{{{\theta }_{B}}}S_{{{A}_{H}}{{A}_{V}}}^{{{\theta }_{A}}}$. The imperfection of the PBSs and the imperfect quantum efficiency of the detectors are included by applying the lossy channels $\mathcal{L}_{{{B}_{V}}}^{{{\tilde{\eta }}_{{{B}_{V}}}}}\mathcal{L}_{{{B}_{H}}}^{{{\tilde{\eta }}_{{{B}_{H}}}}}\mathcal{L}_{{{A}_{V}}}^{{{\tilde{\eta }}_{{{A}_{V}}}}}\mathcal{L}_{{{A}_{H}}}^{{{\tilde{\eta }}_{{{A}_{H}}}}}$. The coincidence count is then given by

$$\begin{eqnarray}\begin{array}{rcl} P_{{{\theta }_{A}}{{\theta }_{B}}}^{CC} & = & {\rm Tr}\left[ {{{\hat{\rho }}}^{{{\gamma }_{{{A}_{H}}{{A}_{V}}{{B}_{H}}{{B}_{V}}}}}}{{\left( \hat{I}-|0\rangle \langle 0{{|}^{\,\otimes 2}} \right)}_{{{A}_{H}}{{A}_{V}}}}{{\left( \hat{I}-|0\rangle \langle 0{{|}^{\,\otimes 2}} \right)}_{{{B}_{H}}{{B}_{V}}}} \right] \\ {} & = & 1-\frac{4}{\sqrt{{\rm det} ({{\gamma }_{{{A}_{H}}{{A}_{V}}}}+I)}}-\frac{4}{\sqrt{{\rm det} ({{\gamma }_{{{B}_{H}}{{B}_{V}}}}+I)}}+\frac{16}{\sqrt{{\rm det} ({{\gamma }_{{{A}_{H}}{{A}_{V}}{{B}_{H}}{{B}_{V}}}}+I)}}, \\ \end{array}\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} \sqrt{{\rm det} ({{\gamma }_{{{A}_{H}}{{A}_{V}}}}+I)} & = & 4(1+{{\tilde{\eta }}_{{{A}_{H}}}}\mu )(1+{{\tilde{\eta }}_{{{A}_{V}}}}\mu ), \\ \end{array}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}\begin{array}{rcl} \sqrt{{\rm det} ({{\gamma }_{{{B}_{H}}{{B}_{V}}}}+I)} & = & 4(1+{{\tilde{\eta }}_{{{B}_{H}}}}\mu )(1+{{\tilde{\eta }}_{{{B}_{V}}}}\mu ), \\ \end{array}\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}\begin{array}{rcl} \sqrt{{\rm det} ({{\gamma }_{{{A}_{H}}{{A}_{V}}{{B}_{H}}{{B}_{V}}}}+I)} & = & 8\left[ \left\{ 1+({{\tilde{\eta }}_{{{A}_{H}}}}+{{\tilde{\eta }}_{{{B}_{H}}}}-{{\tilde{\eta }}_{{{A}_{H}}}}{{\tilde{\eta }}_{{{B}_{H}}}})\mu \right\}\left\{ 1+({{\tilde{\eta }}_{{{A}_{V}}}}+{{\tilde{\eta }}_{{{B}_{V}}}}-{{\tilde{\eta }}_{{{A}_{V}}}}{{\tilde{\eta }}_{{{B}_{V}}}})\mu \right\} \right. \\ {} & {} & +\left\{ 1+({{\tilde{\eta }}_{{{A}_{H}}}}+{{\tilde{\eta }}_{{{B}_{V}}}}-{{\tilde{\eta }}_{{{A}_{H}}}}{{\tilde{\eta }}_{{{B}_{V}}}})\mu \right\}\left\{ 1+({{\tilde{\eta }}_{{{A}_{V}}}}+{{\tilde{\eta }}_{{{B}_{H}}}}-{{\tilde{\eta }}_{{{A}_{V}}}}{{\tilde{\eta }}_{{{B}_{H}}}})\mu \right\} \\ {} & {} & \left. +({{\tilde{\eta }}_{{{A}_{H}}}}-{{\tilde{\eta }}_{{{A}_{V}}}})({{\tilde{\eta }}_{{{B}_{H}}}}-{{\tilde{\eta }}_{{{B}_{V}}}})\mu (\mu +1){\rm cos} 2({{\theta }_{A}}+{{\theta }_{B}}) \right]. \\ \end{array}\end{eqnarray} \tag{ 42 }$$

The minimum and maximum count rates required in (38) are for example obtained by $P_{{\rm min} }^{CC}=P_{0\,0}^{CC}$ and $P_{{\rm max} }^{CC}=P_{0\,\pi /2}^{CC}$ [18].

Zoom In Zoom Out Reset image size

Figure 6 shows the comparison with the experiment in [18] and our theoretical model. Note that the photon-pair generation rate of the Sagnac loop source should be $2{{\lambda }^{2}}$ instead of ${{\lambda }^{2}}$ since the SPDC crystal is pumped twice (clockwise and counter-clockwise directions). In the experiment, the overall efficiencies for the horizontal polarization are measured to be ${{\tilde{\eta }}_{{{A}_{H}}}}={{\tilde{\eta }}_{{{B}_{H}}}}=0.1$ (again, SNSPDs are used as detectors and thus we neglect the effect of dark counts). The transmittance of the vertical polarization at the PBS is not measured and thus we use it as a fitting parameter varying from 0.007 to 0.011 ($\leqslant 0.01$ is guaranteed for a typical commercial PBS). The theoretical estimate agrees with the experimental result.

Zoom In Zoom Out Reset image size

Again, it is worth further simplifying the closed form given above. By setting ${{\tilde{\eta }}_{{{A}_{V}}}}={{\tilde{\eta }}_{{{B}_{V}}}}=0$ and ${{\tilde{\eta }}_{{{A}_{H}}}}={{\tilde{\eta }}_{{{B}_{H}}}}=\eta$, the closed form of the visibility is reduced to

$$\begin{eqnarray}&&{{V}_{{\rm ent}}}=\frac{1+\mu }{1+3\mu +2\eta (2-\eta ){{\mu }^{2}}}.\end{eqnarray} \tag{ 43 }$$

In the limit of weak pumping ($\mu \ll 1$), we have

$$\begin{eqnarray}&&{{V}_{{\rm ent}}}\approx 1-2\mu ,\end{eqnarray} \tag{ 44 }$$

regardless of the detection efficiency η. Those results agree with the previous analyses in [2, 12, 14] ¹ .

The last application in this section is the concatenated entanglement swapping (CES), where we demonstrate how to apply our method to the QIP protocols with multiple use of entangled sources and achieve a drastic improvement of the computation time compared to the previous approach. Entanglement transmission over long distance is limited by the loss and noises in channel and detectors. A quantum repeater can overcome this limit but it is still challenging to implement it on large scale with the current technology. Concatenation of entanglement swapping (also called quantum relay) is known as another protocol to extend the distance of the entanglement distribution albeit with a resource overhead that is exponential with the distance (which is usually observed as an exponential decrease of the success probability) [34]. However, since entanglement swapping has been demonstrated with the current technology (see [33] and references therein), it is still interesting to see the practical performance of the CES protocol in theory.

The protocol considered here is similar to the CES model discussed in detail in [15] where the state vector evolution of the two-mode squeezed state generated from SPDCs are calculated and the authors showed explicit expression of the density matrix and the detection probabilities that are used to derive the visibility of the EPR interference of the distributed entanglement numerically. Figure 7(a) illustrates the schematic of the protocol. Entangled photon-pairs generated from the entangled photon-pair sources (EPSs) are swapped by the Bell measurement consisting of a 50/50 beam splitter, two PBSs and four on-off type photon detectors (figure 7(b)). Entanglement swapping succeeds when photons are detected by particular two detectors in the Bell measurement for example at q and s in figure 7(b). For more details of the protocol see [15] and the references therein.

Zoom In Zoom Out Reset image size

Figure 7(c) is our model which is equivalent to figure 7(a) but the orthogonal polarization modes are explicitly illustrated by different lines. For our EPS, we assume the use of the Sagnac loop SPDC entangled source discussed in the previous subsection while the method can be applied to any other sources not even necessarily entangled in polarization modes (such as time-bin entanglements). For the system imperfections, we follow the assumptions used in [15], i.e. losses in channels and detectors and the dark counts at detectors exist whereas the mode matching at beam splitters and the PBS devices are assumed to be perfect. The polarization modes are labeled by numbers from left to right (see figure 7(c)). Also in the figure, for simplicity, only two detectors are illustrated for each Bell measurement where the simultaneous clicks at these detectors correspond to the successful CES. At the left and right ends, two polarizers (see figure 5) are placed to measure the EPR interference visibility of the swapped state.

We now calculate the joint detection probability of all detectors illustrated in figure 7(c) as a function of two polarizer angles ${{\theta }_{A}}$ and ${{\theta }_{B}}$. Let ${{m}_{{\rm BM}}}$ be the number of Bell measurements (the number of concatenation), i.e. there are ${{m}_{{\rm BM}}}+1$ Sagnac loop sources (EPSs). The covariance matrix for the quantum state generated from the far left Sagnac loop SPDC is $\gamma _{1234}^{{\rm SL}}$ as discussed in (37). The total quantum state generated from ${{m}_{{\rm BM}}}+1$ Sagnac loops is thus simply given by ${{\gamma }^{{\rm SL}\;\,\;\oplus \;({{{\rm m}}_{\ {\rm BM}}}+1)}}$, more precisely,

$$\begin{eqnarray}{{\left[ \begin{array}{ccccccccccccccc} \gamma _{\pm }^{{\rm SL}}(\mu ) & 0 & \cdots \\ 0 & \gamma _{\pm }^{{\rm SL}}(\mu ) & {} \\ \vdots & {} & \ddots \\ \end{array} \right]}^{\,\oplus 2}},\end{eqnarray} \tag{ 45 }$$

where ${{\gamma }^{{\rm SL}}}\equiv {{\left( \gamma _{\pm }^{{\rm SL}}(\mu ) \right)}^{\,\oplus 2}}$ and thus $\gamma _{\pm }^{{\rm SL}}(\mu )$ correspond to the two quadrature components of ${{\gamma }^{{\rm SL}}}$. Following the discussion in [10, 15], all transmission loss is included in the detector efficiency. We also assume that all transmission channels have the same loss and all detectors are identical. Therefore, we first apply the beam splitting operations at the Bell measurements. The beam splitters combine modes $4x-1$ and $4x+1$ for horizontal polarizations and modes $4x$ and $4x+2$ for vertical polarizations where $x=1,...,{{m}_{{\rm BM}}}$. Thus the symplectic matrix applied to ${{\gamma }^{{\rm SL}\;\,\;\oplus \;({{{\rm m}}_{\ {\rm BM}}}+1)}}$ is

$$\begin{eqnarray}&&{{S}_{{\rm BM}}}={{\mathop{\prod }\limits_{x=1}^{{{m}_{{\rm BM}}}}}^{{}}}S_{4x-1,4x+1}^{1/2}S_{4x,4x+2}^{1/2}.\end{eqnarray} \tag{ 46 }$$

The two polarizer operations at the end of the concatenation are also given by applying the beam splitter symplectic matrix

$$\begin{eqnarray}&&{{S}_{{\rm AB}}}=S_{1,2}^{{{\theta }_{A}}}S_{4{{m}_{{\rm BM}}}+3,4{{m}_{{\rm BM}}}+4}^{{{\theta }_{B}}}.\end{eqnarray} \tag{ 47 }$$

Let η be the detector efficiency (including the channel transmittance). This is included by applying

$$\begin{eqnarray}&&\mathcal{L}_{{\rm tot}}^{\eta }=\mathcal{L}_{4{{m}_{{\rm BM}}}+4}^{\eta }\cdots \mathcal{L}_{1}^{\eta },\end{eqnarray} \tag{ 48 }$$

to the covariance matrix of the state. In total, we have the covariance matrix

$$\begin{eqnarray}&&\gamma _{{{m}_{{\rm BM}}}}^{CES}=\mathcal{L}_{{\rm tot}}^{\eta }\left( S_{{\rm AB}}^{T}S_{{\rm BM}}^{T}{{\gamma }^{{\rm SL}\,\oplus \,({{{\rm m}}_{\ {\rm BM}}}+1)}}{{S}_{{\rm BM}}}{{S}_{{\rm AB}}} \right).\end{eqnarray} \tag{ 49 }$$

The joint detection probability of all detectors illustrated in figure 7(c) is now obtained by calculating

$$\begin{eqnarray}\begin{array}{rcl} P_{{{\theta }_{A}}{{\theta }_{B}}}^{CES} & = & {\rm Tr}\left[ {{{\hat{\rho }}}^{\gamma _{{{m}_{{\rm BM}}}}^{CES}}}{{\left( \hat{I}-(1-\nu )|0\rangle \langle 0| \right)}_{1}}{{\left( \hat{I}-(1-\nu )|0\rangle \langle 0| \right)}_{4{{m}_{{\rm BM}}}+3}} \right. \\ {} & {} & \times \left. {{\mathop{\prod }\limits_{x=1}^{{{m}_{{\rm BM}}}}}^{{}}}{{\left( \hat{I}-(1-\nu )|0\rangle \langle 0| \right)}_{4x-1}}{{\left( \hat{I}-(1-\nu )|0\rangle \langle 0| \right)}_{4x+2}} \right] \\ {} & = & 1+\mathop{\sum }\limits_{m=1}^{M}\mathop{\sum }\limits_{C(M,m)}^{_{M}{{C}_{m}}}\frac{{{\{-2(1-\nu )\}}^{m}}}{\sqrt{{\rm det} (\gamma _{C(M,m)}^{{\rm CES}}+I)}}, \\ \end{array}\end{eqnarray} \tag{ 50 }$$

where ν is the dark count rate of detectors, M is the total number of clicked detectors, i.e. $M=2{{m}_{{\rm BM}}}+2$, and $\gamma _{C(M,m)}^{{\rm CES}}$ is the submatrix of ${{\gamma }^{{\rm CES}}}$ taking only modes $C(M,m)$. $C(M,m)$ is the m-detector combination from M detectors. For example, when ${{m}_{{\rm BM}}}=1$, the total number of detectors is M = 4 and the detectors are placed in modes 1, 3, 6 and 7. Then $C(4,1)=1$, 3, 6, 7, $C(4,2)=13$, 16, 17, 36, 37, 67, $C(4,3)=136$, 137, 167, 367, and $C(4,4)=1367$. For reference, we describe a detailed structure of $\gamma _{{{{\rm m}}_{\ {\rm BM}}}}^{CES}$ in the appendix.

In figure 8, the EPR interference visibilities are plotted as a function of μ for ${{m}_{{\rm BM}}}=1,...,5$. For comparison, we choose the parameters used in [15], $\eta =0.04$ and $\nu ={{10}^{-5}}$ for the transmittance and dark counts, respectively. In [15], the CES with ${{m}_{{\rm BM}}}=3$ is numerically simulated by a super computer with the photon number truncation at 3 photons in each mode (the authors also reported that it takes 6 hours to get a single plot for the same simulation by a single-core use of a commercially available computer). In figure 8, the curve for ${{m}_{{\rm BM}}}=3$ (brown line) consists of 400 plots that are calculated by a commercial computer with only 10 seconds without any photon number truncation² . Similarly, 100 plots for ${{m}_{{\rm BM}}}=5$ (blue line) require around two minutes with the same computer (it should be noted that our method works only for on-off detectors while the mathematical formula in [15] includes both for on-off and photon number resolving detectors).

Zoom In Zoom Out Reset image size

Our method is useful for estimating the performance of the CES for the long-distance entanglement transmission. In our model the total transmittance consists of the detector efficiency and the channel transmission as $\eta ={{\eta }_{{\rm D}}}\times {{10}^{-\alpha L/10}}$ where ${{\eta }_{{\rm D}}}$ is the detector efficiency, α is the loss coefficient, and L is the distance of the transmission channel. In the following, we consider the photon pairs at 1550 nm and use of a standard optical fiber with $\alpha =0.2$ as a transmission channel. Figure 9 compares the visibility for different ${{m}_{{\rm BM}}}$ where the total distance of the channel is fixed to be 1000 km. The detector efficiency and dark counts are assumed to be ${{\eta }_{{\rm D}}}=0.7$ and $\nu ={{10}^{-5}}$, respectively, could be typical parameters for SNSPDs discussed in previous subsections. For example, for ${{m}_{{\rm BM}}}=1$, the channel is divided into four arms and each arm has 250 km distance. In figure 9, the visibility is plotted for the CES up to ${{m}_{{\rm BM}}}=5$. First it shows that for ${{m}_{{\rm BM}}}=1$, it is completely impossible to distribute entanglement and even for ${{m}_{{\rm BM}}}=2$, the visibility must be around 0.2. Second, for extremely small μ, the figure shows that ${{m}_{{\rm BM}}}$ should be large as much as possible while for the most region of μ, the curve for ${{m}_{{\rm BM}}}=5$ decreases quickly as μ increases and ${{m}_{{\rm BM}}}=3$ or 4 are optimal for most of μ. The result suggests that to construct a practical CES network, the optimal number of concatenation should be carefully chosen by considering a given distance (channel loss) and detector parameters.

Zoom In Zoom Out Reset image size

Finally, it is worth mentioning a general scaling of the computational complexity of our method with respect to the number of modes N (or equivalently number of detectors to be detected). Our method is not efficient in the sense that the computational complexity grows exponentially with respect to N. This is clearly observed in (50). There are two steps to compute (50). First, one has to calculate the covariance matrix including losses and beam splitters (in general any linear optics and even squeezing operations). This step consists of multiplications of square matrices and each of those are polynomial (at worst order of $O({{N}^{3}})$ by direct calculation). Second, the determinants of the (sub)matrices in each term of (50) should be calculated. While determinants can be calculated for example with $O({{N}^{3}})$ by the LU decomposition [35], the number of terms in the second line of (50) increases as ${{2}^{N}}$ since the first line of (50) includes N multiplication of binomial terms $(\hat{I}-(1-\nu )|0\rangle \langle 0|)$. This makes the method inefficient (at least $O({{2}^{N}})$) for simulating large scale networks. Note that this is not only for the CES protocols but a general property for any linear optics networks with SPDCs and photon detectors. However, we should stress that, as shown in this subsection, for an experimentally feasible size of network with current technology or even larger sizes, our method can work as a powerful tool to estimate the performance of the protocols with a practical computational time.