Enhancing entanglement detection of quantum optical frequency combs via stimulated emission

We investigate the performance of a certain nonclassicality identifier, expressed via integrated second-order intensity moments of optical fields, in revealing bipartite entanglement of quantum-optical frequency combs (QOFCs), which are generated in both spontaneous and stimulated parametric down-conversion processes. We show that, by utilizing that nonclassicality identifier, one can well identify the entanglement of the QOFC directly from the experimentally measured intensity moments without invoking any state reconstruction techniques or homodyne detection. Moreover, we demonstrate that the stimulated generation of the QOFC improves the entanglement detection of these fields with the nonclassicality identifier. Additionally, we show that the nonclassicality identifier can be expressed in a factorized form of detectors quantum efficiencies and the number of modes, if the QOFC consists of many copies of the same two-mode twin beam. As an example, we apply the nonclassicality identifier to two specific types of QOFC, where: (i) the QOFC consists of many independent two-mode twin beams with non-overlapped spatial frequency modes, and (ii) the QOFC contains entangled spatial frequency modes which are completely overlapped, i.e., each mode is entangled with all the remaining modes in the system. We show that, in both cases, the nonclassicality identifier can reveal bipartite entanglement of the QOFC including noise, and that it becomes even more sensitive for the stimulated processes.

One of the most striking features of quantum mechanics is quantum entanglement 1,2 , which accounts for the correlations between different parts of a system that cannot be described within the framework of classical physics. The development of the notion of the entanglement has led to the establishment of new branches of physics, e.g., quantum information theory 3 . Apart from its theoretical developments, entanglement has been already experimentally tested and exploited in quantum cryptography 4,5 , quantum communication [6][7][8][9][10] , quantum metrology 11 , quantum information processing 12 , and quantum machine learning 13,14 .
In a real experiment, it is desirable to have a simple, sensitive, and error robust tool for the entanglement identification of a detected QOFC. One of such methods includes a simple measurement of the mean and variance of the field intensity, which can be performed with the help of quadratic detectors and/or spectrometers. Of course, instead of intensities, one has to measure the mean value and variance of the quadratures of the measured fields. Nevertheless, the latter seems more complicated from the experimental point of view, as it requires balanced homodyne detection techniques along with the use of a spatial light modulator that has to shape the spectral profile of a local oscillator. Thus, naturally, one would prefer to resort to some nonclassicality identifiers (NIs) that just include the first-and second-order intensity moments of the measured fields.
For two-mode Gaussian states generated via spontaneous parametric processes, it has been recently shown in ref. 43 , that with the help of stimulated emission and by applying a certain NI, one can conclusively identify the entanglement of such states, by measuring their intensity moments up to the second order. We note, that a recent study in ref. 44 has shown that with the measured variances of displaced quadratures one can reveal nonclassicality of any CV state.
In this work, motivated by the results in ref. 43 , we study a certain NI, which is based on the integrated second-order intensity moments, to show its applicability in identifying bipartite entanglement of the QOFC, which is generated in both spontaneous and stimulation parametric down-conversion processes. In particular, we consider two different scenarios: First, the QOFC consists of many independent two-mode twin beams, i.e., beams with non-overlapped spatial frequency modes. In the second scenario, the QOFC contains completely overlapped entangled spatial frequency modes, i.e., each mode is entangled with all the rest modes in the system. Most importantly, we show that with the help of the stimulation emission, one can enhance the sensitivity of the studied NI in the entanglement detection of the QOFC.
The paper is organized as follows. In Section Theory we briefly review a theory of spontaneous and stimulated down-converted QOFC. There, we also introduce a NI, which is expressed via integrated second-order intensity moments, and which we use throughout the paper. In Section Entanglement identification of QOFC with spatially non-overlapping frequency modes, we apply the NI to the QOFC that has non-overlapping entangled spatial frequency modes, i.e., independent two-mode twin beams. First, we study the performance of the NI for the two-mode case, and later we generalize it to any multimode bipartitions. Additionally, we consider the effect of stimulating fields on the performance of the NI. Section Entanglement identification of QOFC with spatially overlapping frequency modes is devoted to the case when the QOFC contains completely overlapped entangled spatial frequency modes. We show that for such QOFC, the considered NI also proves to be useful for the identification of multimode bipartite entanglement, and that the induced stimulation boosts the performance of a given NI. The conclusions are drawn in Section Conclusions. theory General QOFC generated in spontaneous and stimulated down-conversion processes. First, we briefly review the dynamics of the quantum optical frequency comb generated in spontaneous parametric down-conversion (PDC) process, and that are driven by an intense classical optical frequency comb. Later, we focus on the dynamics of the QOFC that are generated in the stimulated PDC process, where the stimulating fields can also be classical optical frequency combs (COFCs).
The dynamics of the spontaneous PDC process is described by the following Hamiltonian in the interaction picture 45,46 .
I ps i (2) where − E p is the negative-frequency part of the electromagnetic field operator of the pump COFC field, and + E s ( + E i ) is the positive-frequency part of the electromagnetic field operator of the signal (idler) beam 45 ; χ (2) is the quadratic susceptibility of a nonlinear medium. The integration in Eq. (1) is performed over the medium volume V; h.c. stands for Hermitian conjugate.
In the parametric approximation, the pump field, which generates the pairs of entangled photons, can be treated classicaly. Thus, the operator Ê p , in Eq. (1), becomes a c-number. For the COFC pump field, which propagates along the z-axis, the electric-field amplitude can be written as follows 47 . where ω p and k p are the carrier frequency and wave vector of the pump field, respectively, ω r is the angular frequency difference between the teeth of the COFC separated by the time interval ΔT = 2π/ω r . The carrier-envelope-offset phase is denoted by Δφ ceo . The field amplitude A n corresponds to the nth tooth of the COFC, i.e., to the nth frequency of the comb spectrum. For details regarding COFCs, we refer the reader to the appropriate literature, e.g., see refs 48,49 . The operators of the electric fields for both signal and idler beams, which also propagate alone the z-axis, are quantized as follows is the amplitude of the electric field per photon, μ is the frequency-dependent refractive index, and V is the quantization volume.
www.nature.com/scientificreports www.nature.com/scientificreports/ Combining now Eqs (1, 2), and (3), and assuming that the ideal phase-matching conditions are satisfied 45,50 , i.e., ω ω ω ω + = + n p r k k s i , and k p = k s + k i , one arrives at the following Hamiltonian is a coupling constant proportional to both amplitude of the nth tooth of the COFC pump A n , and the nonlinear susceptibility χ (2) , and which is responsible for the coupling between the signal and idler modes with wave vectors k s and k i , respectively. In what follows, without loss of generality, we assume that g k k , s i is a real-valued parameter. The Hamiltonian in Eq. (4) describes the dynamics of the generated QOFC.
If we assume that there are N different spatial frequency modes in a QOFC, then, one can write down the Heisenberg equations for the boson operators, in Eq. (4), as follows where we denoted the matrix exponential as S. Since we consider a system with a finite number of modes, the matrix S can always be presented as a 2N × 2N-dimensional matrix following the Jordan decomposition of the matrix M.
The knowledge of the quantum fields of the QOFC in Eq. (6) allows one to completely characterize QOFC state through the normally-ordered covariance matrix (CM) A  , which for an N-mode state is written as 51 : where A k and A jl are the block 2 × 2 matrices: , : : , To include quantum thermal noise in the system, we employ a standard model based on the superpositions of the signal and noise 46,52 , where the inclusion of this kind of noise, with the mean noise photon-number 〈n〉, affects only the parameters B k in Eq. (7), as B k → B k + 〈n k 〉, and it leaves the other elements of the A  unchanged.
In the case of stimulated PDC, i.e., when the generated QOFC is additionally seeded by stimulating coherent fields, the dynamics of the stimulating fields obeys the same equation of motion as in Eq. (6) for the boson operators, namely: 2 is a complex vector of N stimulating coherent fields, and the matrix S is given in Eq. (6).
With the knowledge of the covariance matrix A  and the vector of stimulating coherent fields Ξ(t), one can easily arrive at the generating function  G , as follows 51 : . The generating function G  allows one to obtain statistical moments of integrated intensities of the fields and also their photon-number probability distribution function. The integrated-intensity moments are obtained from 53 : www.nature.com/scientificreports www.nature.com/scientificreports/ Nonclassicality identifier expressed via intensity moments. One can write down various NIs, expressed in terms of integrated intensity moments of the first and second order. Such an NI can be derived either from a moments-matrix approach or, e.g., from a majorization theory 54 . At the same time, as recent studies indicate, the moments-matrix approach enables finding NIs with better performance than those derived from the majorization theory 38 . Below, we present two possible NIs based on second-order intensity moments, obtained from the moments-matrix approach, for the entanglement identification of bipartite states, i.e., the entanglement between two parts, denoted as signal and idler arms, as follows i  are given in Eq. (12). Whenever E 1 , E 2 < 0, a bipartite state is considered to be nonclassical, particularly, can be entangled.
One of the most important properties of the NIs, E 1 , and E 2 , is that, the quantum detection efficiencies η s and η i of the signal-and idler-beam detectors, respectively, are factorized, i.e., where on the r.h.s. of Eq. (15), the NIs E j are for the ideal case η s = η i = 1. Therefore, in what follows, without loss of generality, we always assume that one uses ideal detectors. As it has been already shown in ref. 51 , the NI E 1 can be used for complete identification of nonclassicality of any mixed two-mode Gaussian state by utilizing a certain coherent displacement to the state. The NI E 2 , as our preliminary analysis has shown, does not possesses this universality. Nevertheless, for particular cases, such as multimode entangled Gaussian states, the NI E 2 can be even more better than the NI E 1 . For instance, the NI E 2 has a much simpler dependence on the number of modes, that can be utilized in more effective state reconstruction techniques based on this NI. Moreover, in the next sections, we show that for multimode QOFCs with either overlapping or non-overlapping spatial frequency modes, the NI E 2 demonstrates a good performance in revealing of bipartite entanglement. Additionally, we show that the stimulation of a noiseless or low noisy QOFC boosts the performance of the NI E 2 . In other words, the NI E 2 increases its negativity with the increasing intensity of stimulating fields. Hereafter, we denote NI E 2 simply as E.
In general, the NI E can describe the nonclassicality of a bipartite state only qualitatively not in quantitative way. In order to relate this qualitative character of the NI E to a quantitative measure, we employ the method used in ref. 54 . Namely, we establish the correspondence between the NI E and the Lee's nonclassicality depth τ, which is a good measure of nonclassicality 55 . The operational meaning of τ is that it determines the amount of thermal noise, one has to add into both arms of a bipartite system, to remove its nonclassicality. When considering a two-mode state, we relate τ to the least negative eigenvalue, taken with opposite sign, of the two-mode covariance matrix A  56 . In this case, the condition τ > 0 is both necessary and sufficient for the nonclassicality of the two-mode state. For the multimode case, when considering multimode bipartitions involving M modes, we refer to τ M as the τ-parametrized NI, τ E M , that can be written as follows 33,54 : which determines the amount of thermal noise τ M , that one also has to add into both signal and idler arms of a bipartite M-mode state to erase the negativity of E M . In other words, the amount of nonclassicality τ M is defined from the condition = τ E 0 M . Importantly, in this case, the condition τ M > 0 is no more necessary but only sufficient for the nonclassicality of a bipartite QOFC state. Since τ M > 0 holds only when E M < 0, according to Eq. (16). But the condition E M < 0 is sufficient but not necessary for the nonclassicality identification. The reason why we resort to the τ M , derived from Eq. (16), and not from the multimode covariance matrix  A , is that for a large number of modes, M ≫ 1, the problem of finding the eigenvalues of a large-size matrix becomes computationally hard. Nevertheless, τ M can serve as a nonclassicality quantifier for bipartite states of the QOFC.

