Highly nonlocal optical nonlinearities in atoms trapped near a waveguide

Nonlinear optical phenomena are typically local. Here we predict the possibility of highly nonlocal optical nonlinearities for light propagating in atomic media trapped near a nano-waveguide, where long-range interactions between the atoms can be tailored. When the atoms are in an electromagnetically-induced transparency configuration, the atomic interactions are translated to long-range interactions between photons and thus to highly nonlocal optical nonlinearities. We derive and analyze the governing nonlinear propagation equation, finding a roton-like excitation spectrum for light and the emergence of order in its output intensity. These predictions open the door to studies of unexplored wave dynamics and many-body physics with highly-nonlocal interactions of optical fields in one dimension.

Optical nonlinearities are commonly described by local nonlinear response of the material to the optical field, resulting in the dependence of the refractive index at point z on the field at the same point, E(z) [1]. Recently however there has been a growing interest in nonlocal nonlinear optics, namely, in mechanisms whereby the refractive index at z depends on the field intensity at different points z in the material [2]. Mechanisms that give rise to such nonlocal nonlinearities include e.g. heat diffusion [3], molecular reorientation in liquid crystals [4] and atomic diffusion [5,6]. This paper discusses a new regime of extremely nonlocal and controllable optical nonlinearities in atomic media, which is achieved by designing the interactions between atoms prepared in an electromagnetically-induced transparency (EIT) configuration.
EIT is associated with the lossless and slow propagation of light pulses in resonant atomic medium subject to coherent driving of an auxiliary atomic transition [7]. Since the early days of EIT, it has been explored as a means of enhancing optical nonlinearities [7][8][9][10]. A particulary effective mechanism for giant optical nonlinearities is provided by dipolar interactions between atoms that form the medium: Since EIT can be described by the propagation of the so-called dark-state polariton [11], which is a superposition of the light field and an atomic spin wave, the inherently nonlocal dipolar interactions between atoms are translated to nonlocal nonlinearities in polariton propagation. In the case of dipolar interactions between Rydberg atoms in free space, most theoretical [12][13][14][15] and experimental [16][17][18][19][20] studies have focused on their remarkable strength, a useful feature in quantum information, whereas their nonlocal aspect has received less attention [21].
The present work rests on two previously explored mechanisms: that of EIT polaritons and that of modi-fied long-range dipolar interactions in confining geometries, such as fibers, waveguides, photonic band structures [22] or transmission lines, that may entail possible giant vacuum-induced forces [23] and long-range entanglement generation [24]. Yet, we show that the combined effect of these mechanisms may allow for unfamiliar highly nonlocal optical nonlinearity. Namely, purely dispersive laser-induced dipolar interactions between atoms coupled to a nano-waveguide with a grating, which can extend over hundreds of optical wavelengths [25,26], are shown to be translated via EIT into highly nonlocal optical nonlinearities.
To this end, we derive a nonlinear propagation equation for the (possibly quantum) EIT polariton fieldΨ along the fiber axis z, that corresponds to a general interatomic potential U (z −z ) and includes lowest-order nonadiabatic corrections, Here the frequency shift and quadratic dispersion coefficients α and C, and group velocity v (specified below) are all real. The effect of the nonlocal interaction appears in the nonlinear detuning operator,δ N L , which adds up to the detuning of the coupling field δ c (Fig. 1a), giving rise to a nonlocal and nonlinear frequency shift and quadratic dispersion. This paper has two linked objectives: first, to elucidate the conditions under which such nonlocal and nonlinear propagation regime may occur; and second, to reveal, in this regime, a roton-like excitation spectrum of optical waves, which we predict to arise under EIT in atomic media trapped near a fiber-Bragg grating. This roton-like narrow-band spectrum is a signature of the tendency of light in this regime to exhibit spatial self-order, namely,  1: (a) EIT atomic configuration: the probe fieldÊ is resonantly coupled to the |g → |e transition whereas the coupling field Ω is coupled to the |d → |e transition with detuning δc. Interaction between the atoms in the |d -level [see (b)] induces its energy shift δNL which is effectively added to the detuning δc. (b) Laser-induced dipolar interactions: the laser with Rabi frequency ΩL and detuning δL operates on the |d → |s transition of all atoms, |s being an additional level, thus inducing a dipolar potential U (z) between pairs of atoms (z apart) populating the state |d [42]. (c) Setup: atoms in the EIT configuration (black dots), illuminated by the fieldsÊ and Ω (thin blue arrow), are trapped at a distance ra from a nano-waveguide (gray cylinder) along z, from z = 0 to z = L. A far-detuned laser ΩL (thick orange arrow), tilted by an angle θL from the z axis, induces long-range interactions between the atoms, mediated by the waveguide modes.
crystal-like correlations, which can be revealed by measuring the photon intensity at the fiber's output. The predicted self-order of light constitutes a new, hitherto unexplored, optical "phase", analogous to the spatial structure of cold atomic media subject to laser-induced dipole-dipole interactions [25,[27][28][29][30][31].

