Four-wave mixing and enhanced analogue Hawking effect in a nonlinear optical waveguide

TL;DR

This paper demonstrates that nonlinear optical waveguides with solitons can produce observable analogue Hawking radiation through four-wave mixing, with resonance effects amplifying photon pair production significantly.

Contribution

It provides a detailed analysis of photon pair production in nonlinear waveguides without the rotating wave approximation, revealing resonance effects that enhance the analogue Hawking effect.

Findings

Photon pairs are spontaneously emitted at a rate of about one per centimeter for a 10 fs soliton.

Resonance effects related to modulation instability amplify the scattering coefficients.

The results are confirmed by numerical simulations of the full nonlinear field evolution.

Abstract

We study in detail the scattering of light on a soliton propagating in a waveguide which has been proposed as an experimental system in which one could observe the analogue Hawking effect. When not applying the rotating wave approximation, we show that the linearized wave equation governing perturbations has the same structure as that governing phonon propagation in an atomic Bose condensate. By taking into account the full dispersion relation, we then show that the scattering coefficients encoding the production of photon pairs are amplified by a resonance effect related to the modulation instability occurring in the presence of a continuous wave. When using a realistic example of a silicon nitride waveguide on a silica substrate, we find that a soliton of duration 10 fs would spontaneously emit about one photon pair for every cm it travels, which makes the effect readily observable. This result is confirmed by numerically solving the equation encoding the Kerr nonlinearity and governing the evolution of the full field (soliton plus perturbations). We discuss the link with previous works devoted to the analogue Hawking effect where the pair creation rates were about six orders of magnitude smaller.

Equations

$$E(t,z)=A(t,z)\exp \left(i{\beta }_{0}z-i{\omega }_{0}t\right),$$

$$-i{\partial }_{z}A=\left(\beta \left({\omega }_0+i{\partial }_t\right)-{\beta }_0\right)A+\gamma {\left|A\right|}^{2}A.$$

$$\beta \left({\omega }_0+i{\partial }_t\right)-{\beta }_0=\sum _{n=1}^{\infty }\frac{{\beta }_n}{n!}{\left(i{\partial }_t\right)}^n,$$

$$-i\left({\partial }_z+{\beta }_1{\partial }_t\right)A=\left(-\frac{1}{2}{\beta }_2{\partial }_t^2+...\right)A+\gamma {\left|A\right|}^{2}A.$$

$$-i{\partial }_{z}A=B\left(i{\partial }_{\tau }\right)A+\gamma {\left|A\right|}^{2}A,$$

$$B\left(\Delta \omega \right)$$

$$N_A=⟨A|A⟩\doteq \int \mathrm{d}\tau {\left|A(\tau ,z)\right|}^2,$$

$$K_A\doteq \int \mathrm{d}\tau \left\{A^{\ast }(\tau ,z)B\left(i{\partial }_{\tau }\right)A(\tau ,z)+\frac{\gamma }{2}{\left|A(\tau ,z)\right|}^4\right\}.$$

$$-i{\partial }_z{A}_0=-\frac{1}{2}{\beta }_2{\partial }_{\tau }^2{A}_0+\gamma {\left|A_0\right|}^2{A}_0,$$

$$A_0(z,\tau )=\sqrt{P_0}\mathrm{sech}\left(\frac{\tau }{{\tau }_0}\right)e^{i\delta {\beta }_{0}z},$$

$$\delta {\beta }_0=\frac{1}{2}\gamma P_0=\frac{-{\beta }_2}{2{\tau }_0^2}>0.$$

$$A(z,\tau )=e^{i\delta {\beta }_{0}z}\left(A_0(\tau )+\delta A(z,\tau )\right),$$

$$-i{\partial }_z\left(\delta A\right)=B(i{\partial }_{\tau })\delta A-\delta {\beta }_0\delta A+2\gamma |A_0{|}^2\delta A+\gamma A_0^2\delta A^{\star }.$$

$$-i{\partial }_{z}w=Ow,$$

$${O}=\left[\begin{array}[]{cc}B(i\partial_{\tau})-\delta\beta_{0}+2\gamma|A_{0}(\tau)|^{2}&\gamma A_{0}^{2}(\tau)\\ -\gamma\left(A_{0}^{\star}(\tau)\right)^{2}&-B(-i\partial_{\tau})+\delta\beta_{0}-2\gamma|A_{0}(\tau)|^{2}\end{array}\right]\,.$$

$$\left(w_1,w_2\right)\equiv {\int }_{-\infty }^{+\infty }d\tau \left(w_{1,+}^{\star }(\tau )w_{2,+}(\tau )-w_{1,-}(\tau )w_{2,-}^{\star }(\tau )\right).$$

$${\partial }_z\left(w_1,w_2\right)=0.$$

$$D_{+}^{\mathrm{as}}(\Delta \omega )\doteq B(\Delta \omega )-\delta {\beta }_0.$$

$$D_{-}^{\mathrm{as}}(\Delta \omega )\doteq -B(-\Delta \omega )+\delta {\beta }_0.$$

$$w_{D,j}^{in/out}(\tau )$$

$${\bar{w}}_{D,j}^{in/out}(\tau )$$

$$\left(w_{D,i}^{in/out},w_{D^{\prime },j}^{in/out}\right)$$

$$\left({\bar{w}}_{D,i}^{in/out},{\bar{w}}_{D^{\prime },j}^{in/out}\right)$$

$$\left(w_{D,i}^{in/out},{\bar{w}}_{D^{\prime },j}^{in/out}\right)$$

$$\hat{w}(z,\tau )$$

$$w_{D,i}^{\mathrm{in}}$$

$${\bar{w}}_{-D,i}^{\mathrm{in}}$$

$$\sum _{j\in J_{+}(D)}{\left|{\alpha }_{D,ij}\right|}^2-\sum _{j\in J_{-}(D)}{\left|{\beta }_{D,ij}\right|}^2=1$$

$$k_{\mathrm{black}}(D)$$

$$k_{\mathrm{blue}}(-D)$$

$${\left|R_D\right|}^2+{\left|T_D\right|}^2\approx 1$$

$${\left|{\alpha }_D\right|}^2-{\left|{\beta }_D\right|}^2\approx 1.$$

$$\frac{{\partial }^{2}N}{\partial D\partial z}=\frac{1}{2\pi }{\left|{\beta }_D\right|}^2.$$

Full text

I Introduction

In 1981, William Unruh pointed out that an analogue version of the Hawking effect could be observed in a stationary flow of a moving medium when the flow velocity crosses the speed of quasiparticle excitations Unruh (1981). Since then, several proposals have been made in various media, including atomic Bose-Einstein condensates (BEC) Garay et al. (2000) and water in a flume Schützhold and Unruh (2002); see Barceló (2019); Barceló et al. (2005) for other examples and for the link with black hole physics. In 2008, the group led by Ulf Leonhardt showed the existence of a fiber-optical version characterizing the scattering of light on a soliton Philbin et al. (2008, 2007), see Belgiorno et al. (2010); Ciret et al. (2016); Ciret and Gorza (2016); Drori et al. (2019) for subsequent experimental works based on this proposal. In these settings, dispersive effects play a key role. Indeed, the stability of the soliton requires the dispersion relation of light propagating in the waveguide to display an anomalous behavior in the spectral region occupied by the soliton Agrawal (2007). Furthermore, the smallness of $\delta n$, the modification of the refractive index induced by the soliton (typically $\sim 10^{-4}$), implies that the scattering coefficients are correspondingly small Robertson (2011) and that one effectively works in a highly dispersive regime where the emitted spectrum is not governed by the “surface gravity” of an analogue horizon Finazzi and Parentani (2012). In contrast to flowing BEC where the calculation of scattering coefficients is rather straightforward as the spatial gradient near the sonic horizon is dominant Leonhardt et al. (2003a); Macher and Parentani (2009a), the smallness of $\delta n$ in the present settings complicates the calculation since it is not a priori clear which approximations are legitimate.

In this paper, we shall clarify this point by computing the scattering coefficients by two different means. As in Leonhardt et al. (2003a); Macher and Parentani (2009a), we shall first work in the background field approximation by solving the linear wave equation describing mode propagation on top of the background, here an optical soliton propagating in a nonlinear waveguide (WG). We shall then numerically solve the nonlinear Schrödinger equation (NLSE) governing the propagation of all the optical fields in the WG, i.e. both the soliton and small perturbations, along lines close to those adopted when numerically studying the analogue Hawking effect in a flowing BEC Carusotto et al. (2008).

In the first approach, the wave equation governing the mode mixing on the soliton is carefully derived by a linearization of the NLSE. In this we follow and generalize the treatment of four-wave mixing (FWM) presented in Yulin et al. (2004). As in the standard Bogoliubov treatment of linear perturbations in atomic BEC Macher and Parentani (2009a) or in polariton systems Marino (2008); Gerace and Carusotto (2012), the wave equation contains a term which mixes configurations of opposite detuned frequencies, and which would be dropped if performing the rotating wave approximation (RWA). With respect to atomic BEC, the novel element here arises from the non-monotonic character of the dispersion relation, which is here expressed as a relation between the detuned wave number (measured with respect to the carrier wave number of the soliton) and the detuned frequency. Crucially, when keeping the fourth-order term in the Taylor expansion of the dispersion relation, the positive- and negative-norm branches cross each other, a condition known as phase matching in optics, see e.g. Wadsworth et al. (2004); Rarity et al. (2005). In the presence of a continuous wave (CW) laser beam, this would induce a modulation instability in a narrow frequency window near the crossing, as was shown in Refs. Abdullaev et al. (1994); Pitois and Millot (2003). When replacing the CW by a soliton, we here show that the crossing leads to a significant increase of the pair production. As a result, we obtain a broad spectrum centered around the frequency of the crossing, and not at the frequency at which one might have expected the analogue Hawking signal to reach its maximum value, namely with one photon in the frequency domain where there is a strong conversion from probe to idler, or (in analogue terms) where there is a group velocity horizon. 111When comparing the spectrum to those found in previous studies of the analogue Hawking effect in waveguides Philbin et al. (2008); Robertson (2011); Finazzi and Carusotto (2013); Bermudez and Leonhardt (2016); Jacquet and König (2018); Drori et al. (2019), there is an overall enhancement factor of $\sim 10^{6}$. This large factor is due to the combination of two effects. The first is the crossing per se, as just described. The second requires a distiction to be made between “soft” and “hard” photon production processes. By “soft” we mean the processes that involve the standard FWM known in the nonlinear optics literature and which conserve the total number of photons (soliton plus perturbations), see our Eq. (5). Instead, “hard” processes stem from the second-order character of the original Maxwell equations, see e.g. Ref. Amiranashvili (2016). It is worth noticing that the above-mentioned works have focused on “hard” processes and ignored “soft” ones. It remains to identify the appropriate treatment in which both “soft” and “hard” processes are taken into account, thereby providing a complete description of the scatttering on the soliton. We are currently working on this.