Entanglement identification of QOFC with spatially non-overlapping frequency modes
In this section, we apply the NI E, given in Eq. (14), to the QOFC that consists of non-overlapping spatial frequency modes, i.e., any spatial frequency mode k s is entangled with only one given mode k i . In other words, this QOFC is comprised of many independent two-mode twin beams. We note that, in general, the down-converted frequency modes constituting QOFC, which are generated by different frequencies of the COFC pump, can overlap. The latter case is considered in the next section. Here, instead, we consider the case when such overlapping does not occur. Such QOFC has already been experimentally realized in ref. 57 , and, for example, (2019) 9:5090 | https://doi.org/10.1038/s41598-019-41545-y www.nature.com/scientificreports www.nature.com/scientificreports/ in cavity-enhanced spontaneous PDC [58][59][60][61][62] . Additionally, to make our analysis simpler, we will first focus on a two-mode case and, then, we will proceed to the multimode scenario.
Two-mode entanglement. Spontaneous PDC. For the QOFC, with non-overlapping spatial frequency modes, which is generated in a spontaneous PDC, the boson operators of the signal and idler modes of a two-mode twin beam, produced by the nth tooth of the COFC pump, according to Eq. (6), acquires the following form In that case, the covariance matrix  A , in Eq. (7), of the whole QOFC is factorized on a set of independent 4 × 4 matrices, each corresponding to a two-mode twin beam. The nonzero elements of a given two-mode covariance matrix read as follows: is the mean photon number of entangled pairs, and 〈n j 〉 is the mean thermal noise photon-number in jth mode.
Combining now together Eqs (18,11) and (12), the NI E, in Eq. (14), can be written as where superscript sp in E sp accounts for spontaneous PDC. The expression in the first bracket, in Eq. (19), is a Fourier determinant of the normally-ordered characteristic function of the two-mode twin beam 53 . Hence, when this determinant is negative, the Glauber-Sudarshan P function, which is the Fourier transform of the normally-ordered characteristic function, fails to be a classical distribution function 45,46 . The latter serves as a definition of the nonclassicality and, therefore, determines the entanglement of the twin beam state. Therefore, whenever a twin beam is entangled, E sp always attains negative values. As such, E sp becomes a genuine NI for the two-mode twin beams. Figure 1 shows the dependence of the NI E sp on the Lee's nonclassicality depth τ. This graph indicates that E sp is a nonclassicality monotone for any mixed two-mode twin beam, i.e., whenever τ > 0, then E sp < 0.
For pure two-mode twin beams, the NI E sp attains a simple form Hence, the more intense is the twin beam the larger is its entanglement and, thus, the greater is the negativity of E sp .
Stimulated PDC. In a stimulated PDC process, the generated twin beam at the output of a nonlinear crystal contains a nonzero coherent part due to the presence of stimulating coherent fields. The stimulation process of the twin beams, generated by a COFC pump, can be realized by another COFC that seeds both signal and idler fields.
The dynamics of such stimulating fields, which stimulate the nth twin beam, as given in Eq. (17), reads according to Eq. (10), as follows Hereafter, for simplicity, we assume that the stimulation process is performed by a seeding COFC that stimulates only the signal field, i.e., ξ i (0) = 0. For pure states, the NI E, then, acquires the following form where the first term accounts for the negativity of E st due to spontaneous emission, and the second term corresponds to stimulated emission. For a given value of τ, E st increases its negative value with the increasing amplitude of stimulating field ξ s (see Fig. 2). This means that, the stronger is the stimulating field ξ s , the more negative is E st . Moreover, as indicated by Eq. (22), E st is independent of the phase of the stimulating field ξ s and it depends solely on the coherent field intensity.
Multimode bipartite entanglement. Now, we apply the NI E, as denoted by E M , for certifying bipartite entanglement of the multimode QOFC, consisting of M independent two-mode twin beams. By performing bipartition of a multimode twin beam such that all the signal modes belong to the signal arm, and all the idler modes to the idler arm, E M , then, can be written as follows a n a n a n a n n n n n n n , 2 for a = s, i, n = 1, …, M. Note that we have assumed a general stimulated PDC process in the derivation of Eq. (24).
If the system is comprised of M copies of the same two-mode twin beam with the same stimulation, one then attains It is clearly seen that the first sum in Eq. (27) is the sum of the Fourier determinants of the normally-ordered characteristic functions of each two-mode twin beam, which is in analogy to Eq. (19). For the symmetric case, E M sp becomes the sum of the nonclassicality monotones of each two-mode twin beam. If the nth two-mode twin beam is entangled, then it contributes to the total negativity of E M sp . Hence, the larger is the number of entangled two-mode twin beams in the system, the larger is the negativity of  Fig. 3(a). Therefore, the larger is the spectral energy of the QOFC, i.e., the larger is the number of the two-mode twin beams, the larger is the nonclassicality depth τ M , and, as a result, the larger is the negativity of E M sp .
Stimulated PDC. Now, we consider stimulated PDC, when each nth signal beam is stimulated in the signal arm by a coherent field ξ s,n . Then, the NI E M , in Eq. (23), for a bipartite M-mode twin beam state is Meaning that, in this case, stimulating fields improve the performance of E M st . As in the two-mode case, this stimulation leads to the enhancement of the NI E M st (see Fig. 4). At the same time, as Fig. 4 shows, E M st becomes very sensitive to noise. Namely, by increasing the intensities of the stimulating fields, E M st becomes more negative, but at the expense of losing the tolerance to larger noise. We note that, although the application of the NI E M in Eq. (23) to the multimode twin beam seems straightforward, in practical situations, to separate the signal and idler modes might be difficult. Thus, the following problem arises: How to perform an appropriate bipartition that E M can detect conclusively the modes entanglement. In this case, one needs to implement all possible bipartitions for E M to reveal the maximal total entanglement of the multimode twin beam state.