The system
Consider a medium of identical atoms in an EIT configuration as in Fig. 1: The atoms are trapped at a distance r a from a nano-waveguide along its longitudinal z axis [32][33][34][35][36][37][38] (Fig. 1c). A strong (external) coupling field with Rabi frequency Ω drives the |d → |e atomic transition with detuning δ c and wavenumber k c (Fig. 1a), whereas a weak (possibly quantized) field with carrier frequency ω 0 , wavenumber k 0 and envelopeÊ is resonantly coupled to the |g → |e transition (ω 0 = ω eg ). Under tight transverse trapping (around r a ) with respect to the transition wavelength, the atomic positions can be characterized by their longitudinal component z [25,32]. We assume the existence of a dipolar interaction U (z) between atoms in the state |d . Then, the energy of level |d of an atom at z is shifted by δ N L ∼ n a dz U (z − z )P d (z ), n a being the atomic density (per unit length) and P d (z ) the occupation of state |d in an atom at z .
We may now explain the physical reasoning that leads to Eq. (1). In Fig. 2 we plot the complex linear sus- FIG. 2: Linear susceptibility of the EIT atomic medium to the probe field as a function of its detuning ∆p [7]. (a) For total detuning of the coupling field ∆c = 0, the absorption Imχ (red dashed line) and dispersion Reχ (blue solid line) are symmetric and antisymmetric, respectively, with respect to ∆p, so that no (real) quadratic dispersion exists for the probe envelope centered around ∆p = 0. (b) For ∆c = 0, Reχ is not antisymmetric, so that quadratic dispersion exists, giving rise to the term ∆cCv 2 ∂ 2 z in Eq. (1), with ∆c = δc + δNL (Fig.  1a).
ceptibility χ of the EIT medium to the probe field as a function of its detuning ∆ p , in the presence of a coupling field detuned by ∆ c [7], which in our case is given by ∆ c = δ c + δ N L (Fig. 1a). When ∆ c = 0 (Fig. 2a), the absorption coefficient Imχ is symmetric with respect to ∆ p whereas the dispersion Reχ is antisymmetric, so that no (real) quadratic dispersion exists for the probe envelope centered around ∆ p = 0. By contrast, when ∆ c = 0 is introduced (Fig. 2b), Reχ is no longer antisymmetric and quadratic dispersion exists, which explains the term ∆ c Cv 2 ∂ 2 z in Eq. (1). However, this comes at the price of non-vanishing losses at ∆ p = 0. For this reason we choose to work in the so-called Autler-Townes regime [7] Ω γ, ∆ c , where γ is the width of the level |e . Then, for ∆ c smaller than the single-atom transparency window, ∆ c Ω 2 /γ, but still larger than γ, the absorption per atom can become negligible while dispersion is still significant, as illustrated in Fig. 2b. This explains the lossless propagation described by Eq. (1) with real parameters α, v, C. As long as the absorption, associated with dissipation due to spontaneous emission at rate γ, is negligible, so are the noise effects of vacuum fluctuations; Eq. (1) then holds in operator form without additional Langevin quantum noise operators.
The formal derivation of Eq. (1) is sketched in the Appendix. Here we merely point out that the polariton fieldΨ is given by the superpositionΨ(z) = cos θÊ(z) − √ n a L sin θσ gd (z), where n a is the number of atoms per unit length in the medium,σ gd (z) is the atomic |g ↔ |d spin at z and θ is the EIT mixing angle. The lossless propagation regime governed by Eq. (1) with α = sin 2 θ, v = cos 2 θc, C = sin 2 θ(2 − 3 sin 2 θ)/|Ω| 2 is obtained by assuming all detunings to be smaller than the EIT transparency window δ tr = Ω 2 /(γ √ OD), where OD is the optical depth of the atomic medium.