To confirm this result, which is obtained from the linearized wave equation, we then perform numerical simulations describing the evolution of the full field (i.e. soliton $+$ perturbations), governed by the NLSE. We first study stimulated processes triggered by sending a probe wave near the phase matching condition. We recover the enhanced photon pair creation rate and establish that the total number of photons is conserved. As in a BEC, the produced quanta originate from particles extracted from the “condensate” (i.e. photons from the soliton). By replacing the probe wave by vacuum fluctuations, we then present (in a preliminary analysis) the whole spectrum and the correlations of spontaneously produced pairs. This is technically achieved by implementing the doublet formalism of Ref. Leonhardt et al. (2003b) at the level of the NLSE.

The paper is organized as follows. In Section II, we first present the nonlinear wave equation and its soliton solution. We then derive its linearized version describing waves scattered on the soliton, and briefly present the realistic example of a silicon nitride waveguide used in the numerical simulations. In Section III, we present the results displaying the above-mentioned enhanced anomalous coefficients by numerically solving for the eigenmodes of the linearized equation. We summarize and conclude in Section IV. In Appendix A, we consider the propagation of quasiparticle excitations on top of a continuous laser beam at the soliton frequency. Appendix B is concerned with numerical integration of the NLSE.

II Settings

In this section we provide the theoretical background behind the two wave equations used in the following sections. We begin with the nonlinear Schrödinger equation (NLSE) describing the evolution of the total field (background + fluctuations). The soliton propagation, which provides the background field, is a solution of this equation when omitting higher order terms in the dispersion relation. The NLSE is then linearized in order to describe fluctuations propagating on top of the soliton. Finally, the dispersion relation of the linearized modes and the scattering coefficients relating asymptotic modes are discussed.

II.1 The nonlinear wave equation

Let $z$ and $t$ be the standard space and time coordinates in the lab frame, in which the waveguide is at rest. Given our purposes, it is appropriate to write the (complex) electric field $E(t,z)$ as a slowly-varying envelope $A(t,z)$ multiplying a carrier continuous wave:

[TABLE]

where $\beta_{0}=\beta\left(\omega_{0}\right)$ is the wave number solution of the dispersion relation of linear waves in the waveguide. 222Conventionally in guided optics the dispersion relation between the wave number $k$ and the angular frequency $\omega$ is written $k=\beta(\omega)>0$. The absence of a square conveys the fact that only forward-propagating waves are considered. Note also that we have selected a single mode (of fixed polarization and with no nodes in the transverse direction) within the waveguide, and we do not consider processes which couple this mode to others.

In this work, we shall use the following effective equation for $A$ (see e.g. Ref. Agrawal (2007) for its derivation):

[TABLE]

Here, $\gamma>0$ is the nonlinear coefficient, while the operator $\beta(\omega_{0}+i\partial_{t})-\beta_{0}$ gives the detuned wave number $\beta-\beta_{0}$ as a function of the detuned frequency $i\partial_{t}=\Delta\omega=\omega-\omega_{0}$ of linear waves. In Eq. (2) several terms have been neglected, such as those encoding linear and nonlinear loss, retardation, etc.; see Refs. Abdullaev et al. (1994); Dudley et al. (2006); Ciret et al. (2016) for works in which some of these terms are discussed.

To prepare the analysis where the carrier wave is replaced by a soliton solution, it is appropriate to introduce a new coordinate system in which homogeneity (i.e. invariance under translations in $z$) will be preserved as an exact symmetry. To this end, we write the operator on the right-hand side as a Taylor series:

[TABLE]

where $\beta_{1}=\beta^{\prime}(\omega_{0})$, $\beta_{2}=\beta^{\prime\prime}(\omega_{0})$, etc., and then bring the term in $\beta_{1}$ over to the left-hand side. This gives

[TABLE]

Note that $\beta_{1}=dk/d\omega$ is the inverse group velocity of the carrier wave. One then introduces a new time coordinate $\tau\doteq t-\beta_{1}z$, so that $\partial_{z}|_{\tau}=\partial_{z}|_{t}+\beta_{1}\partial_{t}$. Eq. (4) thus reads

[TABLE]

where $B\left(\Delta\omega\right)$ is the detuned wave number of linear waves in the $(z,\tau)$ coordinate system as a function of the detuned frequency:

[TABLE]

As it is more appropriate in this work, we shall use the function $B\left(\Delta\omega\right)$ itself rather than its Taylor development. In this respect, however, we should point out that the term in Eq. (3) proportional to $\beta_{4}$ (which is often neglected) plays a key role in the following analysis; see Abdullaev et al. (1994); Pitois and Millot (2003) for works where the narrow resonance induced by $\beta_{4}$ when sending a CW (rather than a soliton) is discussed.

In the sequel, several conserved quantities will play important roles. When dealing with the nonlinear equation (5), for any complex solution $A$ the norm

[TABLE]

is $z$-independent and positive definite. In quantum settings, when properly normalized, this conserved quantity gives the total number of (forward-propagating) photons injected in the waveguide. This conservation law follows from Eq. (5) and prevents the inclusion of the “hard” processes mentioned in footnote 1. It should also be mentioned that there is another conserved quantity, namely

[TABLE]

This real quantity gives the total momentum carried by $A$ in the ($z,\tau$) coordinate system. In the present settings, it acts as a Hamiltonian, the right-hand side of Eq. (5) being equal to $\delta K_{A}/\delta A^{\star}(\tau,z)$ (see Eq. (5.7) in Pitaevskii and Stringari (2003) for the corresponding equation in an atomic BEC). Therefore, the coordinate $z$ parameterizes the evolution of the field configurations in the waveguide, thus acting as a temporal coordinate.

II.2 Silicon nitride waveguides

In this work we consider more specifically the propagation in a silicon nitride (SiN) waveguide. Following the continuous improvement of their fabrication processes, SiN WGs are nowadays one of the best nonlinear integrated platform avaiable thanks to a high nonlinear index, very low linear ($<$1dB/m) and nonlinear propagation losses and a broad transparency window from the visible to the mid-IR  Halir et al. (2012); Epping et al. (2015); Ji et al. (2017). The combine material and geometric dispersion allows for dispersion engineering that enables, together with low Raman scattering, the propagation of ultra-short fundamental solitons that do not suffer from soliton self frequency shift as in silica-based optical fibres Klenner et al. (2016). In Sec. B, we validate the neglect of linear losses in some dedicated numerical simulations based on the NLSE. Nonlinear losses will not be considered in this work.

We work with a realistic example of a rectangular silicon nitride waveguide (of height $0.5\,\mu{\rm m}$ and width $1.2\,\mu{\rm m}$) on a silica substrate. We took into account the profile of monochromatic probe waves in the perpendicular directions to numerically derive (using the software Lumerical) the effective one-dimensional dispersion relation in the longitudinal direction. This dispersion relation is shown in Figure 1. As $\beta(\omega)$ is approximately linear, we present the index $n(\omega)$ defined by $n(\omega)=c\beta(\omega)/\omega$, which varies more significantly. We also present the group index $n_{g}(\omega)$, which gives the group velocity through the relation $n_{g}(\omega)=c/v_{g}(\omega)$. Of key importance for the sequel is that the dispersion is anomalous in the vicinity of $\omega=2\,{\rm PHz}$, by which we mean the group velocity $v_{g}$ increases with $\omega$.

II.3 Solitons as

approximate background solutions

With the exception of plane waves (which are discussed in Appendix A), we are not aware of any stationary solution of Eq. (5) when the dispersion relation is given by that shown in Fig. 1. However, if we restrict the Taylor series of Eq. (3) to second order, i.e. if we take $B(\Delta\omega)=B_{2}(\Delta\omega)\doteq\beta_{2}(\Delta\omega)^{2}/2$, then there are well-known solutions which describe a one-parameter family of solitons Agrawal (2007). That is, when considering the nonlinear equation

[TABLE]

it is straightforward to show that this is solved by

[TABLE]

where

[TABLE]

These can be parameterized by the soliton duration $\tau_{0}$; indeed, the linearized wave equation will be seen to depend only on $\tau_{0}$ and $\delta\beta_{0}$, so that we shall have no need to separately consider $P_{0}$ or the nonlinear coefficient $\gamma$. Such solitons will be fair approximations to exact solutions (so long as their spectral width is not too large, or equivalently, their duration $\tau_{0}$ is not too small), and we shall treat a particular soliton as the background when considering the propagation of fluctuations in the next section. Hence it will play the same role as, say, the inhomogeneous flow over an obstacle in a water flume Lawrence (1987); Schützhold and Unruh (2002); Weinfurtner et al. (2011); Unruh (2013); Coutant and Parentani (2014); Michel and Parentani (2014); Euvé et al. (2016) when considering the scattering of surface waves.

It is here worth mentioning that the inexactness of the soliton as a solution of Eq. (5) will mostly manifest itself by a significant emission of Cherenkov radiation (CR) Akhmediev and Karlsson (1995), i.e. waves with the same detuned wave number as the soliton, in close analogy with the zero-frequency macroscopic undulation found in stationary flows over an obstacle Lawrence (1987); Coutant and Parentani (2014); Euvé et al. (2016). The CR will not enter the calculations of the scattering coefficients describing linear perturbations propagating on top of the soliton, because the soliton will be treated as an exact background solution when linearizing the wave equation. However, it will be explicitly seen when numerically solving the NLSE (see Appendix B).

As is well known, and as indicated by Eq. (11), $\beta_{2}<0$ is a necessary condition for the existence of stable soliton solutions, i.e. the dispersion at the soliton frequency must be anomalous. Since $\beta_{1}=\beta^{\prime}(\omega)$ is the inverse group velocity, this means that the soliton can only exist in a region where the group velocity is increasing with frequency (as shown in the right panel of Figure 1). Importantly, the finite size of the anomalous region implies a particular structure of the dispersion relation, which in turn necessarily leads to the narrow resonance mentioned in the introduction.

II.4

The linear wave equation

As in other treatments in media, we now linearize Eq. (5) around a background solution so as to get the equation governing the evolution of probe waves. In this section, to separate the emission of CR from the scattering of linear perturbations, we assume we have a homogeneous background field $e^{i\delta\!\beta_{0}z}A_{0}(\tau)$ that exactly solves Eq. (5). Perturbations on top of such a solution are conveniently described by the decomposition

[TABLE]

where the $z$-oscillation of the background governed by $\delta\beta_{0}$ has been factored out. As a result, the dispersion relation of asymptotic modes propagating on top of the soliton will explicitly depend on $\delta\beta_{0}$ (see the next subsection). The homogeneity of the background $e^{i\delta\!\beta_{0}z}A_{0}(\tau)$ ensures that the detuned wave number of linear probe waves will be conserved. Its value, which we shall call $D$, therefore plays the role of the conserved (Killing) frequency in the Hawking effect, both in the relativistic calculation Hawking (1975); Brout et al. (1995a) and in its analogue/dispersive version Brout et al. (1995b); Schützhold and Unruh (2002); Leonhardt et al. (2003a); Robertson (2011).

