Correction: Demonstrating quantum properties of triple photons generated by χ3 processes

. Triple-photon generation (TPG) is based on a third-order nonlinear optical interaction, which is the most direct way to produce pure quantum three-photon states. These states can exhibit three-body quantum correlations, and their statistics cannot be reproduced by any Gaussian statistics of coherent sources or optical parametric twin-photon generator, making them potentially useful for quantum information processing tasks such as quantum state distillation, quantum error-correction and universal quantum computing. Furthermore, the generation of entangled photon pairs heralded by the detection of a third photon can be used in advanced quantum communication protocols. We made the ﬁrst experimental demonstration of TPG in 2004 using a bi-stimulation scheme in a bulk KTP crystal, followed by the quantum theory. The new challenges are now to achieve a spontaneous TPG and the corresponding quantum experiments and protocols using oriented ridge KTP waveguides, which ensures both birefringence phase-matching and light conﬁnement. The waveguides are cut by a precision dicing saw. We recently performed their characterization using third-harmonic generation measurements, which showed their good quality. A rate of about 5 triplets per second is expected when pumping a 5-cm-long waveguide with a 5-W 532 nm beam in the CW regime. Such a spontaneous TPG exhibits low rate of triple photons, which makes the certiﬁcation of quantum features hard. In this article, we review our theoretical and experimental work on TPG and the associated quantum modeling. We also develop theoretical tools for the certiﬁcation of quantum features of spontaneous triple-photon states.


Introduction
Twin photons have deeply influenced the history of nonlinear and quantum optics by their wide range of applications and the paradigmatic place they stand in generating new quantum states of light [1]. As regards triplephoton generation (TPG), the story is only starting. TPG is based on the third-order optical nonlinearity, i.e., the third-order electric susceptibility χ (3) [2,3]: it is a process which can directly generate a 3-photon (3P) state. During TPG, three highly correlated photons at energies ω 1 , ω 2 and ω 3 are created in a nonlinear medium from the annihilation of a higher energy photon at ω 0 , with the energy conservation being fulfilled: ω 0 = ω 1 + ω 2 + ω 3 . Four configurations are possible regarding the level of stimulation of TPG: a stimulated a e-mail: benoit.boulanger@neel.cnrs.fr (corresponding author) TPG over three, two or one modes of the triplet that has to be generated, or the spontaneous TPG, where there is no stimulation at all. Figure 1 shows the three last cases. In the two stimulated cases, the generated photons come from the triplets as well as from the residual photons of the stimulation, as shown in Fig. 1a, b.
At the opposite, the spontaneous scheme shown in Fig. 1c allows to generate a pure 3P state, which corresponds to the third-order spontaneous parametric down conversion (SPDC). These three configurations are all interesting from the quantum point of view regarding the context of continuous or discrete variables. Indeed, these three configurations will all generate states, which exhibit quantum features such as entanglement.
In 2004, we made the first experimental demonstration of a TPG. We considered a two-photon stimulation scheme, as shown in Fig. 1a, using a phase-matched bulk KTiOPO 4 crystal [4]. More recently, we proposed a waveguide approach to boost the generation efficiency The three schemes of generation of a 3-photon state at energies ω1, ω2 and ω3 in a third-order nonlinear medium pumped at ω0: a stimulation over two modes, at ω2 and ω3 for example; b stimulation over one mode, at ω1 for example; c no stimulation. The energy levels are described by continuous lines for the matter and a dashed line for the electromagnetic field [5,6]. TPG has also been an active field of research for many groups around the world [7][8][9][10][11][12][13][14][15]. This new corpus has thus opened new exciting opportunities in quantum optics. Now we wish to overcome a new obstacle by aiming at experimentally generating the 3P state of light by mono-stimulated or spontaneous TPG. This is a real tour de force given both the low efficiency of these two configurations and the fact that the efficiency, contrary to the second-order nonlinear process, increases with the injection intensity. From this step, it will be then possible to open new avenues in quantum information. Indeed, TPG can provide a new powerful resource for advanced quantum information protocols. For example, spontaneous TPG can be considered to generate heralded two photons states, which can be used in a qubit amplifier or in a device-independent quantum key distribution protocol [16]. Our choice for the pump and 3P wavelengths is directly conditioned by this application context: λ 0 = 532 nm, which corresponds to the second harmonic of the Nd:YAG laser, and λ 1 , λ 2 , λ 3 ranging from 1500 nm to 1600 nm, that is to say in the telecom range.
This article is organized as follows: we first review our theoretical and experimental work on TPG in Sect. 2 going from bulk to waveguided configurations, followed by the review of the quantum modeling of TPG in Sect. 3 where we additionally implement an easily accessible numerical solution to this problem that has not known analytical solution. Finally, we propose new tools in order to characterize quantum features of TPG in Sect. 4, a difficult and important problem under active investigation in several international teams.

