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

In spontaneous parametric down conversion (SPDC) based quantum information processing (QIP) experiments, there is a tradeoff between the coincide count rates (i.e. the pumping power of the SPDC), which limits the rate of the protocol, and the visibility of the quantum interference, which limits the quality of the protocol. This tradeoff is mainly caused by the multi-photon pair emissions from the SPDCs. In theory, the problem is how to model the experiments without truncating these multi-photon emissions while including practical imperfections. In this paper, we establish a method to theoretically simulate SPDC based QIPs which fully incorporates the effect of multi-photon emissions and various practical imperfections. The key ingredient in our method is the application of the characteristic function formalism which has been used in continuous variable QIPs. We apply our method to three examples, the Hong-Ou-Mandel interference and the Einstein-Podolsky-Rosen interference experiments, and the concatenated entanglement swapping protocol. For the first two examples, we show that our theoretical results quantitatively agree with the recent experimental results. Also we provide the closed expressions for these the interference visibilities with the full multi-photon components and various imperfections. For the last example, we provide the general theoretical form of the concatenated entanglement swapping protocol in our method and show the numerical results up to 5 concatenations. Our method requires only a small computation resource (few minutes by a commercially available computer) which was not possible by the previous theoretical approach. Our method will have applications in a wide range of SPDC based QIP protocols with high accuracy and a reasonable computation resource.