Entanglement identification of QOFC with spatially overlapping frequency modes
In this section, we discuss another type of a QOFC, namely, when the signal mode of a twin beam generated by the nth tooth of the COFC pump spatially overlaps with the signal or idler modes of the other twin beams produced by different or the same OFC teeth. As a result, one cannot simply divide such QOFC into a set of independent two-mode twin beams, as it was the case discussed in Section Entanglement identification of QOFC with spatially non-overlapping frequency modes. Now, we consider the following interaction Hamiltonian where we assume that the coupling strength g for each generated entangled pair is the same and real. As Eq. (30) implies, any spatial frequency mode k s is equally coupled to various spatial frequency modes k i . Meaning that a given k s mode can contain photons that are simultaneously entangled to different modes k i . For this case, when the Hamiltonian in Eq. (30) contains N different spatial frequency modes, the evolution 2N × 2N matrix, in Eq. (5), takes the form = ⊗ M gL L and L 2 is a N × N hollow matrix of ones, i.e., all its elements equal one, except the main diagonal elements which are zero. The elements of the symmetric exponential matrix S = exp(iMt), in Eq. (6), after straightforward but some algebra, can be found as follows for j = 1, …, N, and k = 1, …, N − 1. Having the matrix S, we can immediately obtain the normally-ordered covariance matrix A  in Eq. (7). Thus, by combining Eqs (7) and (32), we obtain the elements of the matrix A  , which read as follows  www.nature.com/scientificreports www.nature.com/scientificreports/ for j, k = 1, …, N. Since the parameter B p,j does not account for mean photon-numbers of pairs anymore, as it was in Section Entanglement identification of QOFC with spatially non-overlapping frequency modes, we will call it simply as the mean photon number of vacuum fluctuations of the spatial frequency mode j. Considering the stimulation process, where each frequency mode of the QOFC is seeded by a coherent field ξ j , the dynamics of these stimulating fields obeys Eq. (10) with the matrix S, given in Eq. (32).
Two-mode entanglement. Spontaneous PDC. For two-mode entanglement of the QOFC generated in spontaneous PDC with spatially overlapped frequency modes, the NI E, after applying Eq. (14), reads as follows where B j = B p,j + 〈n j 〉 with the mean thermal noise photon-number 〈n j 〉 in the jth mode, and B p,j , C j , D jk , and D jk are given in Eq. (33). For simplicity, we will drop the subscripts in Eq. (34), as we are interested only in two-mode states.
For noiseless QOFC, there is one-to-one correspondence between the NI E sp and the Lee's nonclassicality depth τ [see Fig. 5(a)]. The latter means that with the increasing entanglement between two given modes of the QOFC, the negativity of E sp also increases.
In general, E sp in Eq. (34) fails to detect entanglement between two different spatial frequency modes for noisy QOFC. Namely, as Fig. 5(b) indicates, there is a region of nonclassicality and entanglement, where the NI E sp is positive. Nevertheless, as our numerical findings show, for a two-mode state with large nonclassicality, i.e., with large values of τ, the NI E sp is always a monotone of τ. Moreover, for large number of modes, N ≫ 1, generated in QOFC, there is a bound for τ. Namely, whenever τ > 0.5/N, the NI E sp is always negative [see Fig. 5(b)]. In other words, with the increasing N number of the modes in the QOFC, the NI E sp tends to become a genuine monotone of τ. We note that the value τ = 0.5 is a maximal value of the nonclassicality depth, which can be reached by a Gaussian state.  Figure 5. (a) Nonclassicality identifier E sp versus nonclassicality depth τ for any two spatial frequency modes of the noiseless QOFC without stimulation. (b) The same as in panel (a) but for the mixed two-mode state. The total number of spatial frequency modes of the generated QOFC is N = 100, and the mean photon number of vacuum fluctuations and thermal noise photon number in each spatial frequency mode is B p,j , 〈n j 〉 ∈ [0, 1], respectively. In general, E sp is not a monotone of τ. Nevertheless, whenever τ > 0.5/N, E sp is always a monotone of τ. (2019) 9:5090 | https://doi.org/10.1038/s41598-019-41545-y www.nature.com/scientificreports www.nature.com/scientificreports/ For noiseless QOFC, even when seeding either the kth mode that does not belong to a given two-mode state, the negativity of the NI E st increases with the increasing seeding field ξ k (see Fig. 6). Moreover, E st is independent on the phase of the stimulating signal field ξ k .