If the magnitude of $\delta A$ is small enough, it will obey the linearized form of Eq. (5):

[TABLE]

Performing now the RWA implies the dropping of the last term in $\delta A^{\star}$. This is equivalent to assuming no “phase matching” Agrawal (2007), i.e. that the oscillating part of $A_{0}^{2}\,\delta A^{\star}$ is off-shell, lying far from the dispersion relation. Consideration of existent works on the analogue Hawking effect in waveguides Philbin et al. (2008); Finazzi and Carusotto (2013); Bermudez and Leonhardt (2016); Jacquet and König (2018) indicates that this approximation was effectively applied, thereby missing a resonance that shall play an important role in the scattering of linear waves.

Keeping the term in $\delta A^{\star}$, the solutions of the wave equation are naturally doublets of the form $w\equiv\left[w_{+},\,w_{-}^{\star}\right]$. Details are given in Ref. Leonhardt et al. (2003b) in the context of BECs. (Note in particular that $w_{+}$ and $w_{-}$ play roles similar to those of $u_{k}$ and $v_{k}$ in the Bogoliubov-de Gennes formalism.) As in Leonhardt et al. (2003b), the doublet obeys the same equation as that obeyed by $\left[\delta A,\,\delta A^{\star}\right]$, i.e.

[TABLE]

where the differential operator $O$ is

[TABLE]

Note that, for this form of ${O}$ to make sense, $w$ should be understood as a column vector in Eq. (14). We define the conjugate $\bar{w}$ of $w$ as follows: if $w=\left[w_{+},\,w_{-}^{\star}\right]$, then $\bar{w}\equiv\left[w_{-},\,w_{+}^{\star}\right]$. It is easily shown that, if $w$ is a solution of Eq. (14), then so is $\bar{w}$ (i.e. the equation is invariant under the above conjugation). It is also worth noting that, had we not factored out $e^{i\delta\beta_{0}z}$ from $\delta A$ in Eq. (12), the off-diagonal elements of $O$ would be proportional to $e^{\pm 2i\delta\beta_{0}z}$. In this way, the $z$-independence of $O$ as defined above justifies the decomposition in Eq. (12).

Our aim now is to perform second quantization, so that annihilation and creation operators are associated with modes of positive and negative norm, respectively. To this end, we turn to the scalar product which is used to define the norm of the modes. Equation (14) possesses indeed a conserved scalar product: given two solutions $w_{1}=\left[w_{1,+}\,,\,w_{1,-}^{\star}\right]$ and $w_{2}=\left[w_{2,+}\,,\,w_{2,-}^{\star}\right]$, it is given by

[TABLE]

Using Eqs. (14) and (15), one verifies that

[TABLE]

The norm of a doublet $(w,\,w)$ is thus conserved 333It is to be noted that, when considering a more general $z$-dependent background solution, Eq. (17) is still satisfied, as is the Klein-Gordon norm in a time-dependent spacetime geometry., with the two components $w_{+}$ and $w_{-}^{\star}$ giving respectively the positive- and negative-norm content of $w$. 444If $w$ is of the form $\left[\delta A\,,\,\delta A^{\star}\right]$, then it necessarily has zero norm. This justifies considering the full space of complex solutions, so as to have well-defined positive- and negative-norm modes that are necessarily involved in the second quantization of the field. One should also note that the above norm should be distinguished from that of the nonlinear equation given in Eq. (7). Nevertheless, as we shall show below, these conserved norms are intimately related when considering asymptotic configurations.

The non-positive character of the above scalar product characterizes excitations of bosonic fields, and is also found when considering the solutions of the Bogoliubov-de Gennes equation in BEC physics or those of the Klein-Gordon equation in relativistic Quantum Field Theory. The mixing of modes of opposite norm is called anomalous scattering, and corresponds to the creation of pairs in quantized settings. When the initial state is vacuum, one then gets spontaneous pair creation, as on-shell particles (photons) are found in the output channel.

II.5 The two dispersion relations governing “soft” processes

Since the operator $O$ in Eq. (14) in independent of $z$, one can work at fixed detuned wave number $D$ and consider the globally defined doublets $e^{iDz}w_{D,\,j}(\tau)$, where $j$ labels the various linearly independent solutions. To compute the scattering coefficients encoded in these doublets, it is necessary to identify the complete set of asymptotic solutions defined far away from the soliton, i.e. for $\tau\rightarrow\pm\infty$.

In these asymptotic regions, $A_{0}(\tau)$ vanishes and Eq. (14) becomes independent of $\tau$. As a result, the two components of any doublet $w=(w_{+},w_{-}^{*})$ decouple, as is the case in atomic Bose condensates when interactions can be neglected after having opened the trap Tozzo and Dalfovo (2004); Robertson et al. (2017). Hence the asymptotic values of the detuning frequency $\Delta\omega$ are also conserved. Far away from the soliton, both components of doublets will thus be superpositions of plane waves of the form $e^{iDz}\,e^{-i\,\Delta\omega_{j}(D)\,\tau}$ where the various $\Delta\omega_{j}(D)$ will be algebraically related to $D$. Because positive- and negative-norm modes mix with each other when scattered, we must consider both sets of asymptotic modes, and thus two sets of roots: $\Delta\omega_{j,+}^{\rm as}(D)$ for the former and $\Delta\omega_{j,-}^{\rm as}(D)$ for the latter.

These roots are immediately obtained by considering separately the diagonal terms of the matrix $O$ of Eq. (14) with $A_{0}=0$. The upper (positive-norm) components $w_{D,\,j,\,+}(\tau)$ are characterized by detuned frequencies $\Delta\omega_{j,+}^{\rm as}(D)$ which satisfy

[TABLE]

Instead, lower (negative-norm) components $w^{\star}_{-D,\,j,\,-}(\tau)\propto{\rm exp}(-i\,\Delta\omega_{j,-}^{\rm as}(D)\,\tau)$ have their detuned frequencies $\Delta\omega_{j,-}^{\rm as}(D)$ obeying

[TABLE]

One easily verifies that $D_{-}^{\rm as}(\Delta\omega)=-D_{+}^{\rm as}(-\Delta\omega)$, which expresses the fact that negative-norm modes are conjugates of positive-norm modes. Hence the two dispersion relations are related by a $\pi$ rotation around $\left(\Delta\omega=0\,,\,D=0\right)$, or equivalently the point $(\omega=\omega_{0}\,,\,k=\beta_{0}+\delta\beta_{0})$ describing the soliton. Note that this is not the same as in Refs. Philbin et al. (2008); Robertson (2011); Finazzi and Carusotto (2013); Bermudez and Leonhardt (2016); Jacquet and König (2018); Drori et al. (2019), where positive- and negative-norm modes are related by a $\pi$ rotation around the “absolute” origin $(\omega=0\,,\,k=0)$. This key difference distinguishes what we referred to in footnote 1 as “soft” and “hard” photon production processes: “soft” processes (as studied here) involve the creation of collective excitations on top of the soliton (as for phonons in atomic BEC), while “hard” processes involve the spontaneous creation of additional photons. While both the total momentum $K_{A}$ of Eq. (8) and the total number of photons $N_{A}$ are conserved by the “soft” processes we consider, it remains to be seen what are the conserved quantities when including “hard” processes. However, we expect that any $z$-dependence of $K_{A}$ and $N_{A}$ due to the latter will be very small because of their strongly non-adiabatic character, which greatly suppresses the corresponding scattering amplitudes. (The interested reader is invited to consult Chapters 9 and 10 of Ref. Robertson (2011) for typical values of the pair production rates due to “hard” processes, which seem to be in agreement with the recent observations reported in Drori et al. (2019).)

Using again the example waveguide whose dispersion relation is shown in Fig. 1, and choosing for the soliton the carrier frequency indicated by the black dot in that figure with $\delta\beta_{0}=0.7\,{\rm mm}^{-1}$, $D_{+}^{\rm as}(\Delta\omega)$ ($D_{-}^{\rm as}(\Delta\omega)$) is shown in Fig. 2 by continuous and dotted (dashed and dot-dashed) curves. Importantly, the structure of $D_{+}^{\rm as}(\Delta\omega)$ possesses generic features. As we shall explain, these stem from the fact the soliton lives in a finite frequency window where the waveguide exhibits anomalous dispersion. To start the analysis, we first notice that at $\Delta\omega=0$, the positive-norm branch starts at $D_{+}(0)=-\delta\beta_{0}$, and because the soliton lives in a region of anomalous dispersion, this branch necessarily dips down to more negative values for small $\left|\Delta\omega\right|$. On leaving the window of anomalous dispersion, the curvature of the positive-norm branch then flips sign, rising to positive values of $D$ on both sides. As a result, it has a negative minimum for $\Delta\omega<0$ and another for $\Delta\omega>0$. Crucially, since the positive-norm branch is below the negative-norm branch at $\Delta\omega=0$ and above it for large $\Delta\omega$, it is necessarily the case that these two branches cross. In Fig. 2, the crossing occurs at $D_{\rm cross}=-28.4\,{\rm mm}^{-1}$ and $\Delta\omega_{\rm cross}=0.846\,{\rm PHz}$. (Of course there is another crossing for opposite values of $D$ and $\Delta\omega$, but this describes the same phase matching condition.) This crossing induces a modulation instability when sending a CW in the WG, as shown in App. A, and an enhancement of the anomalous scattering coefficient when sending a soliton, as we shall see below.

II.6 The various wave number bands

By examining Fig. 2, one also sees that for any $\Delta\omega$ there are two values for $D$, which correspond to modes of opposite norm. Instead, at fixed $D$ the number of roots $\Delta\omega_{j}$ varies between two and six, depending on the value of $D$. The full range of $D$ is thus split into a series of bands, in each of which the number of roots is fixed.

The domain of interest of these branches can be restricted to $D<0$, as the half plane with $D>0$ contains exactly the same modes (but with the signs of their norms flipped). The presence of the two minima and the particular behavior near $D=0$ due to the wave number shift $\delta\beta_{0}$ mean that there are four domains in total:

• •

For $D<-D_{\rm max,1}$, only the negative-norm branches exist. They are shown in blue dashed (for $\Delta\omega>0$) and yellow dot-dashed (for $\Delta\omega<0$) in Fig. 2, and the scattering matrix is an element of the unitary group $U(2)$. Since the two branches are well separated in frequency, there will not be any significant mixing between them.

• •

For $-D_{\rm max,1}<D<-D_{\rm max,2}$, there are four solutions of the dispersion relation: the two negative-norm plus two positive-norm branches shown in solid black and dotted red. The scattering matrix is now an element of the pseudo-unitary group $U(2,2)$. The merging of the two positive-norm branches for $D$ near $-D_{\rm max,1}$ means we can expect significant mixing between these modes. We can also expect significant mode mixing between the dashed blue and solid black branches near the value of $D_{\rm cross}$. Crucially, since these modes have opposite norm, we have here an enhanced anomalous scattering.

• •

