aa r X i v : . [ phy s i c s . op ti c s ] J un Quantum gravity simulation by non-paraxial nonlinear optics
Claudio Conti Institute for Complex Systems ISC-CNR and Department of Physics,University Sapienza, Piazzale Aldo Moro 5, 00185, Rome (IT) ∗ (Dated: September 17, 2018)We show that an analog of the physics at the Planck scale can be found in the propagation oftightly focused laser beams. Various equations that occur in generalized quantum mechanics areformally identical to those describing the nonlinear nonlocal propagation of nonparaxial laser beams.The analysis includes a generalized uncertainty principle and shows that the nonlinear focusing of alight beam with dimensions comparable to the wavelength corresponds to the spontaneous excitationof the so-called maximally localized states. The approach, driven by the ideas of the quantum gravityphysics, allows one to predict the existence of self-trapped subwavelength solitary waves for bothfocusing and defocusing nonlinearities, and opens the way to laboratory simulations of phenomenathat have been considered to be inaccessible. A rapidly expanding research direction in theoreticalphysics concerns the investigation of a generalized un-certainty principle (GUP).[1–6] Quantum gravity (QG)models predict fuzziness and discretization in the geom-etry, which result in the fact that the uncertainty of thespatial coordinate ∆ X cannot be smaller than a minimalquantity ∆ X min . The simplest GUP [1] reads as∆ X ∆ P ≥ ¯ h β ∆ P ), (1)with ∆ P the momentum uncertainty. Equation (1) im-plies ∆ X ≥ ¯ h √ β ; standard quantum mechanics (QM)is retrieved for β = 0. Various upper bounds to β havebeen proposed. [5, 7] Among the many consequences ofthe GUP, we mention corrections to the black-body radi-ation spectrum, and to the cosmological constant. [8, 9]Analyses of GUP effects in the mechanical vibration ofmacroscopic objects have been reported. [10–13] Theo-retical developments include the existence of a maximalmomentum P . [14] This may be related to another knownfundamental limit to ∆ X , from special relativity (SR),given by the reduced Compton wavelength λ C = ¯ h/mc [15], with m the particle rest mass and c the vacuumlight velocity. λ C is much larger than the Planck scale l P ∼ = 10 − m, which is typically retained as ∆ X min inQG ( λ C ∼ = 10 − m for the electron).In several respects, the analysis of nonlinear wave equa-tions with GUP is un-explored. In QG nonlinear ef-fects are expected to be extremely relevant because ofthe energy levels needed to access the Planck scale, andbecause of the nonlinearity due to background indepen-dence. Nonlinear mechanisms are also known to be due tothe particle-anti-particle production expected at the scaleof the Compton wavelength [15], and nonlinear modifica-tions to Maxwell equations have been predicted in loopquantum gravity [16, 17]. In condensed matter physics,nonlinearity is due to the atom-atom interaction in ultra-cold gases and Bose-Einstein condensates, where investi-gations of the effects of Planck scale physics have beenalso reported [18–20]. Here we show that all the mentioned theoretical devel-opments naturally apply to the description of tightly fo-cused laser beams[21] with size smaller than the wave-length and propagating in local and nonlocal nonlin-ear media.[22, 23] This analogy allows to obtain an ex-plicit expression for the β parameter. When includingnonlinearity, we find that a local, or nonlocal, intensitydependent refractive index perturbation forces the self-trapped beam to acquire a shape correspondent to theso-called maximally localized states in QG. Previouslyinvestigated QG models, as the generalized harmonic os-cillator, directly describe nonlinear waves in the presenceof an highly nonlocal nonlinearity. In addition, the natu-ral discreteness of quantized geometry enables to map acontinuous equation, as the one describing non-paraxialoptical solitons, to an exact discrete model. The tech-niques derived from the GUP framework allow to predictthe existence of sub-wavelength optical beams for bothfocusing and defocusing nonlinearities, and open the wayto the direct experimental tests of the generalized quan-tum mechanics supposed to be valid at the Planck scale. Projective GUP —
The simplest generalized Schr¨odingerequation (GSE) sustaining a GUP is [5] ı ¯ h∂ t ψ = ˆ p m ψ + βm ˆ p ψ , (2)with ˆ p = − ı ¯ h ∇ the standard momentum operator. Thisequation arises when projecting the Helmholtz equationfor the electromagnetic field E with wavenumber k inthe forward direction z , which is the equation that de-scribes light propagation beyond the paraxial approxima-tion (see, e.g, [24] and references therein): i∂ z E + p ∇ + k E = 0. (3)By ψ = E exp( ikz ), Eq.(3) is written as ıλ∂ z ψ = ˆ P z ψ = h − p − ( − ıλ ∇ ) i ψ (4)with λ = 2 π/k the wavelength, and λ = λ/ (2 π ). ∇ isthe gradient with respect to the transverse coordinates( x, y ). By mc ≡ ¯ hω = ¯ hc/λ , with ω = 2 πc/λ , being T the laboratory time, with t = z/c − T , we have ı ¯ h∂ t ψ = ˆ Hψ = mc ˆ P z ψ , (5)which, for small λ , gives Eq.(2) with β = 38 1 m c = 38 (cid:18) λh (cid:19) . (6)It is remarkable that this approach allows to obtain anexpression for the β parameter, which may be general-ized to other fields. Letting G the gravitational constant,and M P = (¯ hc/G ) / the Planck mass, the normalizedcoefficient β = ( M P c ) β is considered, and we have β = 38 M P m = 38 ¯ hcG c λ ¯ h = 38 c λ G ¯ h . (7)When applied to the photon with wavelength λ = 1 µ m,Eq.(7) gives β = 10 . In [5] it has been estimated β < ; we hence observe that, in this analogue of QG,the GUP effects are several order of magnitudes greaterthan those previously estimated. This shows that GUPeffects can be directly observed in the laboratory and alsofurnishes a direct expression for the parameter β . Eq.(7)can be generalized to an higher number dimensions. Maximally localized states —
Limiting our analysis toone dimension x , we have ˆ H = ˆ P / m , beingˆ P = √ mc vuut − s − (cid:18) ˆ pmc (cid:19) ˆ p | ˆ p | , (8)the generalized momentum. When β → P =ˆ p (cid:16) β ˆ p (cid:17) . In addition, ˆ p = ˆ P q − ( ˆ P / (2 mc )) =ˆ P (1 − β ˆ P /
3) + O ( β ). In the momentum representation,waves have finite support in the interval p ∈ [ − mc, mc ],corresponding to P ∈ [ −√ mc, √ mc ]. This approachleads to states with a finite support for P , as those intro-duced by Pedram in [14]. We have (the prime denotingthe derivative)[ ˆ X, ˆ P ] = ı ¯ h q − ( ˆ P / mc ) − ˆ P / m c = ı ¯ hp ′ ( ˆ P ) (9)that, for β →
0, reduces to the KMM model [1][ ˆ X, ˆ P ] = ı ¯ h (1 + β ˆ P ). (10)The GUP (1) arises at the lowest order in β from∆ X ∆ P ≥ |h [ ˆ X, ˆ P ] i| with h ˆ P i = 0. This is equivalent tothe original KMM proposal [1], with ˆ X = ˆ x , ˆ P = ˆ p + β ˆ p being [ˆ x, ˆ p ] = ı ¯ h .Being a minimal length uncertainty, the eigenstates ofthe position operator ˆ X are not physically realizable, as they correspond to ∆ X = 0. In the P − representation ψ ( P ) = h P | ψ i , and we have for the position operatorˆ X → ı ¯ hp ′ ( P ) ∂ P ψ ( P ). (11)Letting ˆ P z = ( cP ) /
2, one has ı ¯ hψ t ( P ) = P m ψ ( P ) (12)so that the P − eigenstates correspond to free particleswith kinetic energy P / m , which reduces to the stan-dard p / m in the small momentum limit β → P − representation is very useful to determine themaximally localized (ML) states satisfying the conditionof minimal ∆ X . Following [25], ML states with h X i = h P i = 0 are given by ((cid:20) ı ¯ h∂ P p ′ ( P ) (cid:21) − µ ) Φ ML ( P ) = 0, (13)with µ = ∆ X , which leads to ( n = 1 , , ... )Φ MLn ( P ) = r ¯ hmc sin h nπ mc p ( P ) + nπ i . (14)The maximal localization corresponds to n = 1, with µ = ∆ X min = ¯ hπ/ mc = λ/ Nonlinearity —
We consider the following nonlinear waveequation ı ¯ h∂ t ψ = ˆ P m ψ + Z κ ( x − x ′ ) | ψ ( x ′ ) | dx ′ ψ , (15)where κ ( x ) is a kernel function. We consider a nonlo-cal nonlinearity (NN) [26] because the interaction lengthmay be much larger than the geometry quantizationlength. NN supports stable regimes, [27] with exact so-lutions in the highly [28] and weakly nonlocal limits [29].Eq.(15) may be generalized to include vectorial effects.The strength of the nonlinearity is determined by thenorm N = R | ψ | dx .We first analyze the highly nonlocal (HN) case with κ ( x ) ∗ | ψ ( x ) | ∼ = N κ ( x ) and the asterisk denoting theconvolution integral. [30] Letting N κ ( x ) = κ N x / ≡ ¯ hc Ω x /
2, with κ and Ω = κ N coefficients of the Tay-lor expansion of κ ( x ), we have ˆ H = ˆ P / m + m Ω x / ψ ( P, t ) = φ ( P ) exp( − ı E ¯ h t ) with P m φ ( P ) − m Ω ¯ h p ′ ( P ) ∂∂P (cid:20) p ′ ( P ) ∂φ∂P (cid:21) = Eφ ( P ). (16)In the small momentum limit Eq. (16) reduces to ∂ φ∂P + 2 βP βP ∂φ∂P + (cid:18) Em ¯ h − P m Ω ¯ h (cid:19) φ (1 + βP ) = 0,(17) ϵ =p e /2 ϵ =e nonlinearity e e n e r g y ϵ harmonicoscillatorground statemaximallylocalizedstate momentum y harmonicoscillatorground state maximallylocalizedstate g r o un d s t a t e F FIG. 1. (Color online) (a) Profiles of the ground state Φ( y )versus the dimensionless generalized momentum y for threevalues of the nonlinearity ε ; (b) energy ǫ versus ε . which is the QG harmonic oscillator in [1]. We writeEq.(16) in dimensionless units, by y = P/ √ mc and ǫ = E/mc − g ( y ) ddy (cid:20) g ( y ) ddy Φ( y ) (cid:21) + y ε Φ( y ) = ǫε Φ( y ) (18)with Φ( ±
1) = 0, g ( y ) = p − y / / (1 − y ), and ε =¯ h Ω / mc .According to the Sturm-Liouville theory [31], as g ( y ) >
0, Eq.(18) admits discrete solutions Φ n ( P ) and a set ofpositive values ǫ n of ǫ , with n = 0 , , , ... the numberof zeros of Φ n ( P ). For ε →
0, in the small nonlinearitylimit, a multiple scale expansion with Y = y/ √ ε andΦ( y ) = Φ( y/ √ ǫ ) gives the standard quantum harmonicoscillator (QHO), with E = ¯ h Ω( n + 1 / n = 0,Φ( y ) = exp ( − εy /
2) and ǫ = ε . For large nonlinearity ε → ∞ , Eq.(18) reduces to Eq. (13) and the states tendto the ML: Φ( y ) = √ πy p − y / / √
2) with ǫ = ε π /
2, i.e., E = π (¯ h Ω) / mc .To verify these limits, Eq.(18) is numerically solved bywriting Φ( y ) as a superposition of Chebyshev polynomi-als [32]. In Fig. 1a we show the ground state solutionΦ( y ) for three strengths of the nonlinearity ε , and in Fig.1b we show the eigenvalue ǫ = ǫ versus ε . The solutionsof Eq.(18) interpolate between the QHO ground stateand the ML state.In the quasi-position space, the solution can be rep-resented by cardinal functions φ C = √ λ sin( x/λ ) / √ πx ,letting x n = πλ (2 n + 1) / ψ ( x, t ) = ∞ X n = −∞ √ πλψ ( x n , t ) φ C ( x − x n ), (19)with √ πλφ C ( x n − x q ) = δ qn the Kronecker symbol. Weuse this expansion in the local case in Eq. (15), with κ ( x ) = χm c δ ( x ),with δ ( x ) the Dirac δ and χ = ± ı ¯ h∂ t ψ = ˆ Hψ − χmc | ψ | ψ . (20)By using Eq. (19), Eq.(20) is mapped to a discrete model.By the Einstein convention, we have ı ¯ hmc dψ n dt = h qn ψ q − χ | ψ n | ψ n , (21)being ψ n ( t ) = ψ ( x n , t ) and h qn = " ˆ P ψ C ( x − x p )2 m c x = x q = δ qp − J [ π ( n − p )]2( n − p ) . (22)The matrix h qn is the cardinal basis representation of theoperator ˆ P / (2 m c ) and is positive-definite. We alsohave, letting ψ n ≡ ψ n , N = πλ ∞ X n = −∞ | ψ n | = πλψ n ( ψ n ) ∗ . (23)Bound states of (21) are ψ n = φ n exp( − ıEt/ ¯ h ) with h qn φ q − χφ n = Emc φ n . (24)To study the ground state solution of Eq.(24), we adopta perturbation expansion in 1 / N , with the ML as leadingorder, φ MLn = r N λ (cid:0) δ n + δ − n (cid:1) . (25)For N → ∞ the solution is written as φ n = φ MLn + 1 N φ (1) n + o ( N − ) (26)and E/mc ≡ ǫ = ǫ ML + N ǫ (1) . By Eq.(24), φ (1) n = (1 − δ n − δ -1 n ) χ √N λ (cid:26) J ( πn ) n + J [( n + 1) π ] n + 1 (cid:27) ,(27) ǫ (1) = − χ/λ , and ǫ = Emc = h + h -10 − χ N λ = 1 − π − J ( π )2 − χ N λ , (28)being h = 1 − π/ >
0, and h -10 = − J ( π ) / <
0. TheSW tends to ML for χ = ±
1, i.e., h φ ML | φ i = 1 + o ( N -1 ).To validate this analysis, we numerically solve Eq.(24) bya Newton-Raphson algorithm. We first consider the case χ = 1 and we show in Fig. 2a the ground state expressedin terms of the discrete values φ n and after Eq.(19).When increasing N the wavefunction tends to a MLstate, with only two non-vanishing coefficients φ and φ -1 corresponding to x = πλ = λ/ x -1 = − πλ = − λ/ ǫ which, for large N ,tends to Eq.(28). We also show in Fig.2b the calculated∆ X/πλ = [ P q (2 q + 1) φ q / P φ q ] / , which tends to theminimal uncertainty ∆ X min /πλ = 1 / χ = −
1, asshown in Fig.2c,d; following Eq.(28), the energy ǫ is pos-itive and the algorithm converges to a localized SW for N > N C ∼ = 0 . λ .We investigate the local SW stability by linearizingEq.(21) by ψ n = ( φ n + η n ) exp ( − ıEt/ ¯ h ), (29)which gives ı ¯ hmc dη n dt = h qn η q − ǫη q + χφ n ( η n + 2 η ∗ n ). (30)We numerically solve Eq.(30) for exponentially divergingsolutions. φ n for χ = − χ = 1we find a real-valued unstable eigenvalue α growing with N . This is analytically verified by η n = η δ n + η -1 δ -1 n and η ± = η ± η -1 , which gives for N → ∞ and φ n ∼ = φ MLn , ı dη ± dt = ( h ± h -10 ) η ± − ǫη ± + χ N λ ( η ± + 2 η ∗± ), (31)Letting η ± = ˆ η ± exp( α ± mc t/ ¯ h ), we have η + = 0, and α − = 4 | h -10 | (cid:20) h -10 + χ N λ (cid:21) . (32)For χ = 1, Eq.(32) gives an unstable mode α − > N , | α − | = [2 J ( π ) N /λ ] / in agreement with the numerical solutions of Eq. (30).This mode exists for N /λ > N th /λ = J ( π ) / ∼ = 0 . χ = − α − is negative, and the mode is stableand corresponds to an energy shift | α − | . As the energy ǫ ± | α − | must be positive this gives the existence bound-ary N /λ > N C /λ ∼ = 2 J ( π ) ∼ = 0 . Discussion—
The study here reported is based on a one-dimensional (1D) model and correspondingly neglectsvectorial effects, which do not occur for 1D linearly po-larized beams. For focusing and defocusing nonlinear-ities, in the limits of highly nonlocal and local nonlin-earity, non-paraxial self-trapped beams can be predictedby the analogy with QG. During evolution, nonlinearitymay enhance, or reduce, the level of localization, and theparticular regime in terms of the amount of nonlocality.Further analysis is needed to determine if the ML statesmay be spontaneously generated by paraxial spatially ex-tended beams. On the contrary, the predicted QG-drivenSW are expected to be accessible when the initial spatialextension is comparable with the wavelength. It may alsohappen that the wavepacket periodically evolves back-and-forth different nonlocal regimes. This correspondsto periodical QG-like effects when the size of the beamreaches the equivalent of the Planck scale.The extension to two-dimensional (2D) wavepacketsrequires additional mathematical concepts, as reported N / (cid:1) x(units of (cid:0) /2) 0 1 2 3-5 0 5 N C / (cid:2) , X ( un i t s o f / ) XX (a) (b)(c) (d) (cid:4) =1 (cid:3) =-1 X min FIG. 2. (Color online) (a) SW profiles scaled to unitary normfor χ = 1, calculated in the cardinal basis after Eq.(19) andthe numerical solution of Eq. (24), for N /λ = 0 .
075 (circles), N /λ = 0 .
134 (triangles), N /λ = 3 .
00 (squares); (b) ǫ and∆ X when increasing the strength of nonlinearity ; (c) as in(a) for χ = − N /λ = 0 .
85 (circles), N /λ = 1 (triangles), N /λ = 3 (squares); (d) as in (b) for χ = − in future work. The 1D configuration, in optics, can beeither achieved by considering elliptical beams, or in aguiding slab geometry. In the former case, transverse in-stabilities may occur that break the 1D regime (as, e.g.,investigated in [33]). However, these instabilities are ab-sent for a defocusing nonlinearity, as in thermal liquids,or require very high power levels and a shaping of the in-put beam (e.g., to impose an unstable periodical modula-tion) to be relevant. Vectorial effects are specifically ab-sent when considering electrostrictive nonlinearities [34],or thermal focusing [35] and defocusing phenomena [36],which are also known to be highly nonlocal. Nonlocality,in turn, is known to filter out unstable mechanisms.[26]Light absorption sustains thermal effects but may alterthe nonlinear dynamics. However, as specifically true inthe non-paraxial regime with tightly focused wavepack-ets, the absorption length L abs can be much larger thanthe nonlinear and diffraction lengths L d . This fact en-ables to reach highly nonlinear regimes (as in [36]) withnegligible role of loss. For example, in a thermal liquid[36] with refractive index n ∼ = 1 .
3, and L abs ∼ = 2 mmmuch longer than L d for an elliptical beam focused inone transverse direction x at the wavelength scale. For λ = 532 nm, L d ∼ = 4 µ m << L abs , and the input beamcan be generated by high numerical aperture objectives.Due to the very pronounced nonlinear response of ther-mal liquids,[36] the QG-SW are expected at power levelsof the order of 100 m W. Conclusions —
In this manuscript, we have studiednonlinear waves in a model based on a generalizedSchr¨odinger equation implying a modified uncertaintyprinciple with minimal length uncertainty and maximalmomentum; the model is obtained by projecting in a for-ward spatial direction the Helmholtz equation. We haveshown that, for local and nonlocal responses, nonlinear-ity forces the solitary wavepackets to reach the maximallocalization. The reported results show that some of theconsequences of the generalization of quantum mechan-ics emerging from the physics at the Planck scale can benowadays tested in the laboratory by using optics andphotonics. Indeed, the theoretical methods used in thegeneralized quantum mechanics are naturally suited forstudying nonlinear optics beyond the paraxial approxi-mation, and may also be extended to tackle the case ofultra-short pulses without the slowly varying envelopeapproximation, or discrete systems.[37–40] These analy-ses may open the road to the first applications of quan-tum gravity physics, as in microscopy, spectroscopy, andnanophotonics, as well as to the use of quantum linearand nonlinear optics for the simulation of the physics atthe Planck scale.We acknowledge support from the Humboldt founda-tion and the CINECA award under the ISCRA initiative,for the availability of high performance computing. ∗ [email protected][1] A. Kempf, G. Mangano, and R. B. Mann, Phys.Rev.D , 1108 (1995).[2] H. Hinrichsen and A. Kempf, J.Math.Phys. , 2121(1996).[3] R. J. Adler and D. I. Santiago, Mod.Phys.Lett. A ,1371 (1999).[4] R. Akhoury and Y.-P. Yao, Phys.Lett. B , 37 (2003).[5] S. Das and E. C. Vagenas, Phys.Rev.Lett. , 221301(2008).[6] K. Sailer, Z. Peli, and S. Nagy, Phys. Rev. D , 084056(2013).[7] S. Das and E. C. Vagenas, Can.J.Phys. , 233 (2009).[8] L. N. Chang, D. Minic, N. Okamura, and T. Takeuchi,Phys. Rev. D , 125028 (2002).[9] Y. Chargui, L. Chetouani, and A. Trabelsi,Comm.Theor.Phys. , 231 (2010).[10] H. Yang, H. Miao, D.-S. Lee, B. Helou, and Y. Chen, Phys.Rev.Lett. , 170401 (2013).[11] I. Pikovski, M. R. Vanner, M. Aspelmeyer, M. S. Kim,and C. Brukner, Nat. Phys. , 393 (2012).[12] J. D. Bekenstein, Phys. Rev. D , 124040 (2012).[13] G. Amelino-Camelia, Phys. Rev. Lett. , 101301(2013).[14] P. Pedram, Phys.Lett. B , 317 (2012).[15] L. D. Landau and E. M. Lifsits, Quantum Elec-trondynamics , 2nd ed., Course of Theoretical Physics(Butterworth-Heinemann, Oxford, UK, 1982).[16] R. Gambini and J. Pullin, Phys.Rev. D , 124021(1999).[17] J. Alfaro, H. A. Morales-Tecotl, and L. F. Urrutia,Phys.Rev. D , 103509 (2002).[18] F. Mercati, D. Mazn, G. Amelino-Camelia, J. M. Car-mona, J. L. Corts, J. Indurin, C. Lmmerzahl, and G. M.Tino, Classical Quantum Grav. , 215003 (2010).[19] F. Briscese, M. Grether, and M. de Llano, Euro-phys.Lett. , 60001 (2012).[20] E. Castellanos, Europhys.Lett. , 40004 (2013).[21] A. Kempf, Phys. Rev. D , 024017 (2000).[22] E. Granot, S. Sternklar, Y. Isbi, B. Malomed, andA. Lewis, Opt.Comm. , 121 (1999).[23] A. Ciattoni, B. Crosignani, P. D. Porto, and A. Yariv,J. Opt. Soc. Am. B , 1384 (2005).[24] N. Bulso and C. Conti, Phys. Rev. A , 023804 (2014).[25] S. Detournay, C. Gabriel, and P. Spindel, Phys. Rev. D , 125004 (2002).[26] W. Krolikowski, O. Bang, N. I. Nikolov, D. Neshev,J. Wyller, J. J. Rasmussen, and D. Edmundson, Journalof Optics B: Quantum and Semiclassical Optics , S288(2004).[27] O. Bang, W. Krolikowski, J. Wyller, and J. J. Ras-mussen, Phys. Rev. E , 046619 (2002).[28] A. W. Snyder and D. J. Mitchell, Science , 1538(1997).[29] W. Krolikowski and O. Bang, Phys. Rev. E , 016610(2000).[30] V. Folli and C. Conti, Opt. Lett. , 332 (2012).[31] E. L. Ince, Ordinary Differential Equations , Dover Bookson Mathematics (Dover, New York, USA, 1956).[32] J. P. Boyd,
Chebyshev and Fourier Spectral Methods (Dover, New York, 2001).[33] V. Folli, E. DelRe, and C. Conti, Phys. Rev. Lett. ,033901 (2012).[34] R. W. Boyd,
Nonlinear Optics , 2nd ed. (Academic Press,New York, 2002).[35] S. Gentilini, F. Ghajeri, N. Ghofraniha, A. D. Falco, andC. Conti, Opt. Express , 1667 (2014).[36] S. Gentilini, N. Ghofraniha, E. DelRe, and C. Conti,Phys. Rev. A , 053811 (2013).[37] S. Longhi, Applied Physics B , 453 (2011).[38] D. Faccio, T. Arane, M. Lamperti, and U. Leonhardt,Classical and Quantum Gravity , 224009 (2012).[39] A. Marini, T. X. Tran, S. Roy, S. Longhi, and F. Bian-calana, Phys. Rev. A , 023840 (2014).[40] S. Finazzi and I. Carusotto, Phys. Rev. A89