Excitation spectrum of polariton waves
Let us assume a continuous wave (CW) polariton, and find the dispersion relation of wave excitations around this CW background, in analogy to the Bogoliubov spectrum of excitations in a Bose-Einstein condensate (BEC) [39]. The CW solution of (1) is ψ(t) = ψ 0 e −i(αδc+npU0)t , with ψ 0 = |ψ 0 |e iφ , n p = α 2 |ψ 0 | 2 being the effective photon density per unit length and U 0 the k = 0 component of the spatial Fourier transform of the potential U k = ∞ −∞ dzU (z)e −ikz . Here we have neglected edge effects by assuming l < z < L − l, l being the range of the potential U (z). The dispersion relation of small fluctuations ϕ(z, t) around the large average CW field Ψ = ψ(t) are found as usual upon inserting Ψ = ψ(t) + ϕ(z, t) into Eq. (1) and linearizing it by keeping the fluctuations ϕ to linear order. Then, introducing the ansatz [39], into the linearized equation for ϕ, u k and v k being cnumber (Bogoliubov) coefficients, and using standard procedures [39], we find the modified Bogoliubov spectrum [see Supplementary Information (SI)] (3) This means that a polariton wave distortion (about the CW solution) with a wavenumber k relative to the carrier wavenumber (inside the EIT medium), oscillates at a frequency (relative to the carrier frequency ω 0 ) The dispersion relation (4) is composed of the detuning αδ c due to the coupling field, the self-phase modulation of the CW component n p U 0 (analogous to the chemical potential in a BEC [39]), the linear dispersion vk, and the modified Bogoliubov excitation spectrum ω k . The spectrum ω k is determined by the interplay between the interaction Fourier transform U k and the "kinetic-energy" quadratic dispersion ω 0 k , which is affected by both the detuning δ c and by the k = 0 component of U k . This interplay is further discussed below for U (z) resulting from laser-induced interactions near a waveguide grating.