Theory
We are considering a practical case where the four interacting waves propagate in the same direction. The corresponding wave vectors are then written k i = k i s, i =(0,1,2,3), with k i = 2π λi n(λ i ) where n(λ i ) is the refractive index of the wave at λ i in the considered direction s. The electric fields are taken linearly polarized and are expressed as E i = e i E i (Z)expj(k i Z) where e i is the unit vector, namely the light polarization, E i (Z) the complex amplitude, and Z the spatial coordinate of the laboratory frame along s.
In order to maximize the derivatives ∂E i /∂Z in Eq. (1), which corresponds to the maximization of energy exchange between the four waves, it is necessary to maximize the amplitude of χ (3) eff and to get Δk = 0, i.e., Equation 2 is called the phase-matching relation that ensures a constructive interference between the nonlinear polarization and the radiated field over the full interacting length Z. This condition also corresponds to the full momentum conservation of the photons in the quantum picture.
Then, the design of an optimal TPG requires to find a material with a high effective nonlinear coefficient, low spatial walk-off angles, and allowing phase-matching at the targeted wavelengths. We identified the biaxial crystal KTiOPO 4 (KTP) as the good platform, under two different technologies: a bulk crystal [4] and a ridge optical waveguide crystal [6]. Then, the strategy is to perform a birefringence phase-matching in both cases.
The phase-matching properties are calculated from the refractive indices over the full transparency range of the crystal. For bulk KTP, the best dispersion equations in the visible and near infrared of the principal refractive indices n j (λ), with respect to the dielectric axes, are [19]: The dispersion coefficients are given in Table 1 for λ given in (µm). From Eq. (3) and Table 1, we identified that angle non-critical phase-matching, i.e., ρ 0 = ρ 1 = ρ 2 = ρ 3 = 0, is possible at a pump wavelength λ 0 = 532 nm when the four interacting waves propagate along the x-axis of the dielectric frame (O,x,y,z), i.e., s O x . According to Eq. 2, it means that the principle refractive indices verify = 0. The corresponding phase-matching curve is given in Fig. 2: it shows that the triplet (λ 1 , λ 2 , λ 3 ) is spread over a broad continuum, from 1000 to 4500 nm, that is to say up to the infrared cutoff of KTP. This specific feature of TPG is due to the fact that there are only two coupled relations, the energy and momentum conservations, for the determination of three unknown values. Note that it is completely different from twin photon generation for which there is a unique couple of solutions for the signal and idler wavelengths once the pump wavelength and the direction of propagation are fixed. Actually, in this well-known case of second-order SPDC, there is the same number of equations as solutions to be determined. An alternative to avoid any wavelength spreading, and so any energy spreading, is to fix the value of the wavelength of at least one photon of the triplet: that can be done by stimulating at one or two of the three wavelengths, as in Fig. 1a, b, respectively. Concerning the later scheme, Fig. 2 shows that there exists a partially degenerate scenario where the two injection wavelengths are equal, i.e., λ 2 = λ 3 = 1681 nm for λ 1 = 1449 nm, which can be interesting from the experimental point of view as it will be seen in the next section. Note that it is also obviously possible to stimulate over the three wavelengths.
The same kind of calculations can be done in the case of ridge KTP crystals, knowing that in that case the refractive index has to be replaced by the effective index of the optical modes. For a given dimension of the square waveguide (d × d), the effective index of the fundamental guided mode is calculated using Comsol for x-, y-and z-polarized light, by considering wavelengths varying between 500 and 3000 nm. Dispersion equations giving the effective indices are retrieved by fitting those numerical data assuming a Sellmeier-like form. The best equation we found is expressed as [20]: where the wavelength λ is expressed in nanometer (nm). The dispersion coefficients depend on the transverse dimension of the optical waveguide and are given in Table 2 for d × d = 6 × 6 µm 2 . The phase-matching curves of Fig. 3 are calculated using Eq. (4) and Table 2 . The infrared limit of 3000 nm is fixed by the range of validity of Eq. (4). Figure 3 shows that phase-matching is possible from 1000 to 3000 nm for the three triplet wavelengths. But we can expect to have a wider range, up to 4500 nm, as it is the case in Fig. 2 for bulk KTP, when using proper dispersion equations beyond 3000 nm. Figure 3 also shows that the partially degenerated configuration (λ 1 = 1550 nm, λ 2 = λ 3 = 1619 nm) is allowed in this ridge. It is important to notice that the walk-off angle is nil along the y-axis, which ensures a perfect spatial overlap between the interacting waves.
The phase-matching curves calculated for other values of transverse section (d × d) have all the same shapes and extensions. They mainly differ by the loca-  Fig. 4 where it also appears that the full degeneracy (λ 1 = λ 2 = λ 3 = 1596 nm) is possible for d × d = 6.12 × 6.12 µm 2 . Figure 4 shows well that the geometry of the waveguide is a fine parameter of tunability.
The system described by Eq. 1 can be analytically solved using the sn(u|m) and cn(u|m) Jacobi elliptic functions [21,22]. Given the boundary conditions in terms of energy at the entrance of the nonlinear medium (Z=0), this resolution is only feasible if there is a stimulation at the three wavelengths of the triplet, or at two wavelengths as in the scheme described by Fig.  1a. In this latter case, the boundary conditions are: E 0 (Z = 0) = 0, E 1 (Z = 0) = 0, E 2 (Z = 0) = 0 and E 3 (Z = 0) = 0. Knowing that the intensity is expressed as I i (Z) = n(λi) 2 ε0 μ0 |E i (Z)| 2 , it comes for the four intensities at the exit of the nonlinear medium of length Z = L: with a = Λ 2 The indices x, y and z stand for the dielectric axes of the crystal  Figure 5 gives the corresponding curves for a propagation along the x-axis of a bulk KTP crystal pumped at λ 0 = 532 nm. The point of partial degeneracy shown in Fig. 2, i.e., (λ 1 = 1449 nm, λ 2 = λ 3 = 1681 nm), is taken for the calculation. In that case, the numerical values of the refractive indices are n(λ 0 ) = n y (λ 0 ) = 1.7902, n(λ 1 ) = n z (λ 1 ) = 1.8182, n(λ 2 ) = n y (λ 2 ) = 1.7345 and n(λ 3 ) = n z (λ 3 ) = 1.8128; from Eq. 3 and Table 1, the third-order nonlinear effective coefficient is χ yzyz given in Ref. [23]. The incident intensities that are used are: I 0 (0) = 250 GW/cm 2 and I 2 (0) = I 3 (0) = 3.25 GW/cm 2 . These values correspond to the intensities considered in the experiments described in section 2.1.2. Figure 5 shows well the periodic character of Jacobi elliptic functions, with several values of the crystal length for which there is the full pump depletion, i.e., I 0 (L) = 0. It is then sufficient to consider the first value, i.e., L = 1.27 cm, for fixing the optimal crys-tal length, which is suited for an experiment leading to the extremum I 1 (L = 1.27 cm) = 91.7 GW/cm 2 .
The generation at λ 1 and the amplification at λ 2 and λ 3 are due to the generation of 3P states. Then from the knowledge of I 1 assuming a Gaussian spatial and temporal shapes under the parallel beam assumption, it is easy to access to N 3P that is the number of triple photons, where τ 1 and w 1 are the full-width pulse duration and radius at 1/e 2 , respectively. In the example of bulk KTP, by taking τ 1 = 88 ps and w 1 = 66 µm for example, we obtain N 3P = 2.57 × 10 15 triplets/pulse. When the TPG efficiency is so weak that the pump and stimulation fields can be considered as constant over L, i.e., E 0 (Z) E 0 (0) = 0, E 2 (Z) E 2 (0) = 0 and E 3 (Z) E 3 (0) = 0, the integration of Eqs. 1 is immediate. For Δk = 0, it leads to: The comparison of the behavior of I 1 (L) given in Fig.  5 with that calculated from Eq. 8 using the same boundary conditions is shown in Fig. 6.
It appears that the undepleted pump and stimulation approximation are valid below a crystal length of about 5 mm, but it underestimates the generated intensity at the level of the first maximum of the Jacobi elliptic function and it overestimates the intensity at longer interacting lengths.