For $-D_{\rm max,2}<D<-\delta\beta_{0}$, all six of the visible branches of the dispersion relation are present. The scattering matrix is an element of $U(4,2)$ since the two new modes have positive norm. Note that they merge when $D\to-D_{\rm max,2}$. We can thus expect significant mixing between them for $D$ slightly above $-D_{\rm max,2}$. We can also expect some significant mixing for $D\to-\delta\beta_{0}$ where the positive-norm branches shown in light dotted bronze and thick solid black merge.

• •

Finally, for $-\delta\beta_{0}<D<0$, there are four solutions of the dispersion relation: the negative-norm dot-dashed yellow and dashed blue branches, and the positive-norm solid green and dotted red branches. The scattering matrix in this range will thus belong to the group $U(2,2)$. This range is not particularly interesting for the scattering of linearized perturbations, because it is narrow in $D$ and because the scattering coefficients are small since the various roots are quite well separated in frequency. It is, however, the range in which Cherenkov radiation is produced; see Sec. B.

II.7 Scattering coefficients

We now turn to the scattering coefficients on the soliton. These describe the asymptotic mode content of the global solutions of Eq. (14), taking as background solution the soliton of Eq. (10). As already explained, the soliton is not an exact background solution, but it is very close to one. The exact solution will not be exactly stationary, in particular because of the emission of Cherenkov radiation. However, the deviations are very small. Taking the soliton as (approximate) background allows us to use the formalism of Section II.4, and obtain the scattering coefficients.

As mentioned below Eq. (15), if a doublet $e^{iDz}w_{D,\,j}(\tau)$ is a globally defined solution of Eq. (14), then its conjugate doublet $e^{-iDz}\bar{w}_{D,\,j}(\tau)$ is also a solution, although its norm has changed sign. Negative (positive) norm doublets will thus always be written with (without) a bar. This notation is in agreement with that used to describe a real scalar field in relativistic second quantized settings, namely negative, i.e. complex conjugated, (positive) norm modes are associated with creation (destruction) operators, and not treated as independent solutions Wald (1994).

The in- and out- doublets are respectively defined as having a single incoming mode (i.e. with a group velocity directed towards the soliton) in the asymptotic “past” (i.e $z\to-\infty$) and a single outgoing mode in the “future”. Because $A_{0}$ vanishes on both sides, the same set of roots characterizes both in- and out-modes. Therefore, when working at fixed $D$, the asymptotic postive- and negative-norm doublets read 555It should be noticed that the modes $w_{D,\,j}$ do not contain the standard relativistic normalization factor of $1/\sqrt{\omega}$. The origin of this is to be found in the quasi-monochromatic approximation for the slowly-varying envelope $A$ used to derive Eq. (2); see Ref. Brainis (2009) for a detailed discussion of this point.

[TABLE]

so that the set of in- or out-modes forms a complete and orthonormal basis with respect to the scalar product (16), i.e.

[TABLE]

The normalization imposed here explains the presence of the square root in Eq. (25) (as also in other media, see Macher and Parentani (2009b); Robertson (2011)). These bases are complete, so that the quantum field can be written as

[TABLE]

where $J_{+}(D)$ and $J_{-}(D)$ respectively represent the sets of available positive- and negative-norm solutions at the given value of $D$. The requirement that $w$ belong to the conjugation-invariant subspace is automatically satisfied since the creation operator $\hat{a}^{\dagger}$ is the hermitian conjugate of the annihilation operator $\hat{a}$ (as is the case for phononic field operators in atomic BEC). Namely, we adopt the standard bosonic commutation relation $\left[\hat{a}_{D,j}\,,\,\hat{a}_{D^{\prime},j^{\prime}}^{\dagger}\right]=\delta_{j,j^{\prime}}\,\delta\left(D-D^{\prime}\right)$ for both in- and out-modes.

Each in-mode, after scattering, becomes a linear superposition of out-modes. We thus have the transformation:

[TABLE]

for positive- and negative-norm incident modes, respectively. The coefficients $\alpha_{D,\,ij}$, $\alpha_{-D,\,ij}^{\star}$, $\beta^{\star}_{D,\,ij}$, and $\beta_{-D,\,ij}$ collectively form the scattering matrix restricted to a particular value of $D$. The normalization (26) implies the unitarity relation

[TABLE]

between the scattering coefficients of positive-norm doublets; an exactly analogous relation applies to $\alpha_{-D,\,ij}^{\star}$ and $\beta_{-D,\,ij}$ governing the scattering of negative-norm doublets. The $\alpha$ coefficients, which multiply out-modes of the same norm as the incident mode, are the standard amplitudes describing elastic scattering. On the other hand, the $\beta$ coefficients multiply out-modes of opposite norm to that of the incident mode and count negatively towards the unitarity relation.

Physically, $\left|\beta\right|^{2}$ gives the mean number (per unit wave number per unit length) of spontaneously created pairs of photons with opposite values of $D$. Returning to the lab frame, and considering the enhanced mode mixing between the solid black and dashed blue branches (with $D<0$ near $D_{\rm cross}=-28.4\,{\rm mm}^{-1}$), this means that the two wave numbers will be given by (see Eqs. (6) and (18))

[TABLE]

where for both solutions we have used the positive-norm branch solution $\Delta\omega^{\rm as}_{+}(D)$. When considering $D=D_{\rm cross}$, we find the simple relation $k_{\rm black}+k_{\rm blue}=2\left(\beta_{0}+\delta\beta_{0}\right)$, which is in agreement with the phase matching condition of Refs. Wadsworth et al. (2004); Rarity et al. (2005) and which tells us that each pair of quanta is extracted from a pair of “condensed” photons of the soliton. As far as we know, the kinematical relationships of Eqs. (30) have not been found in previous works on the analogue Hawking effect in nonlinear optics, since previous works focus on “hard” processes (following the distinction presented in footnote 1).

III Elastic and

anomalous scattering coefficients

In this section, we numerically calculate the scattering coefficients given the example dispersion relation and soliton carrier wave of Figs. 1 and 2. We focus on two particular processes which emphasize the roles played by the $\alpha$ and $\beta$ types of coefficients defined by Eq. (28). The plots have been obtained by numerically solving Eq. (14) at fixed $D$. The techniques used are the same as those described in Robertson and Leonhardt (2014), generalized to the case of a doublet. As in that paper, the wave equation is re-expressed as an integral equation in Fourier space, so that the dispersion term $B(i\partial_{\tau})$ becomes a multiplicative term and the Fourier transform of the squared soliton profile $\left|A_{0}(\tau)\right|^{2}$ becomes the kernel of an integral operator. Upon discretization, this integral equation becomes a matrix equation, which can be solved numerically after suitable regularization is performed. The validity of the results has been verified by checking their agreement with the unitarity relation (29). The interested reader is invited to consult Robertson and Leonhardt (2014) for a detailed description of the techniques employed.

III.1 Elastic scattering: From total transmission to

total reflection

First, let us consider what happens for $D$ near $-D_{\rm max,1}$, where the solid black and dotted red branches merge (see Fig. 2). When the probe interacts with the soliton, there is an effective deformation of the “local” dispersion relation due to the presence of the pulse. In effect, the dispersion relation is tilted such that the extremum at $-D_{\rm max,1}$ is reduced in magnitude (see the left panel of Fig. 8 in Appendix A). At the peak of the soliton, for $A_{0}^{\rm max}=\sqrt{P_{0}}$, this tilting reaches its maximum extent, and we refer to the deformed extremum at the soliton peak as $-D_{\rm min,1}$. 666The value of this extremum is found by computing the eigenvalue of the matrix ${O}$ of Eq. (15) when using the maximal value of $A_{0}$. A similar concept is found when considering non-monotonic subcritical flows over an obstacle in a flume. In that case, the flow velocity reaches a maximal value which determines an extremum $\omega_{\rm min}$ of the “local” dispersion relation on top of the obstacle Michel and Parentani (2014); Robertson et al. (2016). Frequencies below $\omega_{\rm min}$ are essentially transmitted across the obstacle, while those above $\omega_{\rm min}$ are essentially reflected; see Euvé et al. (2015, 2016) for experimental verifications in these settings. The similarity between the present settings and subcritical flows will also be found when considering the behaviour of the norm of some scattering coefficients, see footnote 7. Then, for $-D_{\rm max,1}<D<-D_{\rm min,1}$, there exist no real solutions of the “local” dispersion relation at the center of the pulse which are continuously connected to the probe frequency, and so the probe is (partially) blocked. It reflects from the pulse at a different frequency, whose asymptotic value is determined by the conservation of $D$: the solid black branch is shifted onto the dotted red. This is exactly the situation that was studied in Refs. Philbin et al. (2008); Ciret et al. (2016), where the frequency-shifting from the probe to the “idler” was observed, thereby revealing the presence of a “group velocity horizon” (see footnote 7 for more details, and also Ref. Choudhary and König (2012) for a discussion of the elastic scattering coefficients).

This situation is illustrated in the left panel of Fig. 3. The unit positive-norm incident mode (see Eq. (25)) lives on the dotted red branch, is characterized by $D$, and acts as a probe wave. When the magnitude of $D$ is far below $D_{\rm min,1}$ (indicated by the dashed vertical line in Fig. 3), it is essentially transmitted across the pulse. The transmission coefficient is equal to $1$, with all other scattering coefficients being negligible. However, there is a clear reversal around $D_{\rm min,1}$, with the transmission coefficient dropping to zero for $D$ significantly larger than $D_{\rm min,1}$, and the scattering coefficient describing red-to-black mixing (which is a reflection coefficient) climbing to $1$. In this regime, the incident wave is completely reflected onto the solid black branch. Note that the scattering is essentially a two-mode mixing process between modes of positive norm. The unitarity relation takes thus the form

[TABLE]

for all values of $D$ in the above-discussed intervals. The numerical analysis shows that the deviation from $1$ is bounded by $5\times 10^{-4}$. This is mostly due to there being a small amount of anomalous scattering not visible on the linear scale used here, but is clearly visible in the upper left panel of Fig. 4 below, where the same results are shown on a logarithmic scale.

III.2 Enhanced

anomalous scattering due to phase matching

Let us now turn to the right plot of Fig. 3, which shows the scattering coefficients when the incident wave is a unit norm mode living on the solid black branch of the dispersion relation in Fig. 2.

Firstly, we note that in the region around $D_{\rm min,1}$, the behavior of the scattering coefficients is essentially the same as described above for an incident mode on the dotted red branch, except that the solid black and dotted red curves have been switched. This is exactly as would be expected if the black and red modes are essentially decoupled from all other modes.

The situation is different, however, when $D$ is in the vicinity of $|D_{\rm cross}|\simeq 28{\rm mm}^{-1}$ where the solid black and dashed blue curves cross, see Fig. 2. There we find a significant scattering between these modes. Since they have opposite norm, there is an increase of the transmission coefficient above $1$, i.e. it is associated with amplification of the incident wave. This can be understood from the corresponding unitarity relation, which, since the coupling to other modes is negligible for $D$ near $D_{\rm cross}$, takes the form