Generation of two-mode correlations
The parametric process described by the foregoing modified Bogoliubov theory also entails the dynamic generation of two-mode squeezing, i.e. pairs of entangled polaritons with wavenumbers ±k. The analysis is similar to that of propagation in fibers with local Kerr nonlinearity [40,41]. Upon inserting the expansion of small quantum fluctuations in the longitudinal wavenumber modes k,φ(z) = kâ k e ikz / √ L, into the linearized equation for ϕ(z, t), we obtain coupled first-order differential equations (in time) forâ k (t) andâ † −k (t), whose solution is a dynamic Bogoliubov transformation (SI), The number of entangled pairs generated after propagation time t at wavenumbers ±k can be quantified by the so-called squeezing spectrum, whose optimum is given by (all fields measured at the fiber's output after propaga- where the averaging is performed with respect to the initial probability distribution, e.g. the quantum state, of the polariton (probe) field. For initial zero-mean fluctuations (around the CW solution) with average polariton occupation at mode k, â † k (0)â k (0) = N k δ kk , and vanishing anomalous correlations â k (0)â k (0) = 0, we find (SI) where φ k = arg(µ k ν k ) and we used (1/L) k → dk/(2π). As we shall see below, these correlations may reveal the ordering of a nonlocal system caused by pairgeneration at preferred k-values.

Laser-induced interaction via waveguide grating
We now turn to the specific case of laser-induced interatomic dipolar interaction [42]. Consider another laser with Rabi frequency Ω L which is detuned by |δ L | Ω L from the transition |d → |s , |s being a fourth atomic level (Fig. 1b), and assume that this transition is distinct and separated from the transitions used for EIT ( Fig.  1a) either spectrally or by polarization. We further take the laser to be oriented at an angle θ L with respect to the waveguide (Fig. 1c), and that the nano-waveguide incorporates a grating, i.e. periodic perturbation of the refractive index with period Λ ≡ π/k B (see e.g. Refs. [35,43]). Then, for ω L inside and close to the edge of the photonic bandgap associated with the grating (the probe field's carrier frequency ω 0 being outside the gap), the laser-induced interaction potential between the atoms becomes [25] where l can extend over hundreds of wavelengths [25].
Here k z L = k L cos θ L with k L the laser wavenumber, and U L depends on the laser parameters, atomic transition and effective area at r a (see Appendix). The resulting spatial Fourier transform U k then consists of four Lorentzian peaks of width ∼ 1/l, centered around the spatial beating frequencies ±(k z L − k B ) and ±(k z L + k B ).

Roton and Anti-Roton spectra
Let us focus on the peak of U k around k R ≡ k z L − k B and its effect on the dispersion relation (spectrum), Eq.
(3). We first consider the case of anomalous dispersion, where the signs of ω 0 k and U k are opposite. In analogy to BEC, this describes the case of an attractive potential U k that competes with the "kinetic energy" ω 0 k . Then, for k-values satisfying |ω 0 k | > 2n p |U k |, ω k is real and exhibits a dip around k R , in contrast to the case of a local potential for which U k is independent of k (standard Bogoliubov spectrum) and this feature is absent. This is seen in Fig. 3a, where both the analytical results of Eq.
(3) and numerical simulations of the nonlinear equation (1) (Appendix), are plotted and shown to agree very well. The narrow-band "dip" of this ω k spectrum is in analogy with the roton minimum in He II [44]. It reflects the fact that wave distortions about the CW field with spatial frequencies around k R cost less energy and are hence favorable. This feature implies that the intensity of the polariton field in its ground state would tend to self-order with typical wavenumber k R [44].
Turning to the case of normal dispersion, where the signs of ω 0 k and U k are identical, ω k exhibits an "antiroton" peak around k R (Fig. 3b). This means that distortions of spatial frequencies around k R are costly, so that the system prefers to avoid these spatial variations. This behavior again manifests the tendency of the system to order, since it indicates the spatial distortions that the system is unlikely to be found in, should it be in its ground state.
In order to measure the roton and anti-roton spectra, we first recall the meaning of the dispersion relation ω(k) from Eq. (4): without an interaction, the frequency associated with a wave envelope at wavenumber k traveling inside the EIT medium is given by , so that ω k (together with n p U 0 ) expresses a frequency, or phase velocity, shift due to the nonlinearity. Namely, the wavenumber k describes a spatial eigenmode of propagation, both with and without interaction, with an eigen- (1) (gray dots, see Appendix). Compared with the spectrum of a local interaction (U k independent of k, dashed red line), the roton-like spectrum exhibits a dip in a narrow band of k-values around kR = k z L − kB. (b) Anti-roton peak of the spectrum ω k in a narrow band around kR in the normal dispersion case (identical signs of ω 0 k and U k ). (c) Possible homodyne detection scheme: the input probe field consists of a CW field + perturbation at wavenumber k and frequency ω (0) (k). The field is split before entering the EIT medium (z = 0), so that a local oscillator of ω (0) (k) is formed (lower arm) by e.g. filtering out the CW component. Then, mixing the output signal (z = L) with the local oscillator reveals their phase difference, from which ω k can be inferred (see text).
frequency ω(k) and ω (0) (k), respectively. Now, suppose we let a weak quasi-CW pulse of length L p < L and frequency ω p = ω (0) (k) enter the medium (on top of the strong CW), when the laser and hence the interaction U k are turned off. Since upon entering the medium the field does not change its frequency ω p , we deduce from the dispersion relation in the absence of interaction, ω (0) (k) that the field inside the medium exhibits a perturbation ϕ(z) at spatial frequency k on top of the strong CW. Subsequently, when the entire pulse is in the medium, we immediately (non-adiabatically) turn on the laser Ω L , and hence the interaction U k , so that the temporal frequency of the perturbation ϕ(z) at wavenumber k becomes ω(k) = ω (0) (k) + n p U 0 + (ω k − δ c Cv 2 k 2 ). The frequency shift n p U 0 + (ω k − δ c Cv 2 k 2 ) can be thought of as an extra energy acquired by the mode k due to the interaction energy. Therefore, the spectrum ω k can be inferred from the frequency shift, measurable by homodyne detection of the pulse ϕ(z) that exits the EIT medium (Fig. 3c). A similar procedure was proposed for measuring the tachyon-like spectrum of polaritons in inverted media [46].

Dynamical instability: pair generation
In the anomalous dispersion case, consider now a sufficiently strong interaction such that for a narrow band of k-values around k R , where U k is peaked, the condition 2n p |U k | > |ω 0 k | is satisfied and ω k becomes imaginary. Then, field perturbations around k R become exponentially unstable (Fig. 4a), resulting in parametric amplification and generation of entangled photon pairs in this narrow band of unstable k-values. The strength of the amplified perturbation and generated field is characterized by the magnitude of the coefficients µ k (t) and ν k (t) from Eq. (5), which grow exponentially with propagation time t and are largest for the narrow peak around k R . The resulting squeezing spectrum G k at the output t = L/v (Fig. 4a) may be measured by homodyne detection [40,45].

Dynamical instability: emergence of self-order
Perhaps the most remarkable implication of the extremely nonlocal potential of Eq. (7) is the dynamic formation of order in the system. Consider that at t = 0 there exist fluctuations around the CW, with a spatial spectrum N k = N 0 e −(k/qtr) , i.e. "δ-correlated" noise limited by the EIT transparency window of width δ tr = vq tr . Then, fluctuations at k-values around the peak k R will be parametrically amplified as they propagate through the medium, as verified in Fig. 4b, where the spatial spectrum of the polariton field, N k (t) = â † k (t)â k (t) , at different propagation times t (t = L/v describing the output field) is calculated both analytically and via direct (classical) numerical simulations of Eq. (1). This suggests that the system becomes spatially ordered with a period ∼ 2π/k R , which may be revealed by measuring the correlation function g (2) between the intensities of the output field that arrive at a detector D at the waveguide's end (z = L) and time difference (z − z )/v (scheme from Fig. 3c without the lower local-oscillator arm). We obtain the corresponding g (2) by numerically integrating over k in the classical limit of Eq. (6), where the vacuumfluctuations contributions are neglected. This calculation is compared to g (2) measured via direct numerical simulations of Eq. (1), yielding excellent agreement. Fig. 4c reveals the emergence of order by presenting the intensitycorrelations g (2) at different propagation times t through the medium (t = L/v at the output): The resulting g (2) exhibits oscillations with a period z − z ∼ 2π/(k R ), which persist over a few l. Considering that the interaction range l can reach hundreds or even thousands of optical wavelengths (l ≈ 3027λ L and L ∼ 10l ∼ 0.026 (a) γ k = Imω k is the exponential growth rate of unstable perturbations at spatial frequency k (anomalous dispersion case). It exhibits a narrow peak around kR (blue dotted line), compared to the broadband instability of the local-interaction case (red dotted line). The instability is accompanied by generation of quantum entanglement, characterized by a narrow-band squeezing spectrum G k (blue solid line; G k < 1 quantifies entanglement between ±k photon modes), in contrast to the broadband spectrum for a local interaction (horizontal red solid line). (b) Dynamics of spatial spectrum N k (t) of the field intensity for different propagation times t inside the medium (t = L/v at the output). The emergence of a large peak around kR (and −kR) out of the initial Gaussian perturbation is clearly seen. Results of both the linearized theory (solid lines) and numerical simulations of Eq. (1) (dots) are shown: slight differences at later t-values are attributed to nonlinear corrections. (c) Self-ordering of the field intensity: considering the narrow peaks of N k (t) around ±kR, fluctuations at these k-values become dominant, resulting in ordered intensity correlations g (2) which grow with propagation time t (t = L/v at the output). Excellent agreement between the theory, Eq. (6), (solid lines) and numerical simulations of Eq. (1) (dots) is observed. The correlations oscillate with a period 2π/kR ∼ 6.1 mm and a range of a few l ≈ 2.7 mm (where L = 2.68 cm and v = 4340 ms −1 ), determined by the long-range interaction U (z). This is in contrast to the local-interaction case, where the correlations vanish (g (2) = 1) after a short distance π/qtr ≈ 1.75 mm, which is determined by the bandwidth qtr of the initial fluctuations [see inset: local (red thin line) and nonlocal (blue thick line) cases for the output field at t = L/v]. m in our example, see Appendix), the light intensity clearly becomes ordered due to the long-range interaction U (z), as can also be seen by the comparison to the local-interaction case (inset of Fig. 4c).

Imperfections and scattering
So far we have considered a purely coherent evolution of the polariton field. Here we address three main sources of scattering and loss of field excitations (see also SI). First, the interaction we consider involves the illumination of the atoms by an off-resonant laser, which is accompanied by an incoherent process of scattering of laser photons Ω L from the |d → |s transition to non-guided modes at rate R f s [42]. This process limits the coherence time of theσ gd spin wave and hence that of the polariton to be below ∼ R −1 f s , which in the example of Figs. 3 and 4 is nevertheless much longer than the experiment time L/v (see Appendix).
Secondly, material imperfections and defects in the waveguide grating structure may give rise to scattering of photons off the guided mode. This effect is analyzed for dipolar interactions in Ref. [26], in terms of the socalled cooperativity of the waveguide structure, P , whose square root √ P gives the ratio between coherent (dipolar interactions mediated by waveguide modes, U ) and incoherent processes (scattering to non-guided modes via the imperfections). For a realistic waveguide material of length λ L , P ∼ 10 4 [26]. Since the imperfection-mediated incoherent processes accumulate with the system size, and considering the dipolar interaction range as the effective system length l ∼ 3000λ L (Appendix), we estimate P ∼ 10 4 λ L /l ∼ 3.3, giving rise to an imperfectionmediated scattering rate of R im ∼ |U L |/ √ 3.3 per atom. The effect of R im on the spectrum ω k and the observables mentioned above can in principle be made arbitrary small however, by noting that ω k depends on n p U L whereas R im depends solely on U L . Then, decreasing |U L | while keeping n p U L constant (by increasing the CW power n p ) reduces R im but keeps ω k unchanged. This reflects the fact that R im is a single-polariton loss mechanism and hence independent of n p , whereas the effects discussed above are cooperative.
Finally, we consider losses and decoherence of spinwave excitations of the EIT medium due to the nonvanishing detuning of the coupling field, ∆ c = δ c + δ N L , where here δ N L = n p U 0 /α is the mean-field value ofδ N L [Eq. (1)] that appears in the linearized equation for ϕ. The transmission through the EIT medium of length L is given by e −(1/2)k0χ L , leading to an imaginary frequency −i(1/2)k 0 χ v to be added to ω k , χ = Imχ(∆ p = 0, ∆ c ) being the imaginary part of the EIT susceptibility [7] evaluated for simplicity at the center of the probe pulse (∆ p = 0) and in the presence of the coupling field detuning ∆ c = δ c + δ N L . In the SI we show that the effect of the EIT losses on ω k can be significantly decreased by reducing both n p and δ c .
In principle, all loss terms discussed above should be accompanied by corresponding quantum noise terms, which are not considered here, restricting this discussion to the classical-field case.

Conclusions
This study predicts a new and hitherto unexplored regime of nonlinear optics; namely, that of highly nonlocal interactions between photons in one dimension (1d). These nonlocal optical nonlinearities arise for light propagation inside driven atomic media in the vicinity of a waveguide. We have derived the nonlinear equation that governs light propagation in this regime, and have analyzed it around its CW solution, finding a narrow rotonlike dispersion relation (Figs. 3a,b) and squeezing (entanglement) spectrum (Fig. 4a), and the emergence of order in the field-intensity (Figs. 4b,c), all of which reflect the tendency of the system to self-organize, which in turn results from the long-range interactions between photons.
We stress that the length-scale associated with order, k R = k L cos θ L − k B is not imposed solely by the external grating with period Λ = π/k B , but rather depends also on the wavenumber of the laser k L and its orientation θ L , reflecting the origin of order in the self-interaction of the light field. This is in contrast to order in e.g. an optical lattice, where the atoms are situated at lattice sites determined by the potential imposed by an external laser.
This work opens the way to experimental and theoretical investigations of new nonlinear wave phenomena, especially by venturing beyond the linearized regime and exploring the role of the nonlinear dispersion term δ N L Cv 2 ∂ 2 z . Specific directions of further research may include: (1) Nonlocal nonlinear optics in 1d, concerning the study of solitons, where the lensing effect created by the nonlinear refractive index change is now highly nonlocal, following U (z). (2) Thermalization in 1d. This important issue has been considered both experimentally [48] and theoretically [49] for isolated BEC in 1d. Our proposed scheme may give access to important modifications in the possible route to thermalization affected by the nonlocal character of the designed interactions U (z) and the nonlinear dispersion ("mass"). (3) Effects of non-additivity of systems with long-range interactions [50] may be studied here for quantum/classical optical fields. (4) Nonlocal interactions in other confining geometries, e.g. cavities, may give rise to other unfamiliar nonlocal optical nonlinearities. (5) The intriguing possibility of measuring field amplitudes and not only intensities, via homodyne detection, provides a possible advantage over the study of quantum fields in matter waves (atoms in a BEC), where measurements using a local oscillator are less straightforward. Equation (1) is derived as follows. The field enve-lopeÊ(z) = kâ k e ikz / √ L, with commutation relations [â k ,â k ] = δ kk and hence [Ê(z),Ê † (z )] = δ(z − z ), is assumed to be spectrally narrow and is guided by a transverse mode of the fiber with effective area A at position r a and polarization vector e 0 . The Hamiltonian in the interaction picture is (8) and H F = k ckâ † kâ k , where g = ω 0 /(2 0 A)d · e 0 , d being the dipole matrix element of the |g → |e transition, andσ ij (z) = |i j| for an atom at z, with i, j representing the states {g, d, e}. We write the Heisenberg equations for the operatorsÊ(z),σ ge (z) andσ gd (z) and keep only linear terms inÊ(z), including nonlinearity solely through the potential U (z − z ) [13]. We then transform them into equations of motion for the dark and bright polaritons,Ψ andΦ, respectively, where tan 2 θ = n a g 2 /|Ω| 2 andσ gd (z) = σ gd (z)e i(kc−k0)z e −iδct . Next, we take the first nonadiabatic correction to EIT [11] by inserting the equation of motion forΦ into that ofΨ, and assume all detunings to be smaller than the EIT transparency window δ tr = Ω 2 /(γ √ OD), with OD = (n a /A)Lσ a and σ a the cross section of the |g → |e transition, finally obtaining Eq. (1).

Numerical simulations of Eq. (1)
In order to obtain the ω k spectrum, we perform numerical simulations of the full nonlinear equation (1) using a three-term splitting method, with an initial field comprised of the CW solution ψ and weak perturbations ϕ at a spatial frequency k. From the dynamics of the field, we can extract the temporal oscillation frequency of the field ω(k) leading to the spectrum ω k . The results of the simulation of Eq. (1) agree very well with those of a simpler split-step method simulations, where the nonlinear dispersion termδ N L Cv 2 ∂ 2 z is approximated by its mean-field value, δ N L = n p U 0 /α. For the case of instability (Fig.  4b,c), we begin with the initial intensity fluctuations with a spectrum N k and run simulations with the mean-field nonlinear dispersion term δ N L = n p U 0 /α. Then, after propagation time t, we measure N k (t) = a * k (t)a k (t) and g (2) = I(z)I(z ) /( I(z) I(z ) ), with I(z) = |Ψ(z)| 2 for a classical field Ψ.
We acknowledge discussions with Ofer Firstenberg, Arno Rauschenbeutel, Philipp Schneeweiß, Darrick Chang, Nir Davidson and Tommaso Caneva and the support of ISF, BSF and FWF [Project numbers P25329-Here we wish to present the derivation of the polariton wave excitation spectrum (dispersion relation) of perturbations around the CW solution, Eq. (3) in the main text. We first write the polariton field as a sum of an average CW Ψ (z, t) = ψ(t) and a small fluctuationφ(z, t),Ψ where the average CW solution is that of self-phase modulation, ψ(t) = ψ 0 e −i(αδc+npU0)t with ψ 0 = n p /α 2 e iφ , U 0 = U k=0 , U k being the spatial Fourier transform of the potential U(z), The integral for U 0 that arises in the CW solution runs over z = 0 to z = L around some point z, L 0 dz U (z − z ), and not from z = −∞ to z = ∞ around z = 0 as in the definition (S2). However, for a symmetric potential [U (z) = U (−z)] with range l, and points z well within the medium, i.e. l < z < L − l, the edges z = 0 and z = L have no effect, so that U 0 can be approximated as for all such points z. Thus, we may neglect the edge effects of a sufficiently long medium L l. Inserting Eq. (S1) into the nonlinear equation (1) in the main text and keeping terms only to linear order in the perturbationφ, we obtain This linearized equation is valid as long as the fluctuationsφ around the mean ψ are small, i.e. for In order to find the dispersion relation of the fluctuations it is enough to consider the case of a classical fieldφ → ϕ. We then insert the ansatz u k and v k being c-number coefficients, into the linearized equation (S3), obtaining coupled equations for u k and v k In order to arrive at the equations (S6), we used which again amounts to neglecting edge effects, and where in the last equality we assumed that U (z) is symmetric and real. Eq. (S6) has a nontrivial solution only if the determinant of the matrix on the left-hand side is zero, yielding an extra equation whose solution is the spectrum ω k from Eq. (3).

Dynamic Bogoliubov theory for a nonlocal interaction
In this part we address the quantum description of the dynamics of the fluctuationsφ around the CW solution and arrive at the dynamical Bogoliubov transformation with the coefficients of Eq. (5) in the main text. We expand the quantum fieldφ(z, t) in spatial Fourier modesâ k (t),φ(z) = k (1/ √ L)e ikzâ k , and the transformed modesĉ k (t) as, where the operatorsâ k (t) satisfy the equal-time commutation relations [â k (t),â † k (t)] = δ kk and hence so do the operatorsĉ k (t). Inserting Eq. (S8) into the linearized equation (S3) and neglecting edge effects as in Eq. (S7), we obtain an equation for ∂ tĉk . Then, upon taking its Hermitian conjugate, we end up with coupled equations of motion forĉ k andĉ † −k , By diagonalizing the matrix, we find the solution forĉ k (t) from which we obtain the dynamics ofâ k (t) as the Bogoliubov transformation from Eq. (5) in the main text.

Intensity correlations
Intensity correlations are characterized by the normalized second order coherence function defined by [1] Here we assume that all points of the fieldΨ(z) have experienced the same duration of interaction t upon arrival at the detector. The detector measures at different times different points of the fieldΨ(z), so that the correlations between detector-signals at different times are in fact correlations of the field in space as in g (2) (z, z ) which quantifies the autocorrelation of the field between z and z . Upon expanding the field in its longitudinal Fourier modes, we recall that the k = 0 mode can be approximated by the strong average CW solution, whereâ k (t) is given by Eq. (5) of the main text. Inserting Eq. (S12) into G (1) (z) from Eq. (S11), we find where we have assumed the following statistics of the initial fluctuationsâ k =0 : â k (0) = 0, â † k (0)â k (0) = N k δ kk with N k = N −k , and â k (0)â k (0) = 0. Moving to the second-order coherence, we insert Eq. (S12) into G (2) from Eq. (S11), keeping terms only to second order in the fluctuationsâ k =0 , in accordance with the assumptions (S4), and obtain (S14) where φ k = arg(µ k ν k ). Inserting G (1) and G (2) from Eqs. (S13) and (S14) in Eq. (S15), we note that to lowest order in the fluctuations we can take G (1) ≈ n p /α 2 in the denominator, arriving at (S15) Finally, upon taking the continuum limit (1/L) k → (1/2π) ∞ −∞ dk, we split the integral into its positive and negative k-values and obtain Eq. (6) from the main text.

Imperfections and scattering
The scattering and loss mechanisms discussed in the main text represent single-atom (single-polariotn) effects which may add a loss term of the type −γ lossΨ (z) on the right-hand side of the propagation equation, Eq. (1), and hence an imaginary part to the spectrum ω k , −iγ loss . In Fig. 1a, we plot ω k for the roton and anti-roton cases (blue solid line), compared with −γ loss for the three types of losses mentioned in the main text: 1. Scattering of offresonant laser photons to non-guided modes, γ loss = R f s (gray dashed line); 2. Imperfection-mediated scattering, γ loss = R im ∼ |U L |/ √ 3.3 (orange dashed line); 3. EIT losses due to detuning of the coupling field, γ loss = R EIT = −i(1/2)k 0 vImχ(∆ p , ∆ c ) (black dashed line). Considering the susceptibility of an EIT medium [2], we find, γ being the width of the atomic level |e and OD the optical depth of the medium (see main text). In all figures here, R EIT is evaluated at the center of the pulse (∆ p = 0) and for ∆ c = δ c + δ N L , δ N L = n p U 0 /α being the mean-field value of the nonlinear detuning.
In the roton case in Fig. S1a, using the same physical parameters from the main text ( Fig. 3a and Appendix therein) we observe that the off-resonant scattering R f s (gray dashed line) is negligible whereas R im and R EIT can become comparable to the ω k spectrum and may hence significantly alter the observable effects discussed in the main text. However, as explained in the main text, the effect of R im can in principle be made arbitrary small by decreasing |U L | and keeping U L n p constant, leaving us with R EIT (black dashed line) as a dominant loss mechanism. Reexamining the expression for ω k , Eq. (3) in the main text, we observe that a reduction of both δ c and n p by a factor f , results in a reduction of ω k by the same factor. This linear scaling is in contrast with the nonlinear dependence of R EIT with ∆ c as per Eq. (S16). This scaling difference is an example of a possible method for obtaining coherent effects (represented by ω k ) which are dominant over incoherent EIT losses, as shown in Fig. S1b, where f = 10 is taken.
Proceeding to the anti-roton case in Fig. S1c, we plot −γ loss for the three types of incoherent processes using the physical parameters from the main text ( Fig. 3b and Appendix therein). It is seen that both the EIT losses (black dashed line) and those of off-resonant scattering R f s (gray dashed line) are negligible, leaving us with the imperfection-induced scattering, which in turn can be drastically reduced as explained above.
Finally, considering the case of an instability, where in the lossless case we have Imω k = γ k > 0, losses give negative contribution to γ k , hence degrading the parametric amplification associated with the instability. Here R f s still gives a negligible contribution to γ k , as well as the effect of imperfections: in Fig. S2a we plot the lossless γ k (blue solid line) compared with γ k − R im (black dashed line) obtaining very similar curves. Considering the EIT losses, in Fig.  S2b we plot γ k (blue solid line) compared with γ k − R EIT and observe a significant reduction of the amplification rate. Nevertheless, reducing both δ c and n p by e.g. a factor f = 10 as in the roton case, may, in principle result in EIT losses which are negligible with respect to γ k , as presented in Fig. S2c.