The preparation of non-Gaussian entangled states has attracted considerable interest over the past two decades, as they are critical components in realizing continuous-variable universal quantum computation [1–3] and have demonstrated exceptional performance in a variety of quantum information protocols, including quantum teleportation [4, 5], quantum cloning [6], quantum key distribution [7] and quantum metrology [8]. One of the standard approaches for creating non-Gaussian entangled states is to apply non-Gaussian operations such as photon addition or subtraction to Gaussian states [9–13], and the resultant non-Gaussian entanglement can be witnessed by the negativity of the Wigner function [14]. Another common method for generating non-Gaussian entanglement is higher-order spontaneous parametric down-conversion (SPDC), for instance, the 3rd-order SPDC, in which a pump photon is converted into triplet photons with phase matching condition. The main advantage of the latter over the former is the certainty of the preparation process.

According to whether the triplet photons are degenerate or not, they can be divided into three categories: fully degenerate [15], partially degenerate [16] and fully nondegenerate [17] triple-photon states (TPS). For the fully nondegenerate case, if one or more bright coherent states are seeded into TPS, the resulting state will degenerate into a Gaussian state [18]. It has been shown that the spontaneously generated fully nondegenerate TPS is a Greenberger-Horne-Zeilinger state with super-Gaussian statistics [19], and its entanglement structure has been systematically analyzed by employing different classes of inseparability criteria [20–22]. The concept of nonlinear entanglement was developed in [16] based on the nonlinear correlation features of two modes in partially degenerate TPS (PDTPS). By applying the positive-partial-transposition criterion to the high-order covariance matrices [23], we revealed the mechanism implied by nonlinear entanglement, that is, PDTPS is entangled on a series of 3nth-order covariance matrices and there are competition and coexistence between them.

Recently, three kinds of TPS have been observed in superconducting cavity with ultra-strong nonlinearity [24]. Other potential implementation platforms for these processes include optical fibers [25] and optical waveguides [26]. When light interacts with materials, the high-order nonlinearity in the materials unavoidably induces other nonlinear interactions, such as self-Kerr and cross-Kerr effects [15, 24]. The impact of these effects on the entanglement of PDTPS, however, has not been studied. In this work, we employ the positive-partial-transposition criterion proposed in [23] and the logarithmic negativity to study the influence of these nonlinear effects and the pump brightness on the entanglement of PDTPS. We find that the brighter the pump, the easier it is for the entanglement of PDTPS to transition to higher order covariance matrices. Self-Kerr interactions suppress the partially degenerate 3rd-order SPDC and reduce the threshold of entanglement loaded on the 3rd- and 6th-order covariance matrices. The influence of the cross-Kerr interactions on the entanglement of PDTPS is similar, but it will increase the threshold of entanglement loaded on the 3rd-order covariance matrix.

Let us start our analysis by considering the physical process described by the Hamiltonian [24] where ${\hat{H}}_{I}$ describes the partially degenerate 3rd-order SPDC and κ_0/1/2 is the 3rd-order coupling constant. The annihilation operators $\hat{a}$, $\hat{b}$ and $\hat{p}$ describe the two down-conversion modes and the pump mode, respectively. ${\hat{H}}_{S}$ and ${\hat{H}}_{C}$ are the interaction Hamiltonians describing the self-Kerr and cross-Kerr interactions, respectively.