Experimental demonstration of bi-stimulated triple-photon generation
An example of experimental setup of a TPG stimulated over two modes is shown in Fig. 7 [22]. The nonlinear medium was a 2-cm-long bulk KTP crystal cut along the x-axis. In order to minimize the number of stimulation beams, we choose the partially degenerated case, i.e., λ 2 = λ 3 , the two corresponding waves being orthogonally polarized.  It was necessary to use very intense pump and stimulation beams, typically several GW/cm 2 because of the weakness of the amplitude of the third-order nonlinearity. This is why we used the picosecond regime. The pump beam at λ 0 =532 nm was the second-harmonic of a 5 Hz picosecond Ekspla SL312-P Nd:YAG laser. The stimulation beam at λ 2 = λ 3 was generated by a homemade optical parametric oscillator emitting at 1665 nm, which exactly corresponds to the experimental phasematching, i.e., (λ 1 = 1474 nm, λ 2 = λ 3 = 1665 nm), that had been determined before from the pioneer experiment thanks to a tunable source for the stimulation beam at λ 2 = λ 3 [4]. Note that these experimental phase-matching wavelengths are very close to the calculated values, i.e., (λ 1 = 1449 nm, λ 2 = λ 3 = 1681 nm) according to Fig. 2, the differences being due to a small inaccuracy of the refractive indices that are used for the calculation. Using a pump intensity I 0 (L = 0) = 250 GW/cm 2 and a stimulation intensities I 2 (L = 0) = I 3 (L = 0) = 3.25 GW/cm 2 , which corresponds to the intensity values taken for plotting the curves of Figs. 5 and 6, we obtained I 1 (L = 2cm) = 0.85 GW/cm 2 at λ 1 = 1474 nm, which corresponds to N 3P = 3.7 × 10 13 triplets/s according to Eq. (7) [22]. The calculation using Fig. 5 at L = 2 cm gives N 3P = 6.16 × 10 13 triplets/s. This small difference with the measurement has been explained by the Kerr effect due to the high intensities that are used [24]. Note also that according to Figs. 5 and 6, a crystal with the optimal length L = 1.27 cm would lead to I 1 (L = 1.27 cm) = 91.7 GW/cm 2 at λ 1 = 1474 nm and so to N 3P = 2.57 × 10 15 triplets/s. The triplets are mixed with residual photons at λ 2 and λ 3 , their numbers corresponding simply to those at the entrance of the KTP crystal, i.e., N 2 (L = 0) = N 3 (L = 0) = 4.25 × 10 14 photons/s, the number of incident pump photons being N 0 (L = 0) = 6.15 × 10 15 photons/s. The use of bulk KTP as described above allows the beams to be strongly confined over a limited length ranging around one centimeter. Actually, it corresponds to the typical value of the Rayleigh length associated with the focusing conditions that are considered. In order to overcome this limitation, which is a crucial point in the case of a spontaneous TPG, we proposed in 2018 to explore the feasibility of a new technology taking advantage of both birefringence phase-matching and confinement. The idea is to use ridge waveguides where the direction of propagation is along a phasematching direction of a KTP crystal. By this way, the pump and triplet waves can exhibit the same propagation modes so that the overlap will be optimal. KTP ridge waveguides are fabricated using a technique based on precise diamond blade dicing [6]. A picture of such a waveguide is shown in Fig. 8. It is then a step index waveguide, the upper and side faces being in contact with air and the lower face being coated with a silica layer. The transverse average section was found to be of about 38 µm 2 , nonconstant along the ridge axis: it corresponds to an average square waveguide of side d = 6.17 µm. We performed a preliminary validation of this technology by achieving high-efficiency third-harmonic generation (THG: ω + ω + ω → 3ω) in the waveguide depicted in Fig. 8 [6]. THG is particularly interesting since it is the exact reverse of TPG that is degenerated in energy, i.e., 3ω → ω + ω + ω. As a consequence, their phase-matching properties are exactly the same. But the advantage of THG is that the conversion efficiency is higher by several orders of magnitude, leading to a much easier way to study phase-matching of TPG. The experiments were carried out using a TOPAS optical parametric generator to deliver the pump beam, with a pulse duration of 15 ps, a repetition rate of 10Hz, and a wavelength that is tunable around 1600 nm. By measuring the third-harmonic (TH) intensity as a function of the fundamental wavelength, we found that the phase-matching was achieved at λ ω = 1594 nm [6]. The calculated value is λ ω = 1594.2 nm from Eq.
(2), with λ 0 = λ 3ω and λ 1 = λ 2 = λ 3 = λ ω , and using Eq. (4) and Table 1. These two values are very close and assuredly within the accuracy of measurement that is at ±1 nm. In these phase-matching conditions, we found that the TH energy conversion efficiency was E3ω Eω = 3.4 % when the fundamental energy is E ω = 2 µJ [6]. The calculation gives 2.6 %, which is a little bit lower than the measurement, but the two values are sufficiently close for a validation of the magnitude of the nonlinear coefficient that is excited here, i.e., χ [6]. Thus, these preliminary measurements of THG allow us to prepare at best the design of spontaneous TPG experiments in KTP ridge waveguides thanks to the validation of the dispersion equations of both the effective index and third-order nonlinear coefficient. The waveguides of the current generation are 3 cm long.