[TABLE]

Here, $\alpha_{D}$ is the amplitude of the transmitted wave on the solid black branch, while $\beta_{D}$ is the amplitude of the outgoing wave on the dashed blue branch. It is clear that a non-zero $\beta_{D}$ requires $\alpha_{D}$ to have a magnitude larger than 1.

This anomalous two-mode mixing is responsible for the analogue Hawking effect. In the absence of enhancement, it reaches a maximum value near $D_{\rm min,1}$, see e.g. the dashed blue curve in the upper left panel of Fig. 4. This behavior was also found numerically in the wave blocking process of counter-propagating surface waves on a subcritical flow Michel and Parentani (2014); Robertson et al. (2016). The novelty here is that there exists a crossing point at $D_{\rm cross}$ where the relevant opposite-norm modes have exactly zero separation in Fourier space. This allows a great enhancement of the corresponding $\beta_{D}$-coefficient, to the extent that it becomes visible even on a linear (rather than logarithmic) scale. This is the main result of the present work.

III.3 Behavior of scattering coefficients on a log scale

Because most of the scattering coefficients are much smaller than $1$, it is appropriate to plot their squared magnitudes on a logarithmic scale. On the top row of Fig. 4, we represent the same scattering coefficients as in Fig. 3. On the left plot, we discover that there is a small anomalous coefficient with a maximum value $\simeq 5\times 10^{-4}$ (see the dashed blue curve). On the right plot, we observe that the anomalous coefficients extend over the entire range of $D$, with a minimum value of $\simeq 3\times 10^{-5}$. We also observe the subdominant scattering coefficients involving modes living on the solid green and dotted bronze branches (shown in a lighter line style).

On the bottom row, on the left plot, we represent the magnitudes of the scattering coefficients when sending a unit norm incident mode living on the dashed blue branch. The three other curves (namely solid black, dotted red, and light dotted bronze) describe the magnitudes of anomalous scattering coefficients involving modes with opposite norms. We notice that the solid black (dotted red) curve is similar to the dashed blue curve in the upper right (left) plot; a similar approximate symmetry was observed when studying a $U(2,1)$ $S$-matrix in atomic BEC (see Eq. (D8) in Macher and Parentani (2009a)) 777As already mentioned (see below Eq. (32)), the anomalous coefficient which is not enhanced by phase matching, namely the dashed blue curve in the upper left panel (or the dotted red curve in the bottom left panel), reaches its maximum value near $D_{\rm min,1}$, as was found for the anomalous coefficient governing the analogue Hawking effect in a subcritical flow over an obstacle; see the behavior of $\left|\beta_{\omega}\right|^{2}$ in the upper left plot of Fig. 5 in Robertson et al. (2016). This common behavior can be understood from the presence of a “group velocity horizon” (i.e., a turning point) Barceló et al. (2005) for quasiparticles, here photons with $D\in\left[D_{\rm min,1},D_{\rm max,1}\right]$ living on the black or red branch. Based on this common behavior, we believe that the anomalous “soft” processes we are studying should be considered as a new version of the analogue Hawking effect in nonlinear optics. . For completeness, on the lower right plot we represent the magnitudes of the scattering coefficients when sending in a mode living on the light dotted bronze curve. One mainly observes a linear mode mixing describing reflection and transmission between that branch and the corresponding modes on the light solid green one, as could have been expected from the behavior on the dispersion relation near $D_{{\rm max},2}$.

III.4 Spontaneously emitted spectra for various soliton durations

In quantum settings, when the initial state is vacuum, the squared magnitudes of anomalous scattering coefficients give the rates of spontaneous emission of photon pairs. More precisely, the number of pairs emitted per unit $D$ (the conserved wave number) per unit $z$ (the direction, conjugate to $D$, along which the wave equation is invariant) is Corley and Jacobson (1996)

[TABLE]

Integrating $\left|\beta_{D}\right|^{2}/2\pi$ over all $D$ thus gives the total emission rate (per unit length) of photon pairs. We now further consider the behavior of the enhanced anomalous coefficient governing the mixing of modes living on the solid black and dashed blue branches.

In Fig. 5, we represent on a linear scale $\left|\beta_{D}\right|^{2}$ for three different soliton durations: $\tau_{0}=10\,{\rm fs}$, $20\,{\rm fs}$ and $40\,{\rm fs}$. The emission spectrum is clearly sensitive to the value of $\tau_{0}$, and integrating over $D$ shows that the total emission rates are, respectively, $7.5\times 10^{-2}\,{\rm mm}^{-1}$, $8.7\times 10^{-3}\,{\rm mm}^{-1}$ and $1.1\times 10^{-3}\,{\rm mm}^{-1}$. That is, the total emission rate is roughly proportional to $\tau_{0}^{-3}$. This can be understood as follows. In Figure 8 in Appendix A is shown the imaginary part of the “instantaneous” dispersion relation at the peak of the soliton, which corresponds to unstable modes and is responsible for the strong anomalous mode mixing leading to the spontaneous emission here described. Figure 8 indicates that both the width and height of the unstable part of the spectrum vary roughly as $\tau_{0}^{-2}$. Its integral over $\Delta\omega$ thus varies as $\tau_{0}^{-4}$, and gives a photon production rate per unit $\tau$ per unit length. We may thus multiply by an effective interaction time to get the production rate per unit length. This interaction time is just the pulse duration $\tau_{0}$, leading to the total photon production rate being proportional to $\tau_{0}^{-3}$.

Since it is frequency (rather than wave number) that is measured at the output of the waveguide, it is appropriate to conclude this analysis by translating the spectra of Fig. 5 into numbers of photons emitted per unit $\Delta\omega$ per unit $z$. Because the group velocities of the relevant modes (namely, those on the solid black and dashed blue branches) significantly differ, the two emission powers are quite different. Their integrals over $\Delta\omega$, however, should agree, and they will be exactly the same as the integral of $\left|\beta_{D}\right|^{2}/2\pi$ over $D$. Indeed, it is this property – the invariance of the integral – that defines the emission spectrum in frequency space:

[TABLE]

The two spectra are shown in Fig. 6, for the same soliton durations used in Fig. 5. In the left panel is shown the emission spectrum on the (positive-norm) solid blue branch (i.e. that with $D>0$ and $\Delta\omega<0$ in Fig. 2), which is very narrow in frequency due to the large group velocity of this branch. The right panel shows the emission spectrum on the (positive-norm) solid black branch (with $D<0$ and $\Delta\omega>0$ in Fig. 2), which is by contrast wide in frequency and correspondingly smaller in amplitude. Interestingly, when using the shortest soliton duration $\tau_{0}=10\,{\rm fs}$, the values of the detuning at the maximal value of these two spectra are not exactly opposite to each other, as one might have expected from a naive use of the mode matching condition obtained when sending a CW in the WG. Indeed, for $\Delta\omega<0$ (along the solid blue branch) one finds that the maximum is located at $\Delta\omega=-0.84\,{\rm PHz}$, whereas for $\Delta\omega>0$ (along the solid black branch) it is localized at $\Delta\omega=0.78\,{\rm PHz}$. A very similar offset is found when integrating the nonlinear equation in App. B, see the two narrow peaks on the right lower plot of Fig. 9.

Finally, it is worth mentioning that the integrated power ($=7.5\times 10^{-2}\,{\rm mm}^{-1}$) obtained with a soliton duration of $\tau_{0}=10\,{\rm fs}$ means that each soliton propagating in a WG of 1cm will produce in the mean $0.75$ pair of entangled pairs with one photon living on the (positive-norm) solid blue branch (with $\Delta\omega<0$) while its partner lives on the (positive-norm) solid black branch (with $\Delta\omega>0$).

In Appendix B we numerically solve the nonlinear wave equation (5) forward in $z$ by sending both the soliton configuration of Eq. (10) and some small perturbation. In doing so, we include all effects encoded in the nonlinear equation, thereby generalizing what we just did by solving the linearized equation (14). The first outcome is that the above results are all recovered to a very good approximation. In addition, the new simulations display the consequences of the fact that the soliton we just used is not an exact solution of Eq. (5), since it exactly solves Eq. (9). Finally, since the new simulations deal with the full field, we are able to check the (nearly exact) conservation of the total number of photons, see $N_{A}$ of Eq. (7), which implies that each created pair of photons (either stimulated by sending a probe, or spontaneously produced from amplification of vacuum fluctuations) is accompanied by a corresponding decrease of the number of photons in the soliton.

IV Summary and conclusions

In this work, we studied the scattering of linear perturbations by a soliton propagating in a nonlinear optical waveguide (WG). We started by considering the wave equation governing the propagation of the full field (soliton + perturbations). As is usually done in these settings, the field is described by a slowly-varying envelope multiplying a given carrier wave. Hence the Fourier components of the envelope are characterized by the detuned frequency and wave number. This equation has been simplified in two respects. First, only configurations co-propagating with the soliton are described; hence it is odd in the wave number $k>0$. Given the smallness of the spatial gradients (fixed by the nonlinear Kerr index and the soliton power), this neglect is well justified. Second, we neglect the terms responsible for linear and nonlinear losses, retardation and self-steepening. None of these effects should significantly modify the main results obtained with the simplified equation. In fact, we believe our results are robust because they rely on a crossing of two branches of the dispersion relation, often referred to as a phase matching condition. This crossing induces a modulation instability when sending a continuous carrier wave in the WG (as is reviewed in App. A), and an enhancement of pair creation processes when sending a soliton (as shown in Sec. III).

For concreteness, the waveguide was taken to be a rectangular silicon nitride waveguide on a silica substrate. This type of WG offers two advantages, namely the possibility of engineering the effective longitudinal dispersion relation (which we compute numerically by taking into account the transverse properties of the WG), and a low loss rate which can be safely ignored. Importantly, the dispersion is anomalous in a finite frequency window. This is relevant in two respects. First, it allows soliton solutions which do not significantly disperse (as verified in App. B), and which play the role of the background when studying scattering processes. Second, when considering the positive- and negative-norm branches of the waveguide dispersion relation defined relative to the carrier of the soliton, the finite extension of the anomalous window necessarily gives the above-mentioned crossing of the two branches of the dispersion relation, as clearly illustrated in Fig. 2.