For the 3rd-order SPDC described by ${\hat{H}}_{I}$, the unperturbed Hamiltonian of the system is ${\hat{H}}_{0}={\hslash }{\omega }_{{ \mathcal P }}{\hat{p}}^{\dagger }\hat{p}+{\hslash }{\omega }_{{ \mathcal A }}{\hat{a}}^{\dagger }\hat{a}+{\hslash }{\omega }_{{ \mathcal B }}{\hat{b}}^{\dagger }\hat{b}$, where ${\omega }_{{ \mathcal A }/{ \mathcal B }/{ \mathcal P }}$ is the frequency of the corresponding mode. Since $[{\hat{H}}_{0},{\hat{H}}_{I}]=0$, the operators ${\hat{U}}_{0}(\theta )=\exp ({\rm{i}}{\hat{H}}_{0}\theta /{\hslash })$ and ${\hat{U}}_{1}(t)=\exp (-{\rm{i}}{\hat{H}}_{I}t/{\hslash })$ are commutative, where θ represents the angle. We define the higher-order quadrature operators ${\hat{q}}_{{ \mathcal A }}^{k}=({\hat{a}}^{k}+{\hat{a}}^{\dagger k})/2$ and ${\hat{p}}_{{ \mathcal A }}^{k}={\rm{i}}({\hat{a}}^{\dagger k}-{\hat{a}}^{k})/2$ $[{\hat{q}}_{{ \mathcal B }}^{l}=({\hat{b}}^{l}+{\hat{b}}^{\dagger l})/2$ and ${\hat{p}}_{{ \mathcal B }}^{l}={\rm{i}}({\hat{b}}^{\dagger l}-{\hat{b}}^{l})/2]$ for Alice's (Bob's) subsystem, satisfying the commutation relations $[{\hat{q}}_{i}^{k},{\hat{p}}_{i}^{k}]={\rm{i}}{\hat{f}}_{i}^{k}$ $(i={ \mathcal A },{ \mathcal B })$, whereWithout loss of generality, we assume that the initial state of the system is ${\hat{\rho }}_{i}={\hat{\rho }}_{\mathrm{th}}^{{ \mathcal A }}\otimes {\hat{\rho }}_{\mathrm{th}}^{{ \mathcal B }}\otimes {\hat{\rho }}_{\mathrm{coh}}^{{ \mathcal P }}(| {\alpha }_{{ \mathcal P }}\rangle \langle {\alpha }_{{ \mathcal P }}| )$, i.e. the pump is the coherent state $| {\alpha }_{{ \mathcal P }}\rangle $ and both modes ${ \mathcal A }$ and ${ \mathcal B }$ are in the thermal state ${\hat{\rho }}_{\mathrm{th}}^{{ \mathcal A }/{ \mathcal B }}={\sum }_{n}(1-{{\rm{e}}}^{-{x}_{{ \mathcal A }/{ \mathcal B }}}){{\rm{e}}}^{-{{nx}}_{{ \mathcal A }/{ \mathcal B }}}| n\rangle \langle n| $, where ${x}_{{ \mathcal A }/{ \mathcal B }}={\hslash }{\omega }_{{ \mathcal A }/{ \mathcal B }}/({k}_{{\rm{B}}}T)$ and k_B is the Boltzmann constant. Using the commutativity of the operators ${\hat{U}}_{0}(\theta )$ and ${\hat{U}}_{1}(t)$, we have where ${\zeta }_{1}=k{\omega }_{{ \mathcal A }}\theta $ and ${\zeta }_{2}=l{\omega }_{{ \mathcal B }}\theta $. Note that the operator ${\hat{U}}_{0}(2\pi /{\omega }_{{ \mathcal P }})$ acting on the initial state ${\hat{\rho }}_{i}$ gives itself, that is, ${\hat{U}}_{0}^{-1}(2\pi /{\omega }_{{ \mathcal P }}){\hat{\rho }}_{i}{\hat{U}}_{0}(2\pi /{\omega }_{{ \mathcal P }})={\hat{\rho }}_{i}$. Substituting this result into equation (3) and after some mathematical calculations, it can be simplified as where ${\vartheta }_{1}=2k\pi {\omega }_{{ \mathcal A }}/{\omega }_{{ \mathcal P }}$ and ${\vartheta }_{2}=2l\pi {\omega }_{{ \mathcal B }}/{\omega }_{{ \mathcal P }}$. Without loss of generality, we adopt the experimental parameters reported in [24], i.e. ${\omega }_{{ \mathcal A }}=6.1$ GHz, ${\omega }_{{ \mathcal B }}=4.2$ GHz and ${\omega }_{{ \mathcal P }}=14.5$ GHz. Obviously, neither $\sin {\vartheta }_{1}$ nor $\sin {\vartheta }_{2}$ is equal to 0. If $(\cos {\vartheta }_{1}-\cos {\vartheta }_{2})=0$, then the condition $2\pi (k{\omega }_{{ \mathcal A }}+l{\omega }_{{ \mathcal B }})/{\omega }_{{ \mathcal P }}=2m\pi $ must be satisfied, where m is a positive integer. In other words, we can always get $\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle $ = 0 if l ≠ 2k. This means that the PDTPS only have quantum correlations of order 3n under this class of observables.