Quantum description
The starting point of the quantum description of triple-photon generation is the interaction Hamiltonian describing the nonlinear process, given bŷ whereâ † 1 ,â † 2 andâ † 3 refer to the creation operators corresponding to the three modes andâ † 0 the creation operator of the pump mode. The nonlinear coefficient κ is proportional to the effective third-order nonlinear susceptibility. The evolution of the quantum system is given in the Heisenberg picture, where the operators follow which reduces to using the Hamiltonian described by Eq. (9). This set of equations is equivalent to the classical ones given in Eq. (1). However, from the quantum mechanics point of view, they have no known analytical solution. We have instead considered numerical methods using the QuTiP package [25,26]. This approach is convenient as long as a small number of photons are considered in order to keep the computation time reasonable, due to the representation of the states and the operators in the Hilbert space of Fock states. To have a significant effect, it is thus necessary to compensate the low pumping field with a higher nonlinear interaction efficiency κ. Indeed, considering the results obtained in [23] and reported in section B, we can infer κ|α 0 | ∼ 1.75×10 −6 1, a rather negligible value. As we are limited in the photon number used in the quantum approach, the expected effects during the interaction described by Eq. (10) will not be observed. In our approach and to gain insights on the quantum properties of TPG, we will instead consider κ = 0.02, corresponding to χ (3) ∼ 1 × 10 −18 m 2 V −2 , a value far above the third-order nonlinearities we can reach in present materials. With this value, a pump field with an average 100 photons is sufficient to observe significant evolution of the quantum system. In the following, we will thus consider TPG evolution under Hamiltonian given by Eq. (9) with a pump field in a coherent state |α 0 , where |α 0 | 2 = n 0 = 100. Such a coherent state has a Poisson distribution and can be expanded in the Fock-basis up to 200 with very good fidelity. As the generated modes are far weaker, we have used a Fock representation with an expansion up to 50, keeping the calculation times reasonable. We start our analysis with the double seeding configuration, where two of the triplet modes are excited with a coherent state containing in average n = |α 0 | 2 = 5. The quantum system is then in the initial state |ψ in = |α 0 , 0, α, α . Next, we consider the single seeding configuration |ψ in = |α 0 , 0, 0, α , and finally the spontaneous triple-photon generation for which the initial state is |ψ in = |α 0 , 0, 0, 0 . Before proceeding further, it is worth noting here that the double stimulation case describes indeed a displacement operator acting on the mode which is initially in vacuum, under the assumption of a classical undepleted pump in the Hamiltonian Eq. (9). In the single stimulation case and under the classical pump approximation too, the Hamiltonian Eq. (9) reduces to the one describing the well-known spontaneous-parametric down conversion. Figure 9 shows the evolution of the different modes computed numerically considering the three initial states. In the travelling wave configuration, the space evolution can be inferred from the time evolution of Eq. (10) by multiplying the time by the speed v of the propagating modes in the nonlinear material. For simplicity, we considered in our analysis that v = 1 and we integrated Eq. (10) up to t = 1. It is very interesting to notice that the amplification of the mode which is initially in vacuum depends on the number of seeded modes. Indeed in the case of the double seeding configuration, at t = 1, the mean photon number for mode 1 is n = 1.76, whereas it is lower, n = 0.27, for modes 1 and 2 in the single seeding regime. These values represent the number of triple photons generated during the interaction. Indeed, we have checked that for the excited modes, the added photon number is exactly Δn = 1.76 at modes 2 and 3 for the double seeded case, and Δn = 0.27 at mode 3 for the single seeded interaction.
A similar analysis with a weaker seeding, |α| 2 = 1, shows that the generated mean photon numbers are even smaller, indicating that the generation rate of triple photons depends not only on the strength of the pumping excitation, but also on the number and strength of the modes excited prior to the interaction. This is unlike the behavior of optical twin-photon para- metric amplification. In the spontaneous configuration, the mean photon number generated is about n 0.04, even lower, and it obviously only depends on the pump intensity and the nonlinear coefficient κ. This is the reason why, despite the fact that more than 10 14 triplets are generated in the double seeded configuration as observed in [22], the spontaneous TPG emission has not yet been reported. Figure 9 shows also the associated quantum fluctuations of the different triple-photon modes. For each mode, we have calculated the variances Δq 2 and Δp 2 of the two quadraturesq = (â +â † ) andp = i(â −â † ). For coherent states and vacuum, which are at the shot noise limit, the variances are Δq 2 = Δp 2 = 1. Our analysis shows that in the case of TPG, the three modes have excess noise and all the variances are greater than unity. Moreover, in the case of the double seeded interaction, the quantum fluctuations are not equally distributed, being larger along theq than thep quadrature. Now, we consider the case of a triple-seeded configuration where the three modes are excited by a coherent state with a mean photon number |α| 2 = 5. In this particular situation, we should also consider the relative phase between the pump and the triple modes given by Δφ = |ψ 1 + ψ 2 + ψ 3 − ψ 0 |. For sake of simplicity, we take φ 0 = 0 and φ 1 = φ 2 = φ 3 = φ. Our analysis reveals two interesting situations: an amplification for φ = π/2, 7π/6 and 11π/6, where the pump photons are converted into triplets; an attenuation for φ = π/6, 5π/6 and 3π/2, which corresponds to a regime where the triple photons are converted back into pump photons. These behaviors are depicted in Fig. 10, representing in the left plot the evolution of mean photon number for φ = π/2 (blue) and for φ = π/6 (red). The inset is the mean photon number at t = 1 as a function of the phase in a polar plot. In the plot on the right, we have reported the quantum fluctuations of the two conjugate quadraturesq andp. The blue curves represent the variances in the case of amplification, and the red ones correspond to the attenuation case. Like in the partially seeded configurations, the quantum fluctuations are always above the shot noise limit. Moreover, as depicted in the contour plot, they exhibit a phase dependence, similar to the one of the mean photon number. The behavior of the fully seeded TPG is comparable to the phase-sensitive twinphoton parametric interaction. It has a dependency to the relative phase between the pump and the generated modes and it exhibits noisy modes, which makes them robust against optical losses. The differences between the twin-photon and triple-photon states are related to the dependence on the seeding intensity and more importantly to the non-Gaussian nature of their statistics. This point will be addressed in the next section.
The generation of triple photons fulfills two conditions: the energy and momentum conservations, as explained in Sect. 2. In the case of spontaneous emission, where only the pump is fixed through its optical frequency ω 0 and its k 0 wave-vector, we end up with two equations with three unknown parameters: ω 1 , ω 2 and ω 3 . The system has thus an infinity of triplet solutions as already mentioned in Sect. 2.
We should thus reconsider our quantum approach by taking into account this broadband generation. A more suitable approach in this case is to consider the space evolution under the nonlinear momentum where where we have replaced the frequency ω 3 of the third photon by ω 3 = ω 0 − ω 1 − ω 2 using the energy conservation condition. In the following, a reasonable assumption is to consider that the pump spectral bandwidth is very narrow in comparison with the bandwidth of the triple photons, especially if a CW laser is used. We can also consider that the pump is strong and undepleted and can be treated as a complex classical amplitude. We obtain the following evolution equations where A 0 is the real amplitude of the pump electric field. Moreover, considering a weak interaction, which is reasonable for spontaneous TPG, we can solve Eq. (14) to the first-order of the Baker-Hausdorff expansion, which giveŝ where ψ = Γ(ω 0 , ω, ω t )LA 0 sinc(ΔkL/2)e +iΔkL/2 . It is now easier to calculate the mean photon number defined as n = ψ in |â †â |ψ in . With the initial state in vacuum, |ψ in = |0, 0, 0 , we obtain: This indicates that the mean photon number at mode i = 1, 2, 3 at frequency ω is the integral over all the contributions of triple photons that fulfill the phasematching condition and the energy conservation. For example, in the KTP waveguides, these two conditions are fulfilled over almost the full transparency window of the nonlinear crystal as shown in Fig. 3. Figure 11 shows the spectral density distribution |ψ(ω 0 , ω, ω t )| 2 as a function of λ 1 and λ 2 . It is obtained for the KTP waveguide considered in Fig. 3 and using the corresponding dispersion relations of the effective index Eq. 4.
Integrating |ψ(ω 0 , ω, ω t )| 2 over λ 2 -axis, we obtain the mean photon number n 1 (ω, L) per second as a function of λ 1 . Similar calculations can be hold to compute n 2 (ω, L) and n 3 (ω, L) . The expected mean photon number is higher than the one calculated by the single mode model described by the Hamiltonian given by Eq. (9). The total mean photon number that we expect is further estimated by integrating Eq. (16) over the full spectrum, i.e., Indeed, considering the full phase-matching bandwidth of the KTP ridge waveguide as shown in Fig. 3, i.e., from 1 to 3 µm, and a CW pump power of 5 W at 532 nm, we estimate N 3P (L = 5cm) 4.9 triplets/s. Such a pump level is a reasonable target regarding the expected improvements of the losses and surface quality of the waveguides. Note that the rate of triplets remains the same in the pulsed regime than in the CW regime if the average power is kept at the same value. But the number of triplets per pulse will depend on both the pulse duration and repetition rate. For example a pump of 5W at 532 nm with a pulse duration of 11.3 ps and a repetition rate of 88 MHz, gives N 3P (L = 5cm) = 4.9 triplets/s, corresponding to 5.6×10 −8 triplets/pulse.