Multimode bipartite entanglement.
For a bipartite state that contains M modes in both signal and idler arms, the applied E M , given in Eq. (14), takes the following form We note that, compared to Eqs (23, 36) has a more complicated form due to the simultaneous presence of autoand cross-correlations in both arms denoted as W 1 and W 2 . Since each of those arms contains M modes which are also entangled. In other words, each term in Eq. (36) consists of a sum of different single-and two-mode integrated intensity moments, given in Eq. (35). For a symmetric system, i.e., when all the modes are statistically equivalent, the terms in Eq. (36) can be simplified as follows Spontaneous PDC. In the case of the spontaneous PDC, the negativity of the NI E M sp is increasing with the increasing number M of the modes involved in a given bipartition (see Fig. 7). This means, that by inserting another pair of the spatial frequency modes into the bipartition, one boosts the sensitivity of E M sp in the entanglement detection of a given state.
Stimulated PDC. For a stimulated QOFC, the NI E M st again enhances its sensitivity to detect bipartite entanglement (see Fig. 8). But for larger stimulating fields, the NI E M st becomes less resistant to noise (see Fig. 8). Note that, as in the two-mode case, in order to boost the performance of E M st , it is not necessary to stimulate the measured fields. It is already enough to seed only one of all the N modes of the QOFC, which does not belong to a given bipartition, in order to make E M st more negative. www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusions
In this study, we have shown the usefulness of the nonclassicality identifier, given in Eq. (14), to detect the bipartite entanglement of the QOFC generated in both spontaneous PDC and stimulated PDC processes. This NI is expressed via integrated second-order intensity moments of the detected optical fields which makes it a convenient and powerful tool for the experimental detection of the entangled modes in QOFCs. We have considered two different cases where a QOFC was comprised either by spatially non-overlapping or completely overlapping frequency modes. We have demonstrated that in both cases the NI displays a good performance in revealing bipartite entanglement for noisy QOFC. Most importantly, with the help of strong stimulating fields, one can sufficiently increase the efficiency of a given NI to reveal the entanglement of QOFCs, but at the expense of a higher sensitivity to thermal noise. st versus nonclassicality depth τ M for a bipartition where each part contains three spatial frequency modes for the case of the stimulated noiseless QOFC. The stimulation occurs only in one spatial frequency mode that does not belong to a given bipartition. The amplitude ξ of the stimulating field is: ξ = 0 (yellow solid curve), ξ = 10 (green dashed curve), ξ = 50 (blue dotted curve), ξ = 100 (violet dash-dotted curve). The total number of spatial frequency modes of the generated QOFC is N = 100, and the mean photon number of vacuum fluctuations in each mode is B p,j ∈ [0, 1].