Scalable multimode entanglement based on efficient squeezing of propagation eigenmodes
HHow to create entanglement among 100 spatial modes with a single photonic lattice
David Barral, ∗ Kamel Bencheikh, Juan Ariel Levenson, and Nadia Belabas Centre de Nanosciences et de Nanotechnologies C2N, CNRS,Universit´e Paris-Saclay, 10 boulevard Thomas Gobert, 91120 Palaiseau, France (Dated: May 14, 2020)Multimode entanglement is an essential resource for quantum networks. Continuous-variableencoding of quantum information in the optical domain has recently yielded large temporal andspectral entangled states instrumental for quantum computing and quantum communication. Weintroduce a protocol for the generation of spatial multipartite entanglement in a monolithic photonicdevice: the array of quadratic nonlinear waveguides. We theoretically demonstrate the generation oflarge multipartite entangled states in the spontaneous parametric downconversion regime. Our pro-tocol is remarkably simple and robust as it does not rely on specific values of coupling, nonlinearityor length of the sample.
Optical networks play a key role in our everyday lifeas the substrate of a long-range communication grid: theinternet. One goal of the blooming quantum technologiesis the development of a quantum internet: an informationweb with unparalleled capabilities with respect to its cur-rent classical counterpart where information will be pro-cessed by quantum computers, transmitted in quantumsecure channels, and routed towards quantum end nodes[1]. A must-have of the quantum internet is multipartiteentanglement, where information is strongly correlatedbetween the distributed nodes which compose the net-work [2]. In addition, multipartite entanglement is theresource of a number of protocols in quantum communi-cation, quantum sensing and quantum computing [3–5].Sources of multipartite entanglement are thus required inquantum networks and, particularly, light-based sourcesat telecom wavelengths are favored due to the currentavailability of large optical fiber networks.Quantum information can be encoded in variablesthat can take a continuous spectrum of eigenvalues –continuous variables (CV)– [6]. In the optical domain,the fluctuations of the field quadratures can be used ascarriers of quantum information [7]. A number of table-top experiments have demonstrated CV quantum net-works in the spatial, frequency and temporal domains[8–13]. The transfer to real technologies is howeverfar from feasible with bulk-optics systems. Scalability,stability and cost are issues that only well-establishedtechnologies like integrated and fiber optics can over-come [14]. Generation of two-mode CV entanglementthrough bulk-integrated hybrid approaches has been ex-plored [15, 16] and, remarkably, the first demonstration ofa fully-on-chip source of CV bipartite entanglement hasbeen recently proposed [17]. Nevertheless, the extensionof that bulk-optics-based scheme –sequential squeezingand entanglement– to larger number of modes is verydemanding with current technology. In this paper wepresent a simple and practical protocol for the generationon chip of spatial multipartite entanglement of sponta- ∗ [email protected] neous parametric downconverted (SPDC) light based ona currently-available technology: the array of χ (2) non-linear waveguides (ANW) [18–21].Harnessing the ANW for any purpose including multi-partite entanglement is challenging: the continuous inter-play between nonlinearity and evanescent coupling pre-vents in general analytical solutions [22, 23]. This in-terplay also results in an extreme sensitivity to phys-ical parameters such as the device length and linearand nonlinear coupling strengths. Here, we analyticallydemonstrate multipartite entanglement among the SPDCmodes generated in the ( N + 1) / N . Our methodis based on efficient build-up of the zero supermode ofthe array when pumping with a flat profile. This proto-col is remarkably simple and robust as it does not relyon specific values of coupling, nonlinearity or length ofthe sample. Below we introduce the ANW, calculate an-alytical solutions for a flat pump profile and use thesesolutions to demonstrate multipartite entanglement.The ANW consists of N of identical χ (2) waveguides.In each waveguide, an input harmonic field at frequency ω h is downconverted into a signal field at frequency ω s .Pump-signal waves phase matching is produced only inthe coupling region and is set to produce degenerateSPDC light. The generated signal fields are then cou-pled through evanescent tails in contrast to the pumpfields which are not coupled due to a higher confinementin the waveguides. The physical processes involved aredescribed by the following system of equations [24] d ˆ A j dz = iC ( f j − ˆ A j − + f j ˆ A j +1 ) + 2 iη j ˆ A † j , (1)where ˆ A = ˆ A N +1 = 0, f = f N = 0 and j = 1 , . . . , N is the individual mode index. ˆ A j ≡ ˆ A j ( z, ω s ) is a monochromatic slowly-varying amplitudeannihilation operator of signal (s) photons related tothe mode propagating in the j th waveguide where[ ˆ A j ( z, ω ) , ˆ A † j (cid:48) ( z (cid:48) , ω (cid:48) )] = δ ( z − z (cid:48) ) δ ( ω − ω (cid:48) ) δ j,j (cid:48) . η j = g α h,j is the effective nonlinear coupling constant correspond-ing to the j th waveguide, with g the nonlinear constantproportional to χ (2) and the spatial overlap of the signal a r X i v : . [ qu a n t - ph ] M a y (cid:1) = (cid:2) (cid:1) = (cid:2) (cid:1) = (cid:2) (cid:1) = (cid:2) (cid:1) = (cid:2) FIG. 1. Sketch of the supermodes related to an array of linear waveguides with a homogeneous coupling profile (cid:126)f = (cid:126) N = 5. The horizontal axis stands for the individual modes. The propagation constants corresponding to each supermode are λ = {√ C , C , , − C , −√ C } . k ≡ l = 3 is the zero supermode. k (cid:54) = l are the side supermodes. and harmonic fields in each waveguide, and α h,j a strongcoherent pump field. The nonlinear coupling constant η j can be tuned by means of a suitable set of pump phaseand amplitude at each waveguide. C j = C f j is the lin-ear coupling constant between nearest-neighbor modes j and j + 1, with C the coupling strength and f j the cou-pling profile. z is the coordinate along the direction ofpropagation.We are interested here in studying CV features of theSPDC light generated in the ANW. The natural vari-ables in this framework are the quadratures of the opti-cal field. The quadratures ˆ x j , ˆ y j , where ˆ x j = ( ˆ A j + ˆ A † j )and ˆ y j = i ( ˆ A † j − ˆ A j ) are, respectively, the orthogonal am-plitude and phase quadratures corresponding to a signaloptical mode A j . The system of equations (1) can berewritten in terms of the quadratures as d ˆ ξ/dz = ∆( z ) ˆ ξ ,where ∆( z ) is a 2 N × N matrix of coefficients andˆ ξ = (ˆ x , ˆ y , . . . , ˆ x N , ˆ y N ) T . In general, this equation –or Equation (1)– can be solved numerically for a specificset of parameters C j , η j and N , or even analytically if N is small. Numerical or low-dimension analytical so-lutions do not provide much physical insight when in-creasing N . Remarkably, we have identified a case whereanalytical solutions are available for any N . This is thecase of a flat pump profile, i.e. pumping all the waveg-uides with η j = | η | . In this case the eigenmodes of thesystem are propagation eigenmodes –supermodes– [25].These eigenmodes form a basis and are represented byan orthogonal matrix M with real elements M k,j . Theindividual and propagation supermode bases are relatedby ˆ ξ S,k = (cid:80) Nj =1 M k,j ˆ ξ j . Figure 1 shows the supermodesof an array of 5 waveguides with a homogeneous cou-pling profile (cid:126)f = ( f , . . . , f N ) = (cid:126)
1. The supermodes areorthonormal (cid:80) Nj =1 M k,j M m,j = δ k,m , and the spectrumof eigenvalues is λ k . The solution of the propagation inthis basis can be written as ˆ ξ S,k ( z ) = S ( z ) ˆ ξ S,k (0), where S ( z ) = diag { S ( z ) , . . . , S N ( z ) } and S k ( z ) = (cid:18) cos( F k z ) − e − r k sin( F k z ) e r k sin( F k z ) cos( F k z ) (cid:19) , (2)with r k = (1 /
2) ln [( λ k + 2 | η | ) / ( λ k − | η | )] and F k = (cid:112) λ k − | η | . For typical coupling strengths and pumppowers found in quadratic ANW the condition | λ k | > | η | is fulfilled. This regime is the relevant one for entan-glement since, as the nonlinear interaction surpasses the linear coupling, the SPDC light tends to be more andmore confined in the waveguide where is created and thenthe ANW acts only as a group of individual squeezers[26]. We consider F k ∈ R in the remainder of the article.Notably, the analytical solution of Equation (2) is gen-eral for any evanescent coupling profile (cid:126)f , any number ofwaveguides N and any propagation distance z .The quantum states generated in ANW are Gaussian.The most interesting observables are then the second-order moments of the quadrature operators, properly ar-ranged in the covariance matrix V ( z ) [27]. For a quan-tum state initially in vacuum, the elements of the covari-ance matrix V ( z ) can be straightforwardly obtained fromEquation (2) using V ( z ) = S ( z ) S T ( z ). The covariancematrix elements in the supermode basis are then V ( x S,k , x
S,k ) = [cosh ( r k ) + sinh ( r k ) cos (2 F k z )] e − r k ,V ( y S,k , y
S,k ) = [cosh ( r k ) − sinh ( r k ) cos (2 F k z )] e + r k ,V ( x S,k , y
S,k ) = sinh ( r k ) sin (2 F k z ) . (3)The squeezing ellipses for each supermode vary alongpropagation. Maximum (respectively null) squeezing areobtained periodically at distances L k = (2 n − π/ (2 F k )(respectively L (cid:48) k = nπ/F k ) that are different for each k th supermode, with n any positive integer. The maxi-mum value of squeezing available in the k th supermodeis e − r k = ( λ k − | η | ) / ( λ k + 2 | η | ).Instrumentally for our protocol, waveguide arrays withan odd number of identical waveguides present a super-mode with zero eigenvalue –the zero supermode– [28].This corresponds to the central supermode with k ≡ l =( N +1) / λ l = 0 as shown in Figure 1. The elementsof the covariance matrix for the zero supermode are thus V ( x S,l , x
S,l ) = V ( y S,l , y
S,l ) = cosh (4 | η | z ) ,V ( x S,l , y
S,l ) = sinh (4 | η | z ) . (4)In contrast to the side supermodes k (cid:54) = l , the zero super-mode noise efficiently builds up at all propagation dis-tances and, notably, for large coupling strength the zerosupermode is quickly dominant over the side supermodes.The reason is that this is the only supermode which isphase-matched along propagation ( λ l = 0).The covariance matrix Equation (3) in the supermodebasis shows the total squeezing available in the downcon-verted signal fields. However, that basis is diagonal andno entanglement is available in that encoding. A sim-ple change of basis takes Equation (3) to the individualmode basis, corresponding to the individual waveguidesoutput, obtaining V ( ξ i , ξ j ) = N (cid:88) k =1 M i,k M j,k V ( ξ S,k , ξ
S,k ) , (5)with ξ ≡ x or y . Hence, the flat pump configuration gen-erates quantum correlations –off-diagonal components ofthe covariance matrix– between any pair i and j of in-dividual modes, and thus full inseparability among indi-vidual modes can be produced.Measuring multipartite full inseparability in CV sys-tems requires the simultaneous fulfillment of a set of con-ditions which leads to genuine multipartite entanglementwhen pure states are involved [29]. This criterion, knownas van Loock - Furusawa (VLF) inequalities, can be eas-ily calculated from the elements of the covariance matrix V . Full N -partite inseparability is guaranteed if the fol-lowing N − VLF j ≡ V [ x j ( θ j ) − x j +1 ( θ j +1 )]+ V [ x j ( θ j + π x j +1 ( θ j +1 + π N (cid:88) m (cid:54) = j,j +1 G m x m ( θ m + π ≥ , (6)where ˆ x j ( θ j ) = ˆ x j cos ( θ j )+ˆ y j sin ( θ j ) are generalized quadra-tures, θ j is the measurement phase corresponding to the j thlocal oscillator (LO), and G , . . . , G N are N real parame-ters corresponding to electronic gains in multimode balancedhomodyne detection (BHD) which are set by optimization. V [ (cid:80) j σ j ξ j ] ≡ (cid:80) j σ j V ( ξ j , ξ j ) + (cid:80) i (cid:54) = j σ i σ j V ( ξ i , ξ j ) with σ j aset of real numbers. (cid:126)θ ≡ ( θ , . . . , θ N ) and (cid:126)G ≡ ( G , . . . , G N )stand, respectively, for the LO phase and gain profiles. Thesetwo families of parameters can be set in order to minimizesuitably the value of Equation (6). Below, we show a remark-ably simple way of generating multipartite entanglement inANW with an odd number of waveguides based on the effi-cient squeezing of the SPDC zero supermode.Firstly, in order to gain insight about the multipartite en-tanglement generated with the flat pump configuration wetackle the limit of large coupling ( C → ∞ ). This limit isnot physical as next-nearest-neighbor coupling should be inthat case included in the model, but it gives us a clear insighton the dynamics of the system as the zero-supermode is thenthe dominant supermode generated in the array. In this limitan asymptotic lower bound on the violations of the VLF in-equalities is obtained for the non-optimized case (cid:126)G = (cid:126)
0. Thecovariance matrix elements Equations (5) for an array withodd number N of waveguides in the limit of large coupling( C → ∞ ) is significantly simplified to V ( x i , x j ) = V ( y i , y j ) → δ i,j + 2 M i,l M j,l sinh (2 | η | z ) ,V ( x i , y j ) → M i,l M j,l sinh(4 | η | z ) . (7)Applying this result into the general expression for the VLFinequalities Equation (6) without optimization ( G m (cid:54) = j,j +1 =0) and using generalized quadratures with θ j = 0 and θ j +1 = − π/ π/
2, we obtainVLF j ( ∞ ,(cid:126) , z ) = 4 − M j,l + M j +1 ,l )+ ( M j,l ± M j +1 ,l ) e | η | z + ( M j,l ∓ M j +1 ,l ) e − | η | z ≥ , (8) (cid:1) = (cid:2)(cid:1) = (cid:3)(cid:1) = (cid:4)(cid:1) = (cid:5)(cid:6) (cid:7)(cid:6) (cid:3)(cid:6) (cid:2)(cid:6) (cid:5)(cid:6) (cid:4)(cid:6)(cid:3)(cid:8)(cid:6)(cid:3)(cid:8)(cid:4)(cid:2)(cid:8)(cid:6)(cid:2)(cid:8)(cid:4)(cid:5)(cid:8)(cid:6) (cid:9) ( (cid:10)(cid:10) ) (cid:1) (cid:2)(cid:3) FIG. 2. Optimized van Loock - Furusawa inequalitiesVLF( C , (cid:126)G, z ) for a small number of involved modes (color).Simultaneous values under the threshold value VLF=4 im-ply CV tripartite entanglement (two inequalities, l=3, dot-ted), quadripartite entanglement (three inequalities –twodegenerate–, l=4, dot-dashed) and pentapartite entanglement(four inequalities –degenerate two by two–, l=5, solid). Wealso show bipartite entanglement for comparison (one in-equality, l=2, dashed). The optimized inequalities in thelimit of large coupling VLF( ∞ , (cid:126)G, z ) are shown in solid gray. C = 0 .
70 mm − . η = 0 .
025 mm − . These are typical valuesfor periodically poled lithium niobate waveguides [17, 19].where the upper (lower) signs corresponds to θ j +1 = − π/ π/
2) and we have introduced the notation VLF ≡ VLF( C , (cid:126)G, z ) for the sake of clarity. The best scenario interms of violation of these inequalities corresponds to the case M j,l = ∓ M j +1 ,l , for which we obtainVLF j ( ∞ ,(cid:126) , z ) = 4 − M j,l (1 − e − | η | z ) < ∀ z > . (9)Particularly, the coefficients of the zero supermode in anarray with homogeneous coupling profile (cid:126)f = (cid:126) M j,l = sin( jπ ) / √ l (Figure 1) [31]. Hence, mapping themode 2 j − j of Equation (8), two solu-tions which maximize the violation of the separability con-ditions are obtained: i) the l odd elements of the zero super-mode satisfy M j − ,l = − M j +1 ,l such that for a LO profile( θ j − , θ j +1 ) = (0 , − π/
2) multipartite entanglement is ob-tained among all the odd individual modes { , , , . . . N } ,and ii) the odd elements of the zero supermode satisfy M j − ,l = M j +1 ,l , thus for ( θ j − , θ j +1 ) = (0 , π/
2) themultimode state is decoupled in two multipartite entangledstates: { , , , . . . } and { , , , . . . } . Thus, the LO phaseprofile acts as an entanglement switch between two multi-mode entangled states. The individual modes propagating inthe odd waveguides are fully inseparable in a measurementbasis and separable in two parties –each with elements fullyinseparable– in the other. The following degenerate violationof the inseparability conditions is obtained in both casesVLF j ( ∞ ,(cid:126) , z ) = 4 ( l −
1) + e − | η | z l < ∀ z, l. (10)Hence, the strength of the violation of the VLF inequalitiesdepends asymptotically only on the number of odd individualmodes l which make up the zero supermode. Moreover, theuse of an optimized gain profile (cid:126)G (cid:54) = (cid:126) (cid:1) = (cid:2)(cid:3)(cid:3)(cid:1) = (cid:4)(cid:3)(cid:1) = (cid:5)(cid:4)(cid:3) (cid:2)(cid:3) (cid:5)(cid:3) (cid:6)(cid:3) (cid:7)(cid:3) (cid:4)(cid:3)(cid:6)(cid:8)(cid:9)(cid:3)(cid:6)(cid:8)(cid:9)(cid:4)(cid:6)(cid:8)(cid:10)(cid:3)(cid:6)(cid:8)(cid:10)(cid:4)(cid:7)(cid:8)(cid:3)(cid:3) (cid:11) ( (cid:12)(cid:12) ) (cid:1) (cid:2)(cid:3) FIG. 3. Optimized van Loock - Furusawa inequalities forlarge number of involved modes in the limit of large cou-pling VLF( ∞ , (cid:126)G, z ). Simultaneous values under the thresholdvalue VLF=4 imply CV multipartite entanglement. η = 0 . − .above result. We have found the following optimized violationof the VLF conditions which depends also on the parity of l (see supplemental material)VLF j ( ∞ , (cid:126)G, z ) = VLF j ( ∞ ,(cid:126) , z )+ − l − l +3)(1 − e − | η | z ) l [2 l ( l − − ( l − l +3)(1 − e − | η | z )] < ∀ z, l odd , − l − − e − | η | z ) l [2 l − ( l − − e − | η | z )] < ∀ z, l even . (11)The correction produced optimizing (cid:126)G is zero for l = 2 , l >
3. It scales as l − in the limit of a largenumber of modes. Therefore, we have demonstrated that ourprotocol always produces multipartite entanglement in ANWin the limit of large coupling.When finite coupling C is taken into account, the sidesupermodes are present in the optimization of the VLF in-equalities. This generates fluctuations around the valueVLF j ( ∞ , (cid:126)G, z ). Figure 2 (color) shows one, two, three andfour inequalities for arrays with, respectively, l =2, 3, 4 and 5propagating modes obtained through minimization of Equa-tions (6) with a suitable gain profile (cid:126)G where we have used theanalytical solutions of Equations (5) [30]. The simultaneousviolation of the inequalities (VLF j < C → C → ∞ and optimized gain profiles (cid:126)G [Equation (11)] for each case.Figure 2 thus demonstrates the validity of Equation (11) as amean value of the real VLF inequalities that can be used toassess the possible entanglement generated in the array.Remarkably, Equation (11) shows that quantum correla-tions are exhibited at any z independently of the number l ofmodes involved. Figures 3 and 4 depict respectively the evolu-tion of multipartite entanglement along propagation for largenumber of modes ( l = 25, 50 and 100), and the relationshipof the asymptotic value of VLF inequalities with the numberof involved modes l for z → ∞ . These figures demonstrate (cid:1) (cid:2)(cid:1) (cid:3)(cid:1) (cid:4)(cid:1) (cid:5)(cid:1) (cid:6)(cid:1)(cid:1)(cid:2)(cid:7)(cid:1)(cid:2)(cid:7)(cid:8)(cid:9)(cid:7)(cid:1)(cid:9)(cid:7)(cid:8)(cid:3)(cid:7)(cid:1) (cid:10)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15) (cid:16)(cid:17) (cid:14)(cid:18)(cid:19)(cid:20)(cid:18)(cid:21)(cid:22)(cid:14)(cid:23) (cid:12)(cid:16)(cid:23)(cid:14)(cid:24) (cid:1) (cid:2)(cid:3) FIG. 4. Optimized van Loock - Furusawa inequalities in thelimit of large coupling and propagation length VLF( ∞ , (cid:126)G, ∞ ).Simultaneous values under the threshold value VLF=4 implyCV multipartite entanglement. η = 0 .
025 mm − the scalability of our protocol. Thus, using our source of mul-tipartite entanglement, fidelities of quantum teleportation inmultimode networks better than what could be achieved inany classical scheme are expected [2]. Noticeably, the asymp-totic violation of the VLF inequalities in the double asymp-totic limit ( C , z → ∞ ) given by VLF j ( ∞ ,(cid:126) , ∞ ) = 4( l − /l is the same as that obtained in second harmonic generation(SHG) when a zero supermode is excited at the input of theANW [31]. Unlike the SHG case where only the zero super-mode is present, here the k (cid:54) = l side supermodes are involvedin the production of entanglement. They increase the viola-tion of the inequalities along z through the use of optimizedgains (cid:126)G as shown in Equation (11). The asymptotic behav-ior exhibited in Figures 2 and 3 appears as a consequence oftracing over the fields present in the even channels [31].Finally, we would like to underline that this configurationis very appealing for the generation of scalable multipartiteentanglement since it relies on coupling C and nonlinearity g within the array, but not on specific values of these parame-ters. Our protocol yields significant and useful entanglementover a wide range of number of modes. Furthermore, ourmethod gives insight on further extension of the possibilitiesof the ANW for a resource-efficient generation of large entan-gled states. For example, trying to phasematch every super-mode of the array. This would represent a compact, efficientand fully scalable way to produce multipartite entanglement.We are currently studying possible solutions in that directionusing for example supermode quasi-phasematching [32]. An-other approach which can improve the entanglement and evengenerate large cluster states for measurement-based quantumcomputing is the optimization of pump and LO profiles. Suit-able optimization of selected tuning parameters indeed turnthe ANW in a versatile spatially-encoded entanglement syn-thesizer [23].This work was supported by the Agence Nationale de laRecherche through the INQCA project (Grants No. PN-II-ID-JRP-RO-FR-2014-0013 and No. ANR-14-CE26-0038),the Paris Ile-de-France region in the framework of DIMSIRTEQ through the project ENCORE, and the Investisse-ments d’Avenir program (Labex NanoSaclay, reference ANR-10-LABX-0035). [1] S. Wehner, D. Elkouss and R. Hanson. Science , 303(2018).[2] P. van Loock and S.L. Braunstein.
Phys. Rev. Lett. ,3482 (2000).[3] M. Hillery, V. Buzek and A. Berthiaume. Phys. Rev. A , 1829-1834 (1999).[4] Q. Zhuang, Z. Zhang and J. .H. Shapiro. Phys. Rev. A , 032329 (2018).[5] R. Raussendorf and H.J. Briegel. Phys. Rev. Lett. ,5188 (2001).[6] S.L. Braunstein and P. van Loock. Rev. Mod. Phys. ,513 (2005).[7] U.L. Andersen, G. Leuchs and Ch. Silberhorn. Laser &Phot. Rev.
4, 337-354 (2010).[8] X. Su, S. Hao, X. Deng, L. Ma, M. Wang, X. Jia, C. Xie,and K. Peng.
Nature Comm. , 2828 (2013).[9] S. Armstrong, M. Wang, R.Y. Teh, Q. Gong, Q. He,J. Janousek, H.-A. Bachor, M.D. Reid and P.K. Lam. Nature Phys. , 167-172 (2015).[10] M. Chen, N.C. Menicucci and O. Pfister. Phys. Rev. Lett. , 120505 (2014).[11] Y. Cai, J. Roslund, G. Ferrini, F. Arzani, X. Xu, C. Fabreand N. Treps.
Nature Comm. , 15645 (2017).[12] J.-I. Yoshikawa, S. Yokoyama, T. Kaji, Ch. Sornphiphat-phong, Y. Shiozawa, K. Makino and A. Furusawa. APLPhotonics , 060801 (2016).[13] S. Takeda, K. Takase and A. Furusawa. Science Advances (5), eaaw4530 (2019).[14] J. Wang, F. Sciarrino, A. Laing and M.G. Thompson. Nature Phot. , s41566-019-0532-1 (2019).[15] G. Masada, K. Miyata, A. Politi, T. Hashimoto, J.L.O’Brien and A. Furusawa.
Nature Phot. , 316-319(2015).[16] M.V. Larsen, X. Guo, C.R. Breum, J.S. Neergaard-Nielsen and U.L. Andersen. NPJ Quant. Inf. , 46 (2019).[17] F. Lenzini, J. Janousek, O. Thearle, M. Villa, B. Haylock,S. Kasture, L. Cui, H.-P. Phan, D.V. Dao, H. Yonezawa,P.K. Lam, E.H. Huntington and M. Lobino. Science Ad-vances (12), eaat9331 (2018).[18] D.N. Christodoulides, F. Lederer and Y. Silberberg. Na-ture (6950), 817-23 (2003).[19] R. Iwanow, R. Schiek, G.I. Stegeman, T. Pertsch, F. Led-erer, Y. Min and W. Sohler.
Phys. Rev. Lett. (11),113902 (2004).[20] F. Setzpfandt, D.N. Neshev, R. Schieck, F. Lederer, A.Tunnermann and T. Pertsch. Opt. Lett. (22), 3589-3591 (2009).[21] A.S. Solntsev, F. Setzpfandt, A.S. Clark, C.W. Wu, M.J.Collins, C. Xiong, A. Schreiber, F. Katzschmann, F.Eilenberger, R. Schiek, W. Sohler, A. Mitchell, Ch. Sil-berhorn, B.J. Eggleton, T. Pertsch, A.A. Sukhorukov,D.N. Neshev and Y.S. Kivshar. Phys. Rev. X , 031007(2014).[22] D. Barral, N. Belabas, L.M. Procopio, V. D’Auria, S.Tanzilli, K. Bencheikh and J.A. Levenson. Phys. Rev. A , 053822 (2017).[23] D. Barral, M. Walschaers, K. Bencheikh, V. Parigi, J.A. Levenson, N. Treps and N. Belabas. arXiv:1912.11154v3 ,(2020).[24] J. Li˜nares, M.C. Nistal and D. Barral. New J. Phys. ,063023 (2008).[25] E. Kapon, J. Katz and A. Yariv. Opt. Lett. (4), 125-127 (1984).[26] J. Fiurasek, J. Perina. Phys. Rev. A , 033808 (2000).[27] G. Adesso, S. Ragy and A.R. Lee. Open Syst. Inf. Dyn. , 1440001 (2014).[28] N.K. Efremidis and D.N. Christodoulides. Opt. Comm. , 345 - 356 (2005).[29] P. van Loock and A. Furusawa. Phys. Rev. A , 052315(2003).[30] S.L.W. Midgley, A.S. Bradley, O. Pfister and M.K. Olsen. Phys. Rev. A , 063834 (2010).[31] D. Barral, K. Bencheikh, N. Belabas and J.A. Levenson. Phys. Rev. A , 051801(R) (2019).[32] D. Barral, N. Belabas, K. Bencheikh and J.A. Levenson. Phys. Rev. A100