Five-partite entanglement generation between two optical frequency combs in a quasi-periodic χ(2) nonlinear optical crystal

We theoretically prove five-partite entanglement can be produced among modes of two simultaneously generated optical frequency combs via second-order nonlinear interaction in a designed periodically poled lithium niobat (PPLN) crystal. An extendible model is proposed to analyze the entanglement characteristics of generated comb modes by applying van Loock and Furusawa criteria. Our proposal provides a potential approach for generating multipartite entangled states, the so-called cluster states, which are the key resources for quantum computation. Moreover, simultaneously generation of two entangled combs can provide much higher efficiency to generate cluster states.

Five-partite entanglement generation between two optical frequency combs in a quasi-periodic χ (2) nonlinear optical crystal Guangqiang He 1,2,3 , Yu Sun 1 , Linxi Hu 1 , Renhui Zhang 1 , Xikun Chen 1 & Jindong Wang 4 We theoretically prove five-partite entanglement can be produced among modes of two simultaneously generated optical frequency combs via second-order nonlinear interaction in a designed periodically poled lithium niobat (PPLN) crystal. An extendible model is proposed to analyze the entanglement characteristics of generated comb modes by applying van Loock and Furusawa criteria. Our proposal provides a potential approach for generating multipartite entangled states, the so-called cluster states, which are the key resources for quantum computation. Moreover, simultaneously generation of two entangled combs can provide much higher efficiency to generate cluster states.
Quantum entanglement is at the central part of applications such as quantum computation and quantum communication, such as quantum teleportation 1-3 , quantum key distribution 4 , quantum secure direct communication [5][6][7] , quantum machine learning 8 , and so on. To experimentally implement quantum computation, scalability is one of the most essential requirements 9 . Thus a so-called one-way quantum computer model 10 , in which quantum computation resource is provided by cluster states 11 , was proposed for the requirement. A cluster state is a multipartite entangled state in which any algorithm can be implemented by one-particle measurements only. Motivated by this discovery, a large number of works related to cluster states have been carried out 9,[12][13][14][15][16][17][18][19][20] . Various generation methods of cluster states are investigated, such as using beta-barium borate (BBO) 13 , cavity quantum electrodynamics (QED) techniques 15 , and optical frequency combs (OFCs) 9, 18-20 . Among these works, resorting to optical frequency combs is quite appealing. Some interesting works have been proposed about entanglement generated in optical combs. it is shown that optical-frequency combs are formed by the interaction between a cavity mode and a continuous-wave two-tone driving laser consisting of a pump field and a seed field via quantum dot-induced strong nonlinearity 21 . A large-scale quantum entanglement between two comb modes has very recently been explored in an interacting semiconductor quantum dot-photonic molecule system 22 .
OFCs have already found its way to applications such as precision spectroscopy, frequency transfer, astronomical spectral calibration, and generation of low-phase-noise microwave and radio frequency oscillators 23 . Initially, optical frequency combs were produced by mode-locked femtosecond laser 23 , it is stable but bulky and complex. For the sake of miniaturization, microresonators have been proposed and experimentally demonstrated to generate optical frequency combs based on the cascaded four-wave mixing (FWM), allowing considerable reduction of complexity, size and power consumption 24,25 . FWM is third-order nonlinear effect and its conversion efficiency is much less than second-order nonlinear effects such as second-harmonic generation (SHG) and sum-frequency generation (SFG). Thus optical frequency combs generated via second-order nonlinearity have drawn growing attention. Iolanda et al. 26 and Ville et al. 27 have shown that two OFCs can be produced simultaneously by cascaded SHG and SFG. In their works, different features and generation regimes are investigated, good spectral quality is also shown. Inspired by their works, we proposed a novel model to analyze the entanglement characteristics of the two combs in this paper. Simultaneously generation of two entangled combs can provide much higher efficiency to generate cluster states, making an important step for quantum computation. Without loss of generality, we simplify the two combs with three modes in each comb, thus there are six modes in total.
The rest of this paper is arranged as follows. In result and discussion, system model and output fluctuation spectra are given as two parts. First, we describe the physical model for the generation of five-partite entanglement, then derive the main values and fluctuations of the output fields in system model. In output fluctuation spectra, the entanglement characteristics of the modes in two frequency combs are investigated. In method, we introduce the principle of designing the model and analyzing the four-partite entanglement.