The vector ${\hat{R}}^{k}={\left({\hat{q}}_{{ \mathcal A }}^{k},{\hat{p}}_{{ \mathcal A }}^{k},{\hat{q}}_{{ \mathcal B }}^{2k},{\hat{p}}_{{ \mathcal B }}^{2k}\right)}^{{\rm{T}}}$ composed of high-order quadrature operators allows us to express the generalized commutation relations compactly as $[{\hat{R}}_{i}^{k},{\hat{R}}_{j}^{k}]={\rm{i}}{{\rm{\Omega }}}_{{ij}}^{k}$ where ${{\rm{\Omega }}}^{k}=\mathrm{diag}({J}_{{ \mathcal A }}^{k},{J}_{{ \mathcal B }}^{2k})$ and ${J}_{i}^{n}=\mathrm{adiag}({\hat{f}}_{i}^{n},-{\hat{f}}_{i}^{n})$. The higher-order covariance matrix V^k is defined as ${V}_{{ij}}^{k}=\langle {\rm{\Delta }}{\hat{R}}_{i}^{k}{\rm{\Delta }}{\hat{R}}_{j}^{k}+{\rm{\Delta }}{\hat{R}}_{j}^{k}{\rm{\Delta }}{\hat{R}}_{i}^{k}\rangle /2$, where ${\rm{\Delta }}{\hat{R}}^{k}={\hat{R}}^{k}-\langle {\hat{R}}^{k}\rangle $ and $\langle {\hat{R}}^{k}\rangle =\mathrm{tr}[{\hat{R}}^{k}\hat{\rho }]$. Using the generalized commutation relation and the basic property of PDTPS $\langle {\hat{R}}^{k}\rangle =0$ [23], we can obtain that the states spontaneously generated by ${\hat{H}}_{I}$ satisfy the uncertainty principleApplying the partial transpose operation to Bob's subsystem, the corresponding V^k becomes ${\widetilde{V}}^{k}$, where ${\widetilde{V}}^{k}={{\rm{\Lambda }}}_{{ \mathcal B }}{V}^{k}{{\rm{\Lambda }}}_{{ \mathcal B }}$ and ${{\rm{\Lambda }}}_{{ \mathcal B }}=\mathrm{diag}(1,1,1,-1)$. If modes ${ \mathcal A }$ and ${ \mathcal B }$ are separable, then the partially transposed ${\widetilde{V}}^{k}$ still satisfies the uncertainty principleThe above criterion is a necessary and sufficient condition for the separability of V^k, as has been proved in [23]. Violating it implies entanglement, which can be equivalent to ${\widetilde{\nu }}_{-}^{k}\lt 0$, where ${\widetilde{\nu }}_{-}^{k}$ is the smallest eigenvalue of ${\widetilde{V}}^{k}+{\rm{i}}\langle {{\rm{\Omega }}}^{k}\rangle /2$.

Before studying the influence of the pump brightness and nonlinear interactions on the entanglement of PDTPS, we first introduce them from a statistical point of view. Since the modes $\hat{a}$ and $\hat{b}$ are nonlinearly coupled in ${\hat{H}}_{I}$, there is no analytic solution to the operator evolution equation in the Heisenberg representation. Here, we utilize the Monte Carlo method to numerically solve the master equation $\dot{\hat{\rho }}(t)=-{\rm{i}}[{\hat{H}}_{I},\hat{\rho }(t)]/{\hslash }$ in the Schr$\ddot{o}$ dinger picture [19, 23]. For simplicity, we assume that the initial state of the system is ${\hat{\rho }}_{i}={\hat{\rho }}_{{ \mathcal A }}(| 0\rangle \langle 0| )\otimes {\hat{\rho }}_{{ \mathcal B }}(| 0\rangle \langle 0| )\otimes {\hat{\rho }}_{{ \mathcal P }}(| {\alpha }_{{ \mathcal P }}\rangle \langle {\alpha }_{{ \mathcal P }}| )$.