Non-classical features of TPG
We now focus on spontaneous TPG, and we consider the non-degenerate case where the Hamiltonian iŝ and the degenerate case characterized by the following Hamiltonian:Ĥ In both cases, we consider a pump with a high photon number so that it can be treated as a classical state. We also define ξ = κα 0 where |α 0 | is the amplitude of the pump. The generation of TPG quantum states has been performed in the double-seeded configuration in [24], and in the single seeded and spontaneous cases in [15]. In those experiments, the authors show strong evidence that the states are generated by a TPG; however, the experimental demonstration of quantum features for such states remains, up to our knowledge, to be done. In this section, we develop the theoretical tools that allow to demonstrate the different quantum features of the states generated by all the possible configurations of TPG, i.e., stimulation over one, two or three modes, as well as no stimulation. We first note that all quantum features can be inferred from the density matrix of the state or equivalently from a phase space distribution. In the present case, we expect 4.9 triplets per second, which makes a full quantum tomography totally out of reach. We thus focus on tools which require a minimum information about the state but are still able to conclude about quantum features. In the first part, we focus on the states generated by the Hamiltonian (18) and propose tools to demonstrate genuine multipartite entanglement (GME) of such states using homodyne measurement. We then focus on the state generated by the Hamiltonian (19) and build tools in order to demonstrate non-classical features of such states.