Results and Discussion
System Model. In this model, two output frequency combs consisted of six modes are generated by cascaded second-order nonlinear processes as shown in Fig. 1. First, a pump with frequency ω 0 generates the beam with frequency ω 3 through a second-harmonic generation. Then, two beams with frequency ω 1 and ω 2 are generated by a down-conversion process. Last, two beams with frequency ω 4 and ω 5 are generated by two sum-frequency processes.
According to the above nonlinear processes, The energy conversion and phase-matching conditions can be written as follows: k i (i = 0, 1, 2, 3, 4, 5) are the corresponding wave vectors of the six output modes with frequency ω i (i = 0, 1, 2, 3, 4, 5). G 1 , G 2 , G 3 , G 4 are the four reciprocals needed to compensate the phase mismatching. The quasiphase-matching sketch is plotted in Fig. 2.
To obtain the reciprocal-lattice vector for corresponding quasiphase mismatch, an optical superlattice (OS) is needed, and the superlattice can be designed using a dual-grid method. Due to Iolanda et al.'s work 26 , we set the six wavelengths with frequency ω i (i = 0, …, 5) at 1064.45, 1065.54, 1063.36, 532.225, 532.50, 531.95 nm, respectively, and the temperature is 39.5 °C. The tiling vectors of the OS are calculated to be 1.739, 1.740, 1.735, Figure 1. Second-harmonic generation (SHG) with cascaded down-conversion (PDC) process gives rise to the frequency components ω 3 , ω 1 and ω 2 , which in turns leads to two sum-frequency processes (SFG1, SFG2) and generate another two subharmonic components ω 4 , ω 5 . Since the structure of the PPLN crystal is obtained, we design a schematic diagram of physical system as shown in Fig. 3. The system is pumped by a continuous-wave pump laser at 1064.45 nm, amplified by a Yb:fiber amplifier. Then the pump beam enters the cavity through a plane coupling mirror, while other mirrors are high-reflectivity. The high-reflectivity plane mirror is mounted on a piezoelectric actuator (PZT) for cavity length control and the designed PPLN crystal is placed between the two curved mirrors. The nonlinear processes give rise to the output frequency combs, which will finally be separated by an arrayed waveguide grating (AWG) and analyzed by the Fabry-Perot (FP) analysis cavities.
Due to the nonlinear processes proposed above, the interaction Hamiltonian H i for the nonlinear processes, the Hamiltonian for the pump beam H p are as follow: where κ i (i = 1, 2, 3, 4) are the dimensionless nonlinear coupling coefficients of the nonlinear processes, â i are annihilation operators of the modes with frequency ω i . An optical oscillator is an open system since it not only exhibits intrinsic scattering loss with a photon decay rate of γ k0 (for mode k), but also couples waves to the coupling waveguide with an external coupling rate of γ kc . In order to describe such an open system, we present the loss and out-coupling terms as: where ρˆ stands for the density matrix of system and γ k = γ kc + γ k0 represents the damping rate of the loaded cavity. Then the output field is determined by the well-known input-output relation given as ref. 28 in which b is the boson annihilation operator for the bath field outside the cavity. As for the system model presented previously, whole procedure could be governed by the following master equation: The above master equation can be converted into the equivalent c-number Fockker-Planck equation in P representation, which may be written as a completely equivalent stochastic differential equation: where where d is given by  In order to obtain the steady-state solutions of the above processes, the noise terms and all the fluctuations can be neglected; thus the equations for the mean values of the fields can be written as: Though the pump threshold can be obtained in analytical solution, it is too complex. So we choose to get the pump threshold by a numerical method. For example, when γ 0 = 0.01, γ 1 = 0.01, γ 2 = 0.01, γ 3 = 0.01, γ 4 = 0.01, γ 5 = 0.01, κ 1 = 0.03, κ 2 = 0.03, κ 3 = 0.01, κ 4 = 0.01, the threshold is ε th = 0.003243. Notice that, when the pump wave power is below the threshold, there would be no steady solution for output waves. Thus we only investigate the entanglement characteristics above the threshold. Since the mean values of the fields in the cavity are obtained, they can be used to linearize the classical motion equations for the fields in the cavity to obtain the equations of the fluctuations of the fields: in which δα δα δα δα δα δα δα δα δα δα δα δα δα  where m 1 and m 2 is  For the validity of linearised quantum-fluctuation analysis, the quantum-fluctuation must be small enough compared with mean values. If the requirement that the real parts of the eigenvalues of −M stay non-negative is satisfied, the fluctuation equations will describe an Ornstein-Uhlenbeck process 29 , for which the intracavity spectral correlation matrix is given by We introduce the quadrature operators for each mode in order to discuss the five-partite entanglement: The output fields is determined by the well-known input-output relations Eq. 8. In particular, the spectral variances and covariances have the general form   30 . In our discussion, we consider Fokker-Planck equation in P representation and then analyse the entanglement condition that van Loock and Furusawa criteria are violated simultaneously. By using the above quadrature definitions, the five-partite criteria is given by    in which g k (k = 0, …, 5) are arbitrary real parameters that are used to optimize the violation of these inequalities.
Considering the frequency component of ω 3 is generated by SHG, it is hardly to entangle with other modes, we only investigate entanglement characteristics among other five mode. According to the symmetry between ω 1 and ω 2 , we choose to investigate S 01 , S 12 and S 14 in our rest analysis.
Output Fluctuation Spectra. According to Eqs [17][18][19][20][21][22], the stable solution is completely determined by three parameters: the total damping rate γ, the coupling coefficient κ, and the pumping power ε, which in turn determines the drift matrix M, the diffusion matrix D, and the intracavity spectral correlation matrix S. In addition, we also conclude that the parameter γ c plays a role in the spectral correlation matrices according to Eqs [30 and 31]. In the following work, we will vary these parameters to investigate the entanglement.
From Fig. 4, we can see that when γ c /γ = 0.05, there is no entanglement between any two frequency modes. As we increase the out-coupling coefficient, the three frequencies in the first comb begin to entangle with each other. When γ c /γ = 0.75, the ω 1 and ω 2 , which are respectively in two different combs, begin to entangle. And eventually when we set the portion to the γ c /γ = 1, the degree of entanglement is the largest compared with other case. Thus, we conclude that the entanglement among output modes increases as the radio γ c /γ increases. This can be explained that the higher portion the coupling coefficient takes, the less consumed entangled pairs are wasted in the internal loss.
In order to investigate the effect the pump power brings to the degree of entanglement, we firstly set the γ 0 = 0 which means no intracavity loss in this part of discussion. With our previous discussion, the variance S i as a function of ω/γ is merely determined by the parameter ε/ε th while choose proper κ and γ. We plot the minimum variance throughout the noise power spectrum as a function of the pump power which has been normalized by ε th in Fig. 5. We plot the minimum variance versus frequency under different pumping power in Fig. 6.
It can be inferred from Fig. 5 that the variance of S 01 and S 14 would first decrease as the pump power increasing and reach their minimum values when ε = 1.5ε th around. Then S 01 and S 14 would ascend with the pump power while S 12 increases as the pump power since the beginning. Considering that S 01 and S 14 are the short slabs of the entanglement model, we conclude that the 1.5ε th is the best pump power in our case. That explains why we choose ε = 1.5ε th in previous investigation.
In conclusion, we propose the theoretical model for the five-entanglement among modes of two optical frequency combs. By solving Fokker-Planck equation in P representation, we analysed the entanglement case where Van Loock and Furusawa criteria are violated. We analytically find that the intensity of entanglement is completed influenced by the ε/ε th , ω/γ, and γ c /γ. The results would offer a new path for the future study for entanglement over optical frequency combs generated via second-order nonlinear interaction.

Method
We design the PPLN crystal using the so-called generalized dual grid method (DGM), which will phase match the four nonlinear processes. In this method, a dual structure, called dual grid, which contains all the topological information required to built the quasi-crystal is constructed and transformed to a quasi-crystal.
We analyze the entanglement case based on the Van Loock and Furusawa criteria. First, we obtain the main equation and transform it into Fokker-Planck equation. After we get the steady-solution of the above equations, we can obtain the equations of the fluctuations of the fields. Then, we calculate the spectral variances and covariances which will be applied into the Van Loock and Furusawa criteria. Finally, the influences of ε/ε th , ω/γ, and γ c /γ on the intensity of entanglement are analysed.
Our analysis method is also suitable for other types of entanglement, such as multi-partite GHZ state and logic-qubit entanglement, as long as the nonlinear processes are given we can derive the Hamiltonian expressions, then we can analyze the entanglement characteristics according to the method we propose above.