The anomalous mode mixing coefficients encoding production of photon pairs were first computed by linearizing the full equation on top of the soliton. Keeping all linear terms, we pointed out that one of them, which would have been dropped if we had performed the standard rotating wave approximation (RWA), mixes modes of opposite detuned frequencies and opposite norm. It is thus necessary to keep this the term. When doing so, we obtained a wave equation which has the same structure as that governing phonon scattering in an atomic Bose condensate. We thus applied the same techniques (based on the use of mode doublets) to compute the scattering coefficients. The numerical calculation of these coefficients was performed for a fixed value of the detuned wave number, since it is a conserved quantity given the stationarity of the background. We first recovered the expected elastic mode conversion interpolating from total transmission to total reflection which appears when the modes are blocked by the soliton. Such conversion, which occurs in a small interval of detuned wave numbers, has been described and observed in several works. Focusing on anomalous scattering coefficients which induce a parametric amplification (i.e. pair creation processes in quantum settings), we observed a significant enhancement of the coefficient in the vicinity of the crossing, see Fig. 6. As can be seen in the two plots of this figure, the detuned frequencies of the produced photons (which are spontaneously produced) are, to a good approximation, opposite to each other. Since their detuned wave numbers are also opposite, one clearly sees that the kinematics of the enhanced anomalous scattering on a soliton is governed by the phase matching condition. It also tells us that the pairs produced by the processes here considered are all extracted from pairs of “condensed” photons belonging to the soliton. When considering the pair production rate, we found it to be about 6 orders of magnitude larger than the results of previous works on the analogue Hawking effect in nonlinear optics. As explained in footnote 1, this discrepancy has a double origin: the crossing of the positive- and negative-norm branches of the dispersion relation, as well as the more adiabatic character of FWM with respect to the “hard” processes previously studied. When considering the scattering on a soliton, we are thus facing two versions of the analogue Hawking effect: the “soft” one here presented, and the “hard” one studied previously. It remains to provide a unified description where both types of anomalous scattering are simultaneously computed.

We then considered the integral over the detuned wave number so as to get the total pair production rate per unit propagation distance of the soliton. We found that it scales with the third power of the inverse duration of the soliton, or equivalently with $P_{0}^{3/2}$ where $P_{0}$ is the peak soliton power. For a pulse duration of $10\,{\rm fs}$ at 971 nm (1.94 PHz), we deduced that, in the mean, one pair will be spontaneously produced for a WG of a length equal to $1.3\,{\rm cm}$. The two photons of the pair are respectively generated around 687 nm ($\Delta\omega=0.8$ PHz) and 1720 nm ($\Delta\omega=-0.845$ PHz) where detectors are commercially available. This puts the effect within measurable range, and would enable the clear demonstration of photon pair production induced by a soliton. The fact that all photons are emitted in pairs suggests using coincidence measurements to greatly increase the signal-to-noise ratio: one would only analyze the data where two photons are detected, and verify that their frequency domains lie within the expected domains shown in Fig. 6. If the signal-to-noise ratio is high enough and the relevant photon modes are initially in their vacuum state, one could envisage measuring the degree of entanglement, as briefly discussed in the last subsection of Appendix B.

This appendix is mainly concerned with the numerical integration of the nonlinear equation governing the full field (soliton + perturbations). Although this treatment completely differs from that used when dealing with the linearized wave with a fixed detuned wave number, we find an excellent agreement of the results. In particular, the squared norms of the scattering coefficients extracted from nonlinear simulations closely agree in profile and magnitude with those obtained in the main text. This is non-trivial, because the evolution of the soliton is now dynamically determined. As a result, one observes a significant amount of Cherenkov radiation, which can be understood from the fact that the soliton is not an exact solution of the nonlinear equation. Morever, because we deal with the full field, we are also able to verify that pair production processes are always accompanied by a corresponding decrease of the number of photons in the soliton. Similarly, we also observe the reduction of the latter due to Cherenkov radiation, as expected from the fact that the first-order equation identically conserves the total number of photons, thereby encoding only “soft” processes in our classification.

Acknowledgments

We thank Florent Michel for useful discussions and for his careful reading of the manuscript. We are grateful to Ulf Leonhardt for useful remarks emphasizing the difference between “soft” and “hard” processes. S.R. thanks LPT (Laboratoire de Physique Théorique), Orsay, where most of the research work was done. S.R. and R.P. were supported by the French National Research Agency through the Grant No. ANR-15-CE30-0017-04 associated with the project HARALAB. S-P.G., C.C. and S.M. acknowledge the support of the Belgian Science Policy (BELSPO) Interuniversity Attraction Pole (IAP) 7-35 Photonics@be, as well as of the Fonds de la Recherche Fondamentale Collective (Grant No. PDR.T.1084.15).

Appendix A Dispersion relation on top of a strong background

In this appendix we consider the case where the background is homogeneous and thus described by a plane wave. Although not original, see Refs. Abdullaev et al. (1994); Pitois and Millot (2003), it is here presented in terms of the doublet formalism so as to clarify its link with the scattering described in the main text. Importantly, at the end of the appendix, we apply the same techniques to determine the effective deformation of the dispersion relation computed at the peak power of a soliton, a notion used to define the quantities $D_{\rm min,1}$ and $D_{\rm min,2}$ used in Figs. 3 and 4. To our knowledge the results presented in Fig. 8 have not previously been presented in the literature.

A.1 Homogeneous background

Starting from Eq. (5), we first consider the solution $A_{0}$ which is indepedent of $\tau$, and which thus satisfies the equation

[TABLE]

This is solved straightforwardly as follows:

[TABLE]

where $A_{0}^{H}$ and $\theta_{0}$ are constant real numbers, and where $\delta\beta_{0}^{H}=\gamma\left(A_{0}^{H}\right)^{2}$ is the nonlinear displacement of the wave number of the carrier wave. Notice that this differs by a factor of 2 from the relation between $\delta\beta_{0}$ and the peak power of a soliton; see Eq. (11). (In the analogy between Eq. (35) and the Gross-Pitaevskii equation, the nonlinear contribution $\delta\beta_{0}^{H}$ to the wave number of the $\tau$-independent background is analogous to the chemical potential of a homogeneous BEC.)

We now proceed as in Sec. II.4, writing the full solution as the sum of the background and a weak perturbation. As in Eq. (12), it is convenient to factor out the $z$-dependent phase of $A_{0}$:

[TABLE]

Then the phase of the perturbation $\delta A$ is relative to that of the background, whereas its amplitude is absolute.888We could, if we wished, also use the relative amplitude in the homogeneous case considered here, as is often done in treatments of phononic perturbations in atomic BEC Macher and Parentani (2009a). In the main text, however, where a localized pulse is considered as background, this would have been problematic as the “perturbation” amplitude would have been much larger than that of the pulse asymptotically. To linear order in $\delta{A}$ and $\delta{A}^{\star}$, the wave equation (5) becomes

[TABLE]

or, in matrix form,

[TABLE]

This is a particular realization of Eqs. (14) and (15), where here $A_{0}(\tau)=A_{0}^{H}$ is defined to be real and its magnitude is such that $\gamma\left(A_{0}^{H}\right)^{2}=\delta\beta_{0}^{H}$.

We may exploit the lack of any explicit $\tau$-dependence in Eq. (39) and perform the Fourier transform working at fixed (conserved) $\Delta\omega$:

[TABLE]

Then Eq. (39) becomes:

[TABLE]

where we have introduced the even and odd parts of $B(\Delta\omega)$, defined as

[TABLE]

This allows Eq. (41) to be written in the form

[TABLE]

where $K_{\Delta\omega}\equiv B_{\rm even}(\Delta\omega)+\delta\beta_{0}^{H}$. This equation is now very similar to that obtained for linear perturbations in a homogeneous BEC. It can be diagonalized by introducing the variables $\varphi_{\Delta\omega}$ and $\varphi^{\star}_{-\Delta\omega}$:

[TABLE]

where

[TABLE]

and the effective shift in the squared wave number

[TABLE]

Finally, this yields the diagonalized matrix equation

[TABLE]

where

[TABLE]

This, then, is the dispersion relation of linear perturbations on top of a constant background. Because $p_{\Delta\omega}^{2}$ of Eq. (46) is not necessarily positive, $D(\Delta\omega)$ is generally complex. Examples of its real and imaginary parts are shown in Fig. 7, for realistic values of $\delta\beta_{0}^{H}$ (namely, $2.5\,\mathrm{mm}^{-1}$ and $5\,\mathrm{mm}^{-1}$). In particular, in the imaginary part of $D(\Delta\omega)$ we recover the standard modulation instability around $\Delta\omega=0$, as well as the additional narrow instability induced by the crossing of the two branches of the (unperturbed) dispersion relation. It is worth reminding the reader that eigenmodes with a complex frequency necessarily have zero norm, see Leonhardt et al. (2003a); Coutant et al. (2016) and other references therein. This can be understood from the fact that they are formed from a linear combination of positive- and negative-norm modes with equal weights, causing the norm of the eigenmode to vanish.

It is instructive to consider limiting cases of the dispersion relation (48). When $\gamma\left(A_{0}^{H}\right)^{2}=0$, it reduces to $D(\Delta\omega)=B(\Delta\omega)$, as expected; $D(\Delta\omega)$ and $-D(-\Delta\omega)$ are just the detuned wave numbers of $\widetilde{\delta A}$ and $\widetilde{\delta A}^{\star}$, respectively. When $\delta\beta_{0}^{H}\neq 0$, we can consider the small $\Delta\omega$ behaviour using the fact that $B(\Delta\omega)\approx\frac{1}{2}\beta_{2}(\Delta\omega)^{2}$ at small $\Delta\omega$. To lowest order, then, $B_{\rm odd}(\Delta\omega)$ vanishes, and we have

[TABLE]

If $\beta_{2}<0$, then modes with small $\Delta\omega$ have an imaginary $D$, and so either grow or decay exponentially with $z$. On the other hand, if $\beta_{2}>0$, these modes propagate with a constant velocity in the co-moving frame of the carrier wave.

A.2 Local description on top of a soliton

We now use the above analysis to characterize the effective “instantaneous” (i.e. $\tau$-dependent) dispersion relation on top of the soliton. As in the standard WKB treatment, we neglect the gradient of the background and keep only the local value of the intensity $\left|A_{0}(\tau)\right|^{2}$. The difference from the above analysis lies in the fact that, while the wave number $\delta\beta_{0}$ of the background is a constant (now given by Eq. (11)), the intensity $\left|A_{0}(\tau)\right|^{2}$ is a function of $\tau$, and the two quantities are therefore disconnected. Moreover, as noted above, even when restricting our attention to the point at the peak of the soliton, the relationship between $\delta\beta_{0}$ and $\left|A_{0}^{\rm peak}\right|^{2}$ differs from that for a $\tau$-independent background by a factor of 2. As a result, the dispersion relations evaluated at the peak intensity of the soliton slightly differ from those of Fig. 7, both in their real and imaginary parts.

The relevant matrix to be diagonalized is that of Eq. (15) with $i\partial_{\tau}\rightarrow\Delta\omega$, and (as noted above) with $\left|A_{0}(\tau)\right|^{2}$ being treated as constant in order to extract the “instantaneous” dispersion relation. Following exactly the same procedure as for the $\tau$-independent case outlined above, we find again Eq. (48) where now

[TABLE]

At the peak of the soliton, we have $\gamma\left|A_{0}(\tau)\right|^{2}=2\,\delta\beta_{0}$, and this becomes

[TABLE]

Examples of the dispersion relation at the peak of the soliton are shown in Fig. 8 for the same soliton durations used in the main text. The corresponding extrema of the real part of the function $D(\Delta\omega)$ of Eq. (48) define the quantities called $D_{\rm min,1}$ and $D_{\rm min,2}$, which are represented by vertical dotted lines in Figs. 3 and 4. Moreover, the variation of the width and height of the imaginary part of $D(\Delta\omega)$ in the vicinity of $D_{\rm cross}$ gives a good estimate for the $\tau_{0}$-dependence of the total photon emission rate, as discussed in Sec. III.4.

Appendix B

Nonlinear propagation of the full field configuration

In this appendix, we numerically solve the nonlinear wave equation (5) forward in $z$ by sending both the soliton configuration of Eq. (10) and some small perturbation. In doing so, our aim is to include all effects encoded in the nonlinear equation, thereby generalizing and validating what has been done in the main text by solving the linearized equation (14). At a deeper level, we also aim at establishing the conservations laws associated with the nonlinear evolution, which have no counterpart when dealing with Eq. (14). We shall first study stimulated and then spontaneous processes. The various observations are presented in separate subsections so that the reader can easily identify both the agreement with, and the novel elements with respect to, the linearized treatment of the main text.

B.1 Stimulated emission by a probe wave

We start by considering the stimulated case because of the clarity of its outcome. As initial condition, we take the soliton $A_{0}(\tau)$ described in Eq. (10) with duration $\tau_{0}=10\,{\rm fs}$, plus an incoming probe wave $\delta A(\tau)$ living on the solid black branch (with $D<0$ and $\Delta\omega>0$) in Fig. 2, centered on the detuning frequency $\Delta\omega_{\rm probe}=0.78\,{\rm PHz}$ which is chosen to coincide with the maximum of the predicted emission rate (see the right panel of Fig. 6). The relative amplitude of the probe wave is such that the peak power ratio $P_{\rm probe}/P_{\rm soliton}\simeq 7\times 10^{-5}$, and it is described by a narrow Gaussian envelope of width $0.005\,{\rm PHz}$ in $\omega$-space.

B.1.1 Space-time description

To clearly separate the evolution of the field configurations describing the soliton from those describing the probe, the simulation is performed twice using as initial configurations $A^{\rm in}_{\pm}(\tau)=A_{0}(\tau)\pm\delta A(\tau)$. Since the initial soliton profile is exactly the same while the probe wave flips sign, this allows us to extract the part of the field which is linear in the probe wave by taking the half difference of the two runs, while the part of the field which is independent of the probe wave is found by taking their half sum. (Terms of quadratic order and higher in the probe field are negligible, as we have verified by using $A_{0}(\tau)$ as initial condition.)

The difference is represented in light gray in Fig. 9, while the sum is shown in thick black. Both have been normalized such that the peak value of the initial soliton field is $1$, both in real (left column) and in Fourier space (right column). On the left side, we clearly see the emergence of two wavepackets after the interaction with the soliton: one is short in duration with a large amplitude, and represents the transmitted probe wave, while the other is long in duration with a small amplitude, and represents the stimulated negative-norm wave on the solid blue branch (with $D>0$ and $\Delta\omega<0$) of Fig. 2. This can be verified by considering the right plots which represent, on a logarithmic scale, the squared norm of the Fourier component $\tilde{A}(\Delta\omega)$ as a function of the detuning frequency. It is interesting to notice that the small shift (to the left) of the negative detuning frequency of the stimulated emission on the blue branch is here recovered: namely, as was observed in Fig. 5, it is maximal for $\Delta\omega\simeq-0.84\,{\rm PHz}$ while that of the incoming probe is $\Delta\omega_{\rm probe}=0.78\,{\rm PHz}$.

In the plots of the right column of Fig. 9, we also clearly see the production of the “idler” (reflected wave), albeit with a much smaller amplitude, in good agreement with the reflection coefficient (shown in dotted red in the upper right plot of Fig. 4) for $D$ near $|D_{\rm cross}|=28\,{\rm mm}^{-1}$.

B.1.2 Cherenkov radiation

Focusing on the thick black curves which describe the background configurations, we observe, in both columns of Fig. 9, the emission of Cherenkov radiation Akhmediev and Karlsson (1995). The CR is emitted in both directions, as can be understood by noting the opposite slopes at the two $D=0$ roots of the dispersion relation (see Fig. 2). However, since the smaller frequency difference occurs for the root which has a negative group velocity, the CR emission rate should be greater on the left side, and this is borne out by Fig. 9. The larger CR emission occurs at $\Delta\omega\simeq-0.75\,{\rm PHz}$, in good agreement with the $D=0$ intercept. This emission is due to the fact that the soliton configuration injected in the WG is not an exact solution of Eq. (5). However, since the soliton is an exact solution of Eq. (9), the modification $\delta A_{0}(\tau,z)$ of the background soliton solution obeys a “forced” equation with the following structure (see Eq. (13)):

[TABLE]

where $\Delta B(\Delta\omega)\doteq B(\Delta\omega)-B_{2}(\Delta\omega)$ is the modification of the wave number with respect to the quadratic approximation of the dispersion relation mentioned above Eq. (9). (In obtaining this equation, we have also applied the rotating wave approximation, which amounts to neglecting the last term of Eq. (13).) As a result, the CR so produced is described by a coherent (displaced) state. 999This forcing should be clearly distinguished from “soft” squeezing due to FWM discussed in Tran et al. (2011), as well as from its “hard” version discussed in Ref. Rubino et al. (2012). As was noticed in Efimov et al. (2005), the CR always has the same polarization as the soliton, indicating that the dominant contribution is from the forced channel. It is also worth mentioning that forced and squeezed channels exist when considering undulations downstream from an obstacle in a flume Shen (1993), and in transonic flows in atomic BEC Busch et al. (2014).

B.1.3 Local decomposition of the probe field

Since the nonlinearities governed by the last term in Eq. (5) are small, typically $\gamma|A_{0}^{\rm peak}|^{2}/|D_{\rm cross}|\sim 0.05$ for a soliton with $\tau_{0}=10\,{\rm fs}$, it is meaningful to decompose the probe field $\delta A(\tau,z)$ at each $z$ into its “instantaneous” Fourier components $\tilde{\delta A(\Delta\omega,z)}$. One can then integrate the squared magnitude of $\tilde{\delta A(\Delta\omega,z)}$ over the relevant detuned frequency windows (defined by the asymptotic dispersion relation of Fig. 2) to obtain the evolution in $z$ of the transmitted, reflected and stimulated parts of the probe. More precisely, the instantaneous squared norm of modes living on the solid black branch (with $\Delta\omega>0$ and $D<0$) of Fig. 2, called $\left|\alpha(z)\right|^{2}$, is found by integrating $\left|\widetilde{\delta A}\left(\Delta\omega\right)\right|^{2}$ over the interval $\Delta\omega\in\left[0\,,\,1.2\right]\,{\rm PHz}$, and dividing by the initial (i.e. $z=0$) value of this integral. Similarly the reflected $\left|R(z)\right|^{2}$, and stimulated $\left|\beta(z)\right|^{2}$, contributions are obtained by integrating over the relevant domains in $\Delta\omega\in\left[0\,,\,1.2\right]\,{\rm PHz}$, and dividing the results by the initial value of this integral of $\left|\alpha(z)\right|^{2}$. By construction, when evaluated for sufficiently large $z$ so that the scattering has taken place, these three functions give the squared norm of the asymptotic scattering coefficients (integrated over the narrow width of the probe field).

These three functions are shown by solid lines in the left column of Figure 10. One clearly sees that the mode mixing starts for $z\sim 10$mm and essentially ends for $z\sim 25\,$mm, as one might have expected from the space-time evolution of the probe field. To quantitatively verify that this decomposition is meaningful, we represent on the lower left panel of Fig. 10 the quantity $1-(\left|\alpha(z)\right|^{2}+\left|R(z)\right|^{2}-\left|\beta(z)\right|^{2})$ which contains the quadratic combination entering the unitarity relation, see Eq. (29). We note that the maximal value of this deviation $\sim 10^{-8}$ is reached for $z\simeq 19$mm, in the “middle” of the scattering, and is much smaller than both $\left|R(z)\right|^{2}$ and $\left|\beta(z)\right|^{2}$. Moreover, the final value reaches a constant $\sim 2\times 10^{-12}$ which is much smaller than any of the norms presented. 101010Though the coefficient which describes scattering onto the dot-dashed yellow branch of the dispersion relation (that with $\Delta\omega<0$ and $D<0$ in Fig. 2) has been neglected here, this will not account for the final error because the results of the linear treatment for a soliton duration $\tau_{0}=10\,{\rm fs}$ indicate that its squared magnitude reaches a maximum on the order of $10^{-17}$, five orders of magnitude smaller than the final error in the bottom left panel of Fig. 10 and fifteen orders of magnitude smaller than the maximum of the squared amplitude describing stimulation of the blue branch. Concomitantly, we verified that the asymptotic values of each contribution agree quite well with the theoretical predictions of Figs. 3 and 4 at the given value of $D=|D_{\rm cross}|$.

B.1.4 Conservation of photon number and small dissipation

Because we now deal with the total field $A_{0}+\delta A$, we can verify that the creation of photon pairs induced by the probe field is accompanied by a corresponding decrease of the residual number of photons in the soliton. To this end, in the right column of Figure 10 are shown plots related to the norm $N_{A}$ of Eq. (7), whose constancy in $z$ (when applied to the total field) implies conservation of the total number of photons. As for the instantaneous norm content of the perturbation field, the integration is performed in $\Delta\omega$-space, which is split into a “soliton” region ($-0.6\,{\rm PHz}<\Delta\omega<0.6\,{\rm PHz}$) and a “probe” region ($\Delta\omega>0.6\,{\rm PHz}$); the region $\Delta\omega<-0.6\,{\rm PHz}$ is not separately presented. In the upper and middle panels are shown the changes in photon number in these two regions, relative to the initial number of photons in the probe; these are relevant quantities because the scattering is (to a very good approximation) linear in the probe field.

In the upper plot, the solid curve shows the increase in the number of photons in the transmitted part of the probe field. One verifies that the increase (as a function of $z$) of the quantity $N^{\rm probe}(z)/N_{\rm in}^{\rm probe}-1$ closely agrees with $\left|\alpha(z)\right|^{2}-1$ which is represented in the upper left plot.

In the middle plot, the dashed curve shows the total decrease of the number of photons in the soliton, which includes a steady decrease due to the CR emission. To correct for this, we subtract $N^{\rm sol}(z)$ extracted from a simulation where the initial profile of $A_{0}$ is exactly the same while $\delta A$ is set to zero. The result is shown by the solid curve. One observes that the remaining decrease exactly compensates for twice the increase of $N^{\rm probe}$, which is as expected since each photon added to the transmitted probe field is accompanied by a partner in the negative-norm wave. Indeed, we have checked that the corresponding curve for the region $\Delta\omega<-0.6\,{\rm PHz}$, after having subtracted the CR contribution in the same manner as for the soliton, is exactly the same as that showing the increase of $N_{A}$ for the probe (i.e. the solid curve on the upper panel).

The lower right panel of Fig. 10 shows the relative change of $N_{A}(z)$ which gives the total number of photons (soliton $+$ probe). Its final value is of order $10^{-12}$, as was the final value of the deviation shown in the lower left column. We believe that these two residual deviations are due to accumulation of numerical errors. To validate this conjecture, we have increased the numerical step size $\Delta z$ by a factor of $2$ (i.e. from $0.8\times 10^{-3}\,{\rm mm}$ to $1.6\times 10^{-3}\,{\rm mm}$). We observed that the change in the total number of photons is reduced by nearly a factor of $2$, while the deviation from unitarity shown in the lower left panel of Fig. 10) increases by a factor of about $20$.

To complete this study, we have included an overall linear loss in Eq. (5). The resulting quantities are shown by dotted lines in all panels of Fig. 10. The loss rate $\Gamma$ was chosen such that $10\,\%$ of the total energy is lost by the time the soliton has passed through $40\,{\rm mm}$ of the waveguide 111111This corresponds to a loss rate of $0.114\,{\rm dB/cm}$, a value at least 20 times higher than that reported for state-of-the-art silicon nitride waveguides Ji et al. (2017).. In the plots, the loss rate has been corrected for by multiplying all curves by the exponential factor $e^{\Gamma z}$. The residual decrease seems to be due to the fact that, on account of having lost some energy by the time of the interaction between the soliton and the probe, the soliton has slightly decreased in power (increased in duration), so that the scattering coefficients are smaller in magnitude. This interpretation is confirmed by additional simulations in which the initial soliton power is chosen so that, at the time of the interaction, the soliton width is equal to its value in the lossless simulations. In this case, the “instantaneous” scattering coefficients become almost indistinguishable from their values in the lossless case. This establishes the robustness of our results against the introduction of a small linear loss.

B.2 Spontaneous scattering: spectrum and entanglement

We now turn to the numerical simulation of spontaneous emission from the soliton, which is the only physical effect when the modes (besides those describing the soliton) are initially in their vacuum state (leaving aside the emission of CR). We aim at obtaining both the power spectrum characterizing the emitted pairs, as well as the correlations between the two photons in each pair. Since we are effectively dealing with a $U(1,1)$ mode mixing when considering the significant pair production for $D$ near $|D_{\rm cross}|$, the strength of these correlations should be large enough to guarantee that the bipartite state of the emitted pairs is nonseparable Busch and Parentani (2014).

As a preliminary attempt 121212The analysis presented here is work in progress. We are grateful to Florent Michel for recent discussions about this method of handling vacuum fluctuations in nonlinear systems. to characterize the strength of the correlations, we have found it useful to adopt the doublet formalism at the level of the nonlinear equation (5). That is, guided by the structure of Eq. (14), we consider the following pair of nonlinear coupled equations:

[TABLE]

This system reduces to Eq. (5) whenever $A_{+}=A_{-}=A_{0}$. In addition, when writing the doublet $A=[A_{+},A_{-}^{*}]$ as a common background term plus perturbations on top of that background, i.e., $A=[A_{0}+\delta A_{+},A_{0}^{*}+\delta A_{-}^{*}]$, then to linear order in the perturbations the system reduces to Eq. (14). Hence the system here considered correctly describes both the nonlinear evolution (in $z$) of the background configurations $A_{0}(z,\tau)$ and the linear evolution of the perturbations described by the doublet $w=[\delta A_{+},\delta A_{-}^{*}]$.

The expressions on the right-hand side of Eq. (53) are the Euler equations of the following Lagrange density:

[TABLE]

where the new generator of $z$-translations is the real part of (cf. Eq. (8))

[TABLE]

Therefore, for any pair of coupled solutions $(A_{+},A_{-}^{\star})$, the generalized version of the scalar product of Eq. (7), i.e., ${\rm Re}\left\{\int\,{\rm d}\tau\,A_{-}^{*}A_{+}\right\}$, is identically conserved, as follows from the $U(1)$ invariance of $L$ under $A_{\pm}\to A_{\pm}e^{i\phi}$.

To introduce the notion of vacuum fluctuations, we use the following asymmetric initial conditions:

[TABLE]

The field $A_{0}(\tau)$ represents the initial configuration of the coherent background state, which is taken to be the soliton solution of Eq. (10) in the forthcoming simulations (exactly as in the top row of Fig. 9). The initial value of the perturbation is described by the doublet $\bar{w}_{\rm in}=[0,\delta A_{-}^{*}(\tau)]$, which represents a configuration describing vacuum fluctuations. By this we mean the following. The upper component of the doublet identically vanishes, as is the case in the second line of Eq. (25). The lower component is given by 131313In the simulations, we apply effective limits to the integral to avoid exciting unphysical branches of the dispersion relation, which appear when large values of $D$ are folded back into the first Brillouin zone defined by the step size $\Delta z$. We set to zero all modes with $\left|D\right|\gtrsim 40\,{\rm mm}$, corresponding to the integral range $\Delta\omega\sim\left[-1\,{\rm PHz}\,,\,2\,{\rm PHz}\right]$. On top of this, after taking the integral of Eq. (57) to get the initial fluctuations in $\tau$-space, we also apply a smooth (double tanh-shaped) filter so as to remove $\delta A_{-}^{\star}$ in the range $\tau\sim\left[-0.2\,{\rm ps}\,,\,0.2\,{\rm ps}\right]$. This allows us to study only the forward evolution of vacuum fluctuations that are not initially on the soliton; without the filter, we observe a significant transient expulsion of fluctuations from the soliton.

[TABLE]

In this expression, each amplitude $a^{*}_{\Delta\omega}$ is a random complex number governed by a Gaussian ensemble with a vanishing mean value and a variance which is independent of $\Delta\omega$. These therefore encode the creation of one photon in each Fourier mode. 141414Planck’s constant $\hbar$ provides an absolute normalization for the field fluctuations, and an absolute meaning as to the number of photons in any Fourier mode. However, if effects which are nonlinear in the fluctuation $\delta A_{-}$ are negligible, we may scale it up in magnitude for computational convenience, so long as the photon number is calculated relative to its initial value, this being always set to $1$. The contents of this field are defined only probabilistically, and in principle the extraction of mean values requires averaging over many realizations of $\delta A_{-}^{*}$, as is done when using the truncated Wigner approximation Sinatra et al. (2002); Robertson et al. (2018), (TWA).

The advantage of this approach is that, as for the linearized wave equation, anomalous scattering occurs between the two components of the doublet, while normal scattering does not. Therefore, only the spontaneously produced photons will appear in the upper component $A_{+}$, making their emission rates easy to extract. In fact, comparing this treatment to the TWA, the extraction of small numbers of photons is here much simpler, because in the TWA one deals with the expectation value of a symmetrized product of the field operator (that is, an anti-commutator). As a result, to get the emission rate from vacuum, one should remove the contribution (of the commutator) which is present even in the absence of pair creation. When the emission rate is very small, this method requires many realizations, as it involves extracting a small difference between two terms dominated by a contribution which is independent of the number of pairs produced.

B.2.1 Emission spectra

In practice we proceed by first evaluating the Fourier components of $A_{+}(\tau)$ at different values of late $z$. From these, we extract an average creation rate per unit $z$ per unit $\Delta\omega$ that we compare with the predictions of the linearized wave equation already shown in Fig. 6. It should be noted that, for any single discrete value of $\Delta\omega$, the Fourier components oscillate rapidly when varying $\Delta\omega$ since we deal with a single realization of a Gaussian ensemble. These oscillations can be suppressed either by considering an ensemble of initial realizations of $\delta A_{-}$, as one does when using the truncated Wigner approximation Sinatra et al. (2002); Robertson et al. (2018), or, since we use a high resolution in $\Delta\omega$, by binning over many Fourier components. For the sake of efficiency, we adopt the latter approach. The results are shown in Figure 11. Even after binning, the numerical results show quite large oscillations when using a linear scale, but the convergence towards the predictions of the linearized equation becomes evident when plotted on a logarithmic scale. Hence the above method works efficiently in the sense that all Fourier components (and hence all values of the detuning parameter $D$) are dealt with at once.

In Figure 12, we represent the squared magnitude of the positive-norm component $A_{+}(\Delta\omega,D)$, the soliton having been removed using a filter which vanishes for $\tau$ near zero. We have applied a Fourier transform in both $\tau$ and $z$ so as to get the power spectrum of the spontaneously created pairs in the dispersion relation plane, see Fig. 2. To avoid spurious signals which are very wide in $D$, we have applied a regularizing function which vanishes at $z=0$ and $z=40\,{\rm mm}$, thereby removing a strong discontinuity when the data is treated (as we do) as periodic in $z$. One observes that the field configurations here considered only live along the positive-norm branch. This is precisely the sought-for result from having put to zero the initial value of $\delta A_{+}$ when using a doublet $(A_{+},\,A_{-})$ to describe and propagate configurations describing the soliton plus vacuum fluctuations. It is clear that the strongest signals appear on the solid black ($\Delta\omega>0$, $D<0$) and solid blue ($\Delta\omega<0$, $D>0$) branches of the dispersion relation, with a faint signal on the dotted red (idler) branch that decreases with increasing $\Delta\omega$.

B.2.2 Correlations in emitted pairs

Finally, we study the correlations between the photons emitted in pairs. To this end, we consider the connected part of the first-order correlation function

[TABLE]

The use of the two-component formalism means $g_{1}$ is not necessarily symmetric when dealing with a single random configuration to describe vacuum fluctuations. Hence, we define its symmetrized version using a geometric mean:

[TABLE]

To obtain the correlation map shown in Figure 13, we have proceeded in three steps.

• •

First, we take an average over late values of $z$ (over the final $100$ steps of the simulation), which sufficiently suppresses products of modes which do not have equal and opposite $D$.

• •

We then bin neighboring Fourier modes together, suppressing some of the randomness due to the particular realization for the initial value of $\delta A_{-}$.

• •

Lastly, although only one set of random amplitudes for $\delta A_{-}$ is generated, eight different runs are performed: the sign of all fluctuations are flipped; the sign of fluctuations living only on the solid black branch ($\Delta\omega>0$, $D<0$) are flipped; and the sign of fluctuations living only on the solid blue branch ($\Delta\omega<0$, $D>0$) are flipped.

This combination of initial conditions efficiently suppresses spurious correlations which appear for a single realization, but which should not survive in a proper ensemble average taken over many realizations.

The norm of $g_{1}^{\rm symm.}(\Delta\omega\,,\,\Delta\omega^{\prime})$ so calculated is shown using a linear scale on the left panel in Figure 13. On the right panel, we show the complete set of correlations between positive-norm modes with detuning wave numbers $\Delta\omega_{j}^{+}$ which have equal and opposite values of $D$. By direct comparison one finds, as expected, that the strongest correlations occur between the black and blue branches of the dispersion relation. As expected as well, we also see a smaller correlation between the dotted red (idler) and solid blue branches (which also have opposite values of $D$).