Entanglement
A state is said to be GME if it cannot be written as a biseparable state for any bipartition. A general biseparable state is a mixture of states that are product states for some bipartition, i.e., a partition of all modes into two groups. Formally, it is given by: Here, the sum runs over all 3 partitions G 1 |G 2 of the 3 parties where G 1 ∪ G 2 = {1, 2, 3} and G 1 ∩ G 2 = ∅. The probabilities of different partitions sum up to one, G1|G2 p(G 1 |G 2 ) = 1, and ρ G1|G2 is a separable state with respect to the partition G 1 |G 2 . We aim to demonstrate entanglement using homodyne detection, so that we define a general homodyne measurement on the mode i asX θ i = (â † e iθ +âe −iθ )/2 and writeX π/2 i =p i andX 0 i =q i . Our first attempt to describe the quantum properties of triple photons, especially their three-body quantum entanglement, was in 2018 [27] using the existing tools dedicated to characterize the inseparability of multi-body quantum states. Our analysis is based on P. van Loock and A. Furusawa non-separability criterion S [28], relying on the evaluation of the quantum fluctuation of the generalized quadratureŝ where h i and g i are arbitrary real numbers. Thus, the non-separability criterion is easily accessible experimentally using standard balanced homodyne detections. The criterion is defined as S = Δû 2 + Δv 2 . We can show that when S < 2 min(|h k g k | + |h l g l + h m g m |) for any permutation of k, l, m = 1, 2, 3, the triplets exhibit genuine entanglement. However, when S > 2 (|h k g k | + |h l g l | + |h m g m |), the system is completely separable. Surprisingly, we found that the spontaneous triple photons do not exhibit three-body quantum entanglement in the continuous variable (CV) regime, always fulfilling the last inequality. Entanglement has only been predicted in the different seeded cases as depicted in Fig. 12. The different results are obtained following the analysis reported in our work [27]. In fact, TPG is a third-order nonlinear process with non-Gaussian statistics. This statement is obvious in the degenerate case whenâ 1 =â 2 =â 3 and the Hamiltonian (18) becomes (19). It will be shown hereafter in Fig. 16 that the Wigner function exhibits interferences and negativities, which is a clear signature of the non-Gaussian nature of the triple-photon mode a. In the non-degenerate case (Hamiltonian (18)), we have recently shown that even though the Wigner function associated with each mode looks Gaussian, the quadrature probability distribution is super-Gaussian [29]. The S criterion is thus not anymore a relevant parameter to analyze and reveal the entanglement properties of the triple-photon states. Indeed, the S criterion is based  . 12 Evolution of the entanglement criterion S for different interactions TPG configurations. Red: fully seed TPG with a relative phase of π/2. Blue: double-seed TPG. Green: One mode seeding configuration and S calculated for the two unseeded modes. |α| 2 is the mean photon number per second of each seeding beam on the second-order moments, i.e., variances, which are sufficient to describe a Gaussian distribution. When at least one of the triplet modes is seeded, the statistics become again Gaussian, and hence, using Furusawa and Van Loock criterion, S is sufficient. The failure of the S criterion to describe the entanglement in the case of spontaneous TPG due to its super-Gaussian nature pushed us to further explore the entanglement nature of the triplet, in order to claim that the triple photons are indeed entangled despite the results based on S. Hence, we have used the logarithmic negativity defined as E N = ln ||ρ Ti ||, where ρ is the density matrix of the triple photons and ρ Ti is its partial transpose over mode i = 1, 2 or 3 [29,30]. A quantum system is said to be entangled whenever E N > 1. Figure 13 shows the calculated E N in the case of spontaneous TPG, starting from the Hamiltonian Eq. (9) and using a numerical approach to solve the Heisenberg equation. For this simulation, we considered a coherent pump field with a mean photon number n 0 = 10. The logarithmic negativity E N increases with the nonlinear parameter κ|α 0 |, reaching E N 2 for κ|α 0 | = 1, which clearly demonstrates the three-body entanglement of the triple photons. In the same figure, we have also reported the evolution of n 0 and n 3P , respectively, the pump and triplets mean photon numbers.
Even though the logarithmic negativity has the ability to reveal the entanglement of a quantum system, it is very hard to measure in practice, as it requires a full tomography of the state. One instead can use a witness of genuine entanglement, i.e., an observableÔ such that Ô ≤ 0 for all biseparable state. As stated before, the state |ψ 3 associated with non-degenerate TPG has non-Gaussian entanglement, which implies that a witness of entanglement (or GME) for the state |ψ 3 requires at least a third-order field operator. The authors of [31] derive a non-Gaussian witness of GME, which allows the demonstration of genuine entangle- Fig. 13 Red: logarithmic negativity of the spontaneous TPG excited by a coherent pump with a mean photon number |α0| 2 = 10. Black: Evolution of the pump mean photon number. Blue: Evolution of the mean photon number of the generated triplets ment for the state |ψ 3 . Moreover, their witness failed to detect entanglement in the seeded configuration, which highlights the difference between the two states. We briefly review their witness. The authors of [31] show that any biseparable state satisfies | â 1â2â3 | ≤ n 1 n 2n3 + n 2 n 1n3 + n 3 n 2n1 (22) wheren i is the number operator for mode i. The state |ψ 3 does not satisfy (22) and is thus GME. We want to formulate a relaxation of this witness only in terms of local homodyne measurement. In order to formulate our relaxation, we first note that (22) implies that any biseparable state satisfies: By summing (22) and (23) and using the triangular inequality, we end up with which holds for all biseparable states. Interestingly, the left-hand side of this previous equality can be measured using local homodyne measurement , i.e., The termsn inj being hard to measure experimentally, we bound them according to where n 2 i can be measured locally by combining second-and fourth-order moments of phase average field quadrature, that is to say: We define the quantity where . ρT P G(η) is the expectation value of . on the state ρ T P G . We can certify the presence of GME for the state ρ T P G if the quantity is positive. We plot in Fig. 14 w with respect to the efficiency η. We see a violation even for low efficiency, which makes the detection of GME for the state ρ T P G robust against losses when homodyne detection is used.

A hierarchy of quantum features for single mode of light
We focus on the degenerate state |ψ , which we take to be the outcoming state created by the process associated with the Hamiltonian of Eq. (19), and which we assume to be single mode. There exists a hierarchy of quantum features for a single mode of light in Fig. 15. The first layer of quantum features is the nonclassicality, which is defined using the Glauber-Sudarshan (GS) or P distribution [32,33]. Any state ρ can be written in terms of coherent states |α as follows: A state is classical if its GS distribution P (α) can be interpreted as a classical probability distribution [34]. Non-classicality is a necessary feature for a state to exhibit any quantum advantage, for example in quantum metrology where any advantage over classical metrology requires non-classicality. The second layer of quantum features is quantum non-Gaussianity (QNG). A Gaussian state is by definition a state with a Gaussian Wigner function. Hudson's theorem [35] stipulates that all non-Gaussian pure states are Wigner negative; nevertheless, it is not the case for mixed states. Any pure Gaussian state can be generated by the action of the squeezing operator S(s) = eâ 2 s * −â †2 s following by displacement operator Also, any state can be written in terms of Gaussian states, i.e., ρ = p(α, s)|s, α s, α|.
A state is Gaussian if p(α, s) can be interpreted as a probability distribution. Non-Gaussian quantum states are essential to a variety of quantum information processing tasks such as quantum state distillation [36], quantum error-correction [37] or quantum computational speedup [38]. The third layer of quantum features is Wigner negativity. The Wigner function is a representation of a single mode state ρ in terms of the a quasi-probability distribution given as [39]: Wigner negativity is arguably the strongest form of quantum feature for a single mode state and is a necessary condition to perform efficient quantum computing using CV states [40]. We plot in Fig. 16 the Wigner function of state |ψ as a function of its canonical quadratures q and p. We can observe different areas of negativities in the Wigner function, which implies that the TPG state is Wigner negative, non-Gaussian and non-classical. We note that this function is not phase-invariant, and as a consequence the presence of such negativities cannot be detected by only photon counting strategies. The reconstruction of the Wigner function requires a high number of experimental runs, which is not available in the present experiment. Witnesses of Wigner negativity can be systematically derived using hierarchy of semi-definite programs [41], but those witnesses require again a number of measurement runs that is much higher than the one available in our case, and for those reasons we proceed by focusing on non-classicality and non-Gaussianity only.