I. INTRODUCTION
Spontaneous parametric down conversion (SPDC) is one of the most standard tools in photonic quantum information processing (QIP), such as quantum key distribution (QKD) [1], linear optics quantum computation [2] and multi-photon interferometry [3].
Toward implementing higher rate entanglement-based QKD or larger scale QIP protocols, it is important to increase the photon-pair generation rate from the SPDC source such that it provides reasonable coincidence counts of photons in multiple detectors. Though this is possible by simply increasing the pump power into the SPDC crystals, it simultaneously degrades the quantum interference visibility due to the unwanted multi-photon emissions. Therefore, on one hand, it is an important experimental topic how to reduce the effect of multi-photon emissions while keeping the higher generation rates. The experimental progress in this direction has been reported recently [4][5][6][7][8].
On the other hand, in theory, it is desirable to establish a method which fully incorporates the multi-photon emissions and various practical imperfections of the experiment, and is also applicable to complicated optical circuits for various QIP applications. The quantum state generated into the signal and idler modes from an SPDC source is described by a two-mode squeezed vacuum (TMSV): where |n represents an n-photon state, λ = sinh 2 r/(1 + sinh 2 r), and r is the squeezing parameter. Obviously it includes infinitely higher order photons. Also, to simulate experiments precisely, one has to take into account var-ious imperfections such as losses in channels and detectors, dark counts of detectors, mode mismatch between the pulses, and so on. The model including these effects is in principle calculable by directly computing the evolution of the state vector of the TMSVs [7,[9][10][11][12]. However, the drawback of this approach is that the mathematical expression of the state or count rates become very complicated even with a relatively simple setup. This often causes some difficulty in deriving a closed formula for the coincidence counts, which is usually measured in the experiments, and taking into account all practical imperfections for the precise simulation of the experiments. To circumvent such complexity, a method to directly calculate the count rates from ideal source, channel, and nonideal photon detectors, was developed in [13]. It also showed closed formulas for several experimental paramaters including the Hong-Ou-Mandel (HOM) interference and the Einstein-Podolsky-Rosen (EPR) interference visibilities where multi-photon pair effects and detector imperfections are successfully incorporated.
In this paper, we propose an alternative approach to compute the SPDC based QIP experiments that fully involve the multi-photon effects, detector imperfections, and also other practical imperfections, e.g. mode mismatching between the pulses from SPDC sources. Our method is inspired by the so-called continuous variable-QIP (CV-QIP) rather than the photonic qubit (i.e. discrete variable) QIPs. The key idea behind our method is to use the fact that the TMSV is a CV Gaussian state, i.e. whose characteristic function is a Gaussian function, and the linear optics used in QIP (phase shifters and beam splitters, etc) as well as linear imperfections (losses, mode-mismatch, optics imperfections, etc) are also Gaussian operations (see below their definitions) [14]. This implies that even for a photonic qubit based QIP with SPDC, a whole state before the detectors is still a Gaus-sian state. In CV-QIP, it is known that the Gaussian states and their evolution by Gaussian operations are calculated by tracing the covariance matrix of their characteristic functions, which drastically reduce the complexity of the calculation [14].
The single photon detection is a non-Gaussian operation. However, we will show that the Gaussian operation formalism is still valid. This comes from the fact that the standard photon detector used in the SPDC based QIP is a so-called on-off detector or threshold detector, which discriminates only zero or non-zero photons. Such a detector is mathematically described by two operators, one is a vacuum (zero photons) and the other is an identity operator minus vacuum (non-zero photons). Since vacuum is a Gaussian state, one can compute the full process of the protocols in the Gaussian characteristic function formalism (this fact is also well recognized in the CV-QIP, for example see [14][15][16]).
We apply our method to the recent experimental observations on the pump power dependence of the HOM interference visibility [8] and the EPR interference visibility [17] and showing very precise agreements with the experiment and the theory. Also, we are able to show the closed form of these visibilities as a function of the pump power (squeezing parameter) and almost complete practical imperfections (such as mode mismatch and optics and detector imperfections). Our method could be useful to simulate more complicated SPDC based QIP protocols protocols simply and precisely.
The paper is organized as follows. In Sec. II, basic definition, notation, and examples of the characteristic function formalism are reviewed. A reader who is familiar with the formalism could skip this section. In Sec. III, we apply our method into two QIP experiments, the HOM interference and the EPR interference. Section IV is the conclusion.

II. CHARACTERISTIC FUNCTIONS AND GAUSSIAN STATES
In this section, we give a brief review of the characteristic function formalism for Gaussian states and operations (for more details, see [14,18,19] for example).

A. Characteristic function
Let us consider an n-mode bosonic system associated with an infinite-dimensional Hilbert space H ⊗n and N pairs of annihilation and creation operators, ,··· ,n , respectively, which satisfy the commutation relations From these, one may construct the quadrature field operators: It is easy to verify that the commutation relations now translate to [x i ,p j ] = iδ ij . In the n-mode bosonic system, a quantum state with density operatorρ is described by its characteristic function where,Ŵ is a Weyl operator,R = [x 1 , . . . ,x n ,p 1 , . . . ,p n ] T is a 2n vector consisting of quadrature operators, and x = [x 1 , · · · , x 2n ] is a 2n real vector.

B. Gaussian states
Gaussian states. The characteristic function of any Gaussian state is represented by where 2n × 2n matrix γ and 2n vector d are called the covariance matrix and the displacement vector, respectively. For example, coherent state |α is a Gaussian state. Its characteristic function and displacement vector are given by where I is a 2-by-2 identity matrix (also a vacuum state is given by {γ = I, d = 0}). The most important Gaussian state in this paper is the TMSV. Its characteristic function is given by where γ ± (r) = cosh 2r ± sinh 2r ± sinh 2r cosh 2r , and r is the squeezing parameter (note that d = 0). The average photon number per mode for such a state isn = ψ|n ⊗Î|ψ = ψ|Î ⊗n|ψ = sinh 2 r. Partial trace of Gaussian states. The covariance matrix of the reduce state after partial trace is simply given by the submatrix corresponding to the remained system. For example, the covariance matrix of the reduce state of the TMSV ρ S = Tr I [|ψ ψ|] is where cosh 2r = 2n + 1. This corresponds to the covariance matrix of a thermal state with average photon numbern.

C. Gaussian unitary operations
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: 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 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 matrices are given below. Phase shift: Beam splitter on mode A and B: Throughout the paper, we often simplify the description of a block diagonalized matrix like Eq. (13) as The overlap between any two operators O 1 and O 2 is described by their characteristic functions as: A major detection device in the photonic entanglement based QIP is a photon detector, which often discriminate only zero or non-zero photons. Such a device is mathematically described by a set of measurement operatorŝ whereÎ is an identity operator. Similar to the state, we can also define the characteristic function of the measurement operator as For example the probability of obtaining the measurement outcome i = 1 (i.e. detector click) by measuring a single-mode stateρ with χ ρ (x) = exp[− 1 4 x T γx] is given by Of course, the above measurement is an ideal one, i.e. unit efficiency and no dark counts. Note that the detector loss can be modeled by a linear loss channel (see the next subsection) followed by a lossless detector. The dark counts are included by replacing |0 0| with a thermal state, which is also a Gaussian state and thus the measurement operator is also kept as Gaussian.
In the SPDC QIP experiments, we often consider the coincidence photon counts of a multi-mode quantum state. Letρ γ be a density matrix of an m-mode Gaussian state with covariance matrix γ. Then the following formula is useful to calculate the coincidence counts:

E. Imperfections
Linear loss. The optical channel with transmittance t (i.e. 1 − t loss) is modeled by combining the channel with a vacuum environment with a beam splitter of cos 2 θ = t and then trace out the environment. Note that the lossy channel transforms Gaussian states to other Gaussian states and thus is a (non-unitary) Gaussian operation. As an example, consider an application of lossy channel with transmittance t to a single-mode Gaussian state with covariance matrix γ. Let A and E represent the channel and an environment, respectively. Then the covariance matrix of the output stateγ A is given by the submatrix of the covariance matrix 1: (a) Temporal mode mismatch at the beam splitter, and (b) its model based on effective beam splitters. Ai and Bi represent the modes where A and B corresponds to the spatial mode and the subscript are the labels for the temporal modes. Note that our model in (b) is phenomenological and thus is not restricted to the temporal mode mismatch.
where S t AE = S BS AE (θ t ) with θ t = cos −1 √ t. Mode mismatch. Imperfection of mode matching at a beam splitter degrades the interference visibility. Figure 1(a) and (b) depict an example of the mode mismatch (in temporal mode) and its theoretical model, respectively. Let ξ be a phenomenological mode match factor with 0 ≤ ξ ≤ 1 where ξ = 1 represents a perfect mode matching. The effect of mode mismatch is precisely modeled by adding two (effective) beam splitters with transmittance ξ as shown in Fig.1(b).
Suppose we would like to measure the coincidence counts after interfering the states via a beam splitter with transmittance t. Let γ A1B1 be the covariance matrix of the two-mode state before the beam splitter. Then after the beam splitting, a whole state is described by a six-mode covariance matrix where and the matrices in the rhs' are the beam splitter matrices defined in Eq. (13). The coincidence count probability is then given by 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.

III. APPLICATION TO THE SPDC BASED QIP EXPERIMENTS
In this section we apply the characteristic function approach to some examples of the SPDC based QIP experiments.

A. Hong-Ou-Mandel interference of an SPDC source
The Hong-Ou-Mandel (HOM) interference test of the SPDC source is a typical experiment to evaluate the quality of SPDC sources [20]. In the test, 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 with the zero-delay point. Let P CC min be the coincidence count (CC) probability without the delay and P CC mean 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. The visibility of the HOM test is then defined by [21] To see the validity of our method, we simulate the HOM experiment in [8] where the pump power dependence of V HOM is experimentally observed with a standard mode-locked laser (Ti:Sapphire laser with a repetition rate of 76MHz). The experimental setup and the corresponding theoretical model are shown in Fig. 2(a) and (b), respectively.
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 effectively modeled by beam splitters with the transmittance t A,B coupled to the environments E A,B , respectively. Moreover, t B also includes the losses in the controllable delay line (the collimator in Fig. 2(a)) which causes an extra loss. η A and η 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.
The covariance matrix of the system ABE A E B before the beam splitters t A,B are given by where I E A E B is a 2-by-2 identity matrix (E A,B are initially in vacua) and "±" means that sign + applies for the first block diagonalized matrix (x quadrature) and − for the second one (p quadrature). The symplectic matrix of the two beam splitters is The state after the losses is obtained by extracting the submatrix for system AB from (S t T γS t ) ABE A E B : .
(25) Application of a 50/50 beam splitter S 50/50 AB with the mode matching factor ξ is then calculated to be where S MM and S BS are given by Eqs. (20) and (21) with t = 1/2, respectively. Finally photon detection consists of lossy channels corresponding to the detector loss and perfect detectors. The lossy channels with η A and η B , quantum efficiencies of the detectors, are applied by operating After discarding environment modes E 1 . . . F 3 , we have the six-mode covariance matrix which we denote γ A1...B3 .
We are now at the position to derive the coincidence count probability: where γ A1A2A3 and γ B1B2B3 are the submatrices of γ A1...B3 . The determinants in Eq. (27) are explicitly given by where µ = sinh 2 r. The coincidence count with delay, P CC mean is simply obtained from P CC min by setting ξ = 0, i.e. no overlap between the two modes.
The theoretical visibilities derived above and the experimental plots in [8] are compared in Fig. 3 where the parameters were experimentally measured as t A = 0.42, t B = 0.29, η A = 0.68, and η B = 0.70. The photon-pair generation rate is related to the squeezing parameter via λ 2 = µ/(1 + µ). The mode matching factor is not measured and thus given as a fitting parameter. The result shows an excellent agreement between the theory and the experiment. Also the theory curves reveal that the vis- ibility is highly sensitive to the mode matching at the beam splitter. Note that [8] 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 [8] and the theoretical curve in Fig. 3 are exactly identical. However, our method is much simpler and more systematic than the approach in [8] and thus 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 η A = η B = η 1 and t A = t B = ξ = 1. Assuming µ 1, which is reasonable for the photon-pair experiments, the visibility is simply given by It should be noted that Eq. (31) agrees with the one obtained in [13]. Another interesting limit is the ideal case (t A = t B = η A = η B = ξ = 1 and µ 1): which shows that the visibility degradation due to the multi-photon emission is accelerated by the low efficiency of the detectors.

B. EPR interference of a Sagnac loop entanglement source
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 [22] and is now widely used for photonic QIP experiments [17,[23][24][25][26].
Here we compare our model with the EPR interference 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 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 angle constant and rotate the other one. The visibility is given by where P CC max and P CC min are the maximum and minimum count rates in the fringe.
The polarizer (HWP and PBS) effectively works as a beam splitter between the horizontal and vertical polarization modes. 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 ≈ 1 and t V ≈ 0 where t H and t V are the transmittances of the horizontal and vertical polarization modes, respectively (see Fig. 5(a) and (b)). In the following, to simplify the calculation of the setup in Fig. 4(b), we useη H = t H η H andη V = t V η V where η H and η V are the quantum efficiencies of the detector for horizontally and vertically polarized photons, respectively. The covariance matrix of the entangled source in Eq. (34) is first transformed by perfect beam splitter The imperfection of the PBSs and the imperfect quantum efficiency of the detectors are included by applying Sη and F V are the environment modes, and discarding these environments. The coincidence count is then given by where Figure 6 shows the comparison with the experiment in [17] and our theoretical model. In the experiment, the overall efficiencies for the horizontal polarization are measured to beη A H =η B H = 0.1. The transmittance of the vertical polarization at the PBS is not measured and thus we use it as a fitting parameter varying from 0.004 to 0.006 (≤ 0.01 is guaranteed for a typical commercial PBS). The figure shows an excellent agreement and also allows us to estimate the leakage of the vertical photons at the PBS.
Again, it is worth to further simplify the closed form given above. By settingη A V =η B V = 0,η A H =η B H 1, and µ 1, we have which agrees with the analyses in [10,13]. In the ideal situation with weak pumping (η A V =η B V = 0,η A H = η B H = 1, and µ 1, the visibility is approximated to be Interestingly, these two expressions coincide up to the first order of µ.

IV. CONCLUSION
In summary, we have proposed a method to compute the SPDC based QIP experiments theoretically, which fully involves the multi-photon emissions and various experimental imperfections. The key ingredient of our method is an application of the characteristic function formalism which has been widely used in CV-QIPs, in particular, for Gaussian states and operations. We apply our methods to two examples, the HOM interference and the EPR interference experiments, and show excellent agreements with the experiments. Moreover, we provide the analytical expression for the HOM and EPR visibilities that include full multi-photon components and various experimental imperfections. Interesting future works include the application of our method to the simulation of more complicated QIPs in which SPDC sources are used, such as multi-partite entanglement generation [3], entanglement swapping [26], QKD [1] and quantum computation [2].