Figures 1(a) and (b) show the photon number distributions of mode ${ \mathcal A }$ at the linear and log magnitude for interaction strengths ($\xi ={\kappa }_{0}t{\alpha }_{{ \mathcal P }}$) of 0.1 and 0.3, respectively, where modes ${ \mathcal B }$ and ${ \mathcal P }$ have been traced out. For comparison, we also present the photon number distribution of mode ${ \mathcal A }$ in the two-mode squeezed vacuum state generated by ${\hat{H}}_{M}={\rm{i}}{\hslash }{\kappa }_{4}{\hat{a}}^{\dagger }{\hat{b}}^{\dagger }\hat{p}+{\rm{H}}.{\rm{c}}$. We select κ₄ = 2κ₀ because, given the same conditions, the evolution rate of the nonlinear process described by ${\hat{H}}_{I}$ is twice that of the process described by ${\hat{H}}_{M}$. As shown in figure 1(a), mode ${ \mathcal A }$ in the two-mode squeezed vacuum state decays exponentially with increasing photon number (in black), resulting in a straight line in the log scale. The PDTPS, however, exhibits very different traits. With the increase of ξ, the photon number distribution of mode ${ \mathcal A }$ evolves from an arc shape to a ‘∼' shale in the log scale, which is identical to the fully nondegenerate TPS [19]. The photon number distribution of mode ${ \mathcal B }$ in PDTPS is similar to that of mode ${ \mathcal A }$ but on even photon numbers. In figure 1(c), we provide the marginal distribution of mode ${ \mathcal B }$. Unlike the normal distribution of Gaussian states, the PDTPS have higher peaks, longer tails and narrower waists, which are the characteristics of super-Gaussian statistics. Figures 1(d) and (e) show the joint probability distribution of the PDTPS. Starting from the circularly symmetric vacuum state, the PDTPS exhibits extraordinarily complicated correlation features on the standard quadrature after a period of evolution, which is known as nonlinear correlation [16], as opposed to the linear correlation of the two-mode squeezed vacuum.

Different from the 2nd-order SPDC, the amplitude of the pump beam has a non-negligible effect on the 3rd-order SPDC The photon number distribution of the two-mode squeezed vacuum is convergent whether the parametric approximation is used or not. However, for the 3rd-order SPDC, the parametric approximation holds only when ξ is less than a certain threshold [27]. Otherwise, non-physical conclusions, such as photon number divergence, are obtained.

Let us now examine the effect of the brightness of the quantized pump on the entanglement of PDTPS. Figure 2(a) shows the evolution of ${\widetilde{\nu }}_{-}^{1}$ versus ξ for ${\alpha }_{{ \mathcal P }}=\sqrt{10}$, ${\alpha }_{{ \mathcal P }}=\sqrt{25}$ and ${\alpha }_{{ \mathcal P }}=\sqrt{50}$, respectively. The threshold for the existence of non-Gaussian entanglement, which is loaded on the 3rd-order covariance matrix, clearly falls with increasing pump brightness. Reference [16] has shown that the entanglement of PDTPS is enhanced with the increase of the pump brightness. Note that the 3rd-order covariance matrix can only carry part of the correlation information of PDTPS. To fully describe its entanglement properties, we need a series of covariance matrices of order 3n due to its non-Gaussianity [23]. Basically, the larger the amplitude of the pump beam, the easier it is for the non-Gaussian entanglement to jump to higher order covariance matrices. Figure 2(b) demonstrates the evolution of ${\widetilde{\nu }}_{-}^{2}$ with ξ for different ${\alpha }_{{ \mathcal P }}$. The modes ${ \mathcal A }$ and ${ \mathcal B }$ are entangled over the entire parameter interval when the pump brightness is low. However, if ${\alpha }_{P}=\sqrt{50}$, we can plainly see that the entanglement of the PDTPS loaded on the 6th-order covariance matrix vanishes when 0.3 ≤ ξ ≤ 0.5. Next, we investigate the influence of the self-Kerr effect on the entanglement of PDTPS. Theoretically, if there are self-Kerr interactions in the system, then the interaction Hamiltonian is ${\hat{H}}_{I}+{\hat{H}}_{S}$. The evolution operator of the system in the Heisenberg picture can be expressed as ${\hat{U}}_{2}(t)=\exp [-{\rm{i}}({\hat{H}}_{I}+{\hat{H}}_{S})t/{\hslash }]$. Since the operators ${\hat{U}}_{2}(t)$ and ${\hat{U}}_{0}(2\pi /{\omega }_{{ \mathcal P }})$ are commutative ([${\hat{H}}_{I}+{\hat{H}}_{S},{\hat{H}}_{0}]=0$), we can get the conclusion $\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle $ = 0 when l ≠ 2k after similar steps from equations (3) to (4). This suggests that the appearance of the self-Kerr interaction does not change the correlation properties of PDTPS. Figure 3(a) depicts the evolution of ${\widetilde{\nu }}_{-}^{1}$ versus ξ under the conditions κ₁ = 0, κ₁ = κ₀/10, κ₁ = κ₀/5 and κ₁ = κ₀/2, respectively. When ξ is small, the influence of the self-Kerr interactions on the non-Gaussian entanglement loaded on the 3rd-order covariance matrix is negligible. With the increase of ξ, this influence gradually becomes prominent. In comparison to the situation of κ₁ = 0, we can see that the threshold for the existence of non-Gaussian entanglement loaded on the 3rd-order covariance matrix decreases with increasing value of κ₁. Figure 3(b) is the same as figure 3(a), except with ${\widetilde{\nu }}_{-}^{2}$. The evolution of the non-Gaussian entanglement loaded on the covariance matrix of order 6 with ξ for different κ₁ is similar to that of the 3rd-order case. It is worth noting that criterion (6) can only determine whether entanglement exists but cannot quantify it. Here, we quantify the impact of the self-Kerr interactions on the entanglement of PDTPS using the logarithmic negativity ${E}_{N}=\mathrm{ln}\parallel {\varrho }^{{\rm{T}}}{\parallel }_{1}$, where ϱ^T denotes the transpose operation on a subsystem and ∥∥₁ is the trace norm. As shown in figure 3(c), the stronger the self-Kerr interactions, the weaker the entanglement of PDTPS. This shows that the self-Kerr effect not only introduces nonlinear phase shift but also prevents the formation of non-Gaussian entanglement.

We now study the influence of the cross-Kerr interactions on the entanglement of PDTPS. Similarly, in the Heisenberg picture, the evolution operator of the system can be expressed as ${\hat{U}}_{3}(t)=\exp [-{\rm{i}}({\hat{H}}_{I}+{\hat{H}}_{C})t/{\hslash }]$. Since the operators ${\hat{U}}_{3}(t)$ and ${\hat{U}}_{0}(\theta )$ are commutative, the same methods can be used to derive the conclusion $\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle $ = $\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle $ = $\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle =0$ when l ≠ 2k. The entanglement properties of PDTPS are invariant in the presence of cross-Kerr interaction. Figure 4(a) shows the evolution of ${\widetilde{\nu }}_{1}^{1}$ versus ξ for different κ₂ = 0. Surprisingly, the stronger the cross-Kerr effect, the larger the threshold for the existence of entanglement loaded on the 3rd-order covariance matrix, which is obviously contrary to the self-Kerr effect. The non-Gaussian entanglement loaded on the 6th-order covariance matrix, on the other hand, follows a fully different rule, as shown in figure 4(b). Figure 4(c) shows that with the increase of κ₂, the suppressing impact of cross-Kerr interactions on the entanglement of PDTPS becomes more and more evident. In particular, under the same parameters, the cross-Kerr effect is weaker than the self-Kerr effect in suppressing the entanglement of PDTPS.

Finally, we consider the case where both the self- and cross-Kerr interactions are present. Unsurprisingly, the identical result—$\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{q}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{q}}_{{ \mathcal B }}^{l}\rangle =\langle {\hat{p}}_{{ \mathcal A }}^{k}{\hat{p}}_{{ \mathcal B }}^{l}\rangle =0$—can also be reached if l ≠ 2k. This means that the non-zero high-order momenta are already included in the higher order covariance matrices for the class of observables we have chosen. Figures 5(a) and (b) illustrate the evolution of ${\widetilde{\nu }}_{-}^{1}$ and ${\widetilde{\nu }}_{-}^{2}$ versus ξ for different κ₁ and κ₂. If ξ is small, the impact of self- and cross-Kerr effects on the entanglement loaded on the 3rd-order (6th-order) covariance matrix can be ignored. Focusing on the 3rd-order covariance matrix, since the self-Kerr interaction lowers the threshold of entanglement while the cross-Kerr effect increases it, the threshold for the coexistence of these two nonlinear effects is intermediate to that when they exist independently. In particular, these nonlinear interactions reduce the entanglement of PDTPS as κ₁ and κ₂ increase, as shown in figure 5(c).

In summary, we introduced some basic properties of PDTPS and then studied the influences of the pump brightness, self-Kerr and cross-Kerr interactions on their entanglement. We found that the brighter the pump beam, the easier it is for the partially degenerate triple-photon state entanglement to jump to higher-order covariance matrices. The presence of the self-Kerr and cross-Kerr effects alone has exactly opposite effects on the non-Gaussian entanglement loaded on the 3rd-order covariance matrix, with the former lowering the threshold for the presence of entanglement while the latter raising it. The existence of these nonlinear effects does not induce other types of high-order correlations, but reduces the entanglement of PDTPS. These results can further improve our understanding of the process of experimental preparation of PDTPS [24, 28].