Witnessing non-classicality of the degenerate state |ψ
A witness of non-classicality can be represented by an observableŴ together with the maximum expectation value ofŴ on a classical state. The most widely used witness of non-classicality is probably the second-order correlation function g 2 (0) [42] together with its classical bound 1. Several criteria that analyze matrices of moments of annihilation and creation operators have also been derived [43][44][45]. In this section, we focus on the setup of Fig. 17 and derive a witness tailored to the state |ψ and that is based on photodetection events. We consider a single-mode case for the sake of clarity. The results are holding in the multimode case and, since i) the detector cannot distinguish between differ-  Fig. 17 does not allow us to acquire information about the coherence terms of |ψ , then we can consider a single-mode state with the same number of photons than that of the expected multimode state in order to correctly model our experiment (see Appendix A).
In Fig. 17, an input state ρ i of the ith experimental run is split into four spatially separated modes using three balanced beam splitters and sent to four non-photon number resolving (NPNR) detectors. We do not want to make assumptions about the efficiency of the NPNR detectors, so we consider perfect detection and that any inefficiencies can be mapped into the input states ρ i . A NPNR detector with unit efficiency can be modeled by a two element positive operatorvalued measure (POVM) as the set {E • , E • } = {1 − |0 0|, |0 0|}, corresponding to the single detector events click and no click, respectively.
The set of non-classical states is convex (29), so that the set of the different probabilities of clicks achieved by classical states is also convex, by linearity of the trace. We can thus consider a linear combination of probability operators corresponding to the events of click and no click on the detector arrangement, whereP i are the POVM elements of the arrangement, which correspond to the event where i detectors click while the others do not click, as for example: • E sin θ 1 sin θ 2 cos θ 1 sin θ 3 cos θ 2 cos θ 1 cos θ 3 cos θ 2 cos θ 1 , The next step is to compute the maximum expectation value of the observable constructed in Eq. (32) that any classical state can achieve, i.e., We note that Tr Ŵ (θ)ρ c is linear in ρ c , and since any arbitrary classical state can be written as a mixture of pure classical states according to Eq. (29), we find that: Therefore, we can conclude that the maximum of the RHS of (34) is achieved by a pure coherent state. With this simplification, we have that We find that this maximization equivalent to a problem of finding the roots of a third-order polynomial (see Appendix A) and can thus be simply performed analytically. The vectorθ = (θ 1 , θ 2 , θ 3 ) is associated with the weights of the different events defined by the opera-torsP i on the witnessŴ . The witness is now arranged such that if for a given state ρ, it exists one vectorθ for which Tr(Ŵ (θ)ρ) > w c (θ), then the non-classicality of the state ρ can be demonstrated using the setup of Fig. 17.
In practice, the state |ψ will experience losses; we then define the state where U η is a unitary corresponding to a beam-splitter with transmitivity η, referring to the efficiency, and the quantity q(θ) = Tr(ρ T P G (η)Ŵ (θ)) − w c (θ), that will be positive when the witness can detect non-classicality for the lossy state ρ T P G (η). In order to find the optimal witness of the state for a given efficiency, we can compute: We plot in Fig. 18 q opt with respect to the overall efficiency η of the setup for = ξt = 0.01. We find that, in the present case, the optimal q is always found when w c (θ) = 0. We also find a positive q opt for the range of all efficiencies, which proves the detection of nonclassicality of the state |ψ using the setup of Fig. 17 to be robust against losses.
In the last part of this section, we turn to give an estimation of the number of experimental runs that would be necessary to estimate the quantity W s = Tr(ρ T P G (η)Ŵ (θ)). We assume that the POVM ele-mentsP i are independent quantities that are measured N times. At each run, we evaluate a random variable X i , associated withP i , which takes the value 1 when i detectors click and i − 4 do not click, and the value 0 otherwise. An unbiased estimator of W s after N runs is given by: We bound the standard deviation ofW s as follows: where σ Xi = X i (1 − X i ) is the standard deviation of a binary variable. The number of runs that is needed to estimate the value of the witness, with a precision 3 times smaller than the distance to the classical bound, can thus be estimated by finding the number of runs N opt such that: We plot N opt with respect to the efficiency in the inset of Fig. 19. With 4.9 triplets/s and a 5-W laser at 88 MHz, we find that 18 seconds of experiment are enough to certify the non-classicality of the state for an efficiency of 50%.

Demonstrating non-Gaussianity of the degenerate state |ψ
A witness operatorŴ in the same general form of Eqs. (32)(33) associated with the linear optical setup of Fig.  17 proves to be useful to the demonstration of the non-Gaussianity of the one-mode state |ψ as well, meaning that, by analyzing the mean value of this operator on this state, we can assert that |ψ cannot be written as a convex mixture of pure Gaussian states. Witnesses of non-Gaussianity for arbitrary Fock states using linear optics have been derived in [46]. The aim of this section is to derive a witness of non-Gaussianity tailored to the state |ψ where the full knowledge of the probabilities distributions of each detector is used, and not only of two as in [46]. In order to witness such non-Gaussianity, we again need to specialize the witness over the free parameters of the coefficients c(θ). A first step is defining and computing the maximum mean value w g of the witness operator over all Gaussian states, i.e., We note one more time that this expression is linear in ρ G , and thus, we can conclude in the same manner that its maximum is achieved by a pure state and that it is then sufficient to perform the optimization over the set of pure Gaussian states. In this case, as opposed to the witnessing of non-classicality via w c (θ), the expression for w g (θ) is more intricated and the optimization was performed numerically over all states of the form i.e., squeezed coherent states, in whichD α andŜ s are the canonical displacement and squeezing operators. More details can be found in Appendix B. The second and last step is performing a further optimization, now over the set of the coefficient parameters θ, of the difference where ρ T P G is the lossy state given by Eq. (37), in order to find g opt = max θ (g(θ)). Finding a positive value of g opt certifies that the mean value of the witness in state ρ T P G breaks the Gaussian bound and thus that the state |ψ cannot be represented by a mixture of states of the form of Eq. (43), meaning that this state is non-Gaussian.

Fig. 19
Optimized mean values of gopt. The optimization is performed numerically over the set of parameter coefficients θ for a range of efficiencies η and for = 10 −2 . In the upper inset, we plot the expected minimum number of runs Nopt necessary for non-Gaussianity certification. In the lower inset, we zoom in on the region of efficiency η from 0 to 45%, showing that we find positive values of gopt for all efficiencies Figure 19 shows the values of g opt obtained by numerical optimizations performed over a range of efficiencies η of the setup. We find that a violation of the Gaussian bound can be witnessed by the auto-correlation witness of the form of Eq. (32) over the range of all efficiencies, which again demonstrates the robustness of the method against photon losses. In order to estimate the number N opt of experimental runs needed for the certification, we one more time set where σ Ws is given by Eq. (39), but now using the coefficients c i optimized for the non-Gaussianity witness. This estimation is plotted in the upper inset Fig. 19 for the range of efficiencies η.

Discussion and conclusion
In this article, we have reviewed our past work, both experimental and theoretical, from bulk crystal to KTP waveguide and from a simplistic quantum analysis extrapolated from second-order nonlinear process, to a deep understanding of peculiarities of quantum behavior introduced by the third-order nonlinear interaction. We thus highlighted these differences, such as the dependence of the classical and quantum behavior on the seeding intensity, and the emergence and constraints induced by the non-Gaussianity. Then, we presented original theoretical results and fixed one of the remaining important problems: how to experimentally conclude at the quantumness and the non-Gaussianity of the generated triple-photon states engaging a minimal resource budget. We have thus introduced quantum witnesses using the minimum number of photocounting detectors, optimizing thus losses and complexity. The spontaneous generation of triple photon statistics by the third-order nonlinear interaction is still elusive, but we feel now closer. The following step is to demonstrate the ability of these states to exhibit stronger form of quantum correlation such as Wigner negativity or non-locality. It is then possible to use these behaviors in quantum information protocols. The ability of generating heralded pairs of photons can be used in quantum repeater protocols. Moreover, using photon detectors preceded by displacement operation on one of the three modes of the state |ψ 3 , it is then possible to herald a state of the form |00 + |11 that has the ability to violate Bell inequalities [47]. Time-bin genuine entanglement and non-locality can be demonstrated on the state |ψ 3 using Franson-type measurement. Indeed, all the photons can be sent to a single imbalanced interferometer and look at threefold coincidences at the outputs of the interferometer. Then, if the arm length difference of the interferometers is smaller than the coherence length of the pump, it is by principle impossible to tell if all three photons have taken the short or the long arm, effectively generating a three-mode Greenberger-Horne-Zeilinger state (GHZ state). The long-term aim is to add TPG process to the quantum optics toolbox for quantum communication and computation.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.]

Declarations
Competing interests Authors declare to disclose interests that are directly or indirectly related to the work submitted for publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Non-classicality witness
A pure coherent state can be written in the Fock or number state basis as where α is a parameter defining the specific state. The action of balanced beam-splitters at coherent states given by Eq. (A1) is straightforwardly computable. Actually, after a beam splitting, the state is represented by two pure coherent states attenuated as α → α/ √ 2 at each output path, so that the probabilities P i of observing i clicks at the detectors of Fig. 17 can be written as for i = 1, 2, 3, 4. The non-classicality witness will then simply be given by the summation with c i (θ) defined by Eq. (33), and the maximization of w c (θ) of Eq. (34) will consist of a maximization over the sole parameter α. If we first define the change of parameters x = exp(− |α| 2 4 ), the value of α which maximizes Eq. (A3) can be found by computing the roots of so that α opt will be given by ± −4 ln(x), with x adequately picked among the found roots of Eq. (A4). As can easily be seen, due to the form of the P i , Eq. (A4) is, in turn, a third-order polynomial in the new variable x, and thus, its roots can be routinely found. By computing the value of W (θ) at values α associated with each of the found roots, and verifying which of those points consist of W (θ) maxima ( ∂ 2 W (θ) ∂ 2 α < 0), the maximization of Eq. (34) can thus be completely computed analytically.
Let us focus on the multimode case, since we cannot have access to the coherence on the setup. A general multimode classical state can be written as: The POVM associated with the event "i detector click and 4-i does not clicks" can be written as: The maximum of Tr(W (θ)ρ m ) is achieved by pure multimode states for which we have Tr(P im |α k α k | ⊗k ) which is the same expression than (A2) if we take |α| = k i=1 |α k | 2 . Thus, the bound for single mode holds for the multimode case.

Appendix B: Non-gaussianity witness
The expectation value of the probability P i of observing i out of 4 detectors clicking is associated with the projector operatorP i of Eq. (32) acting on the state after the unitary evolution U representing the beamsplitters, and for a general pure Gaussian state |G arriving in the setup of Fig. 17 this can be written as: whereD α andŜ s are the canonical displacement and squeezing operators, with displacement and squeezing parameters equal to α and s, respectively, and withP i given byP whereP 0 = |0 0| is the projector onto vacuum at each detector. We use the fact that any Gaussian state |G can be written as a squeezed coherent state and, without loss of generality, we can choose s to be real while keeping α generally complex. Alternatively, and for convenience of calculation, we can reinterpret Eq. (B1) as a reverse evolution U † acting on theP i operators and compute the expectation value of U †P i U in the incoming general pure Gaussian state |G . In accordance with Eq. (B2), this can be performed by deriving all the operators of the form U † (P 0 ) k (Î) 4−k U , each of those associated to having 0 clicks at k detectors without mention of what is observed at the other 4 − k detectors. This can be done as follows.
For a Fock state |n , the probability of detecting a vacuum state (no click) at a detector placed after two balanced beam-splitters, without mention to what is observed at any remaining detectors, can be written as: For any convex combination of Fock states |n , the associated projector onto vacuum at this detector (U † (P 0 ) 1 (Î) 3 U ) can then be written in terms of the number operator a † a as: The operators associated with the detection of two or three vaccum states after the beam-splitters can be derived in the same manner, and are given, respectively, by:P Finally, the operator associated with all four detections not clicking (P 0000,n = U † (P 0 ) 4 U ) will just be the vacuum operator |0 0|, with associated probability being the modulus squared of the overlap of the incoming Gaussian state with the vaccum, that is G|0 0|G = | 0|G | 2 . As derived in ref. [48], we have that: Again, instead of directly computing the expectation value of such projection operators in a general Gaussian state, in accordance with Eq. (B1), we can apply the reversed squeezing operator to these projectors and compute their expectation values on a coherent state |α . Following the result in Eq. (9) of ref. [49], and the fact that e (ka † a) =:e ((e k −1)a † a) :, where :: represents normal ordering, we havê S † s e ka † aŜ s = e − k 2 (e where S s = sinh (s), C s = cosh (s), T s = tanh (s) and F s = e −k C 2 s − e k S 2 s , and, to be consistent with ref. [49] results, we have to take s to −s. By taking values for the parameter k in accordance with each of the projectors of Eqs. (B4) and (B5), we will have all needed terms of the form S † (s)U † (P 0 ) k (Î) 4−k US(s), and we are ready to compute their expectation values on the coherent state |α : as Eq. (B7) is normal ordered in the field operators, it can be directly applied between α| and |α in Eq. (B1) using α|a †2 |α = α * 2 , α|a 2 |α = α 2 and α|a † a|α = |α| 2 .
Moreover, we can optimize the witness by computing g opt = max θ (g(θ)).