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

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 analog 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\mathrm{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 analog Hawking effect where the pair creation rates were about six orders of magnitude smaller.

Figure 1

Computed dispersion relation (of one particular mode, see footnote 2) in a 500 nm $\times$ 1200 nm silicon nitride waveguide on a silica substrate (see text for more details). On the left is plotted $n\left(\omega \right)=c\beta \left(\omega \right)/\omega$, the index as a function of the frequency; on the right is plotted the group index $n_g\left(\omega \right)=c{\beta }^{\prime }\left(\omega \right)$, again as a function of frequency. The black dot shows the position of the pulse engendering the soliton (at $\lambda =970\mathrm{nm}$ or ${\omega }_0=1.94\mathrm{PHz}$), while the open circle shows where the third derivative ${\beta }_3\left(\omega \right)=0$ (at $\omega =1.86\mathrm{PHz}$). We have adopted this frequency range because, as we shall later see, there will only be significant pair production for $\mathrm{\Delta }\omega =\omega -{\omega }_0\in \left[-1,1\right]$ PHz.

Figure 2

Dispersion relation in the silicon nitride waveguide, using the detuned variables $D$ and $\mathrm{\Delta }\omega$ giving, respectively, the difference in wave number and frequency of the modes with respect to those of the soliton (equal to ${\beta }_0+\delta {\beta }_0$ and ${\omega }_0$), whose carrier (at 1.94 PHz) is indicated by the black dot in Fig. 1. The solid and dotted curves correspond to the positive-norm modes of Eq. (18), while the dashed and dotted-dashed curves show the negative-norm modes of Eq. (19). The line styles and the various colors are used to ease the distinguishing of the different branches, with solid and dashed (dotted and dotted-dashed) curves showing left-moving (right-moving) modes. Thicker line styles have been used to highlight the three modes (and their conjugates) most relevant in this study. The extremal values of the curves are indicated by horizontal dotted lines. Here, their values are $D_{\max ,1}=34.9{\mathrm{mm}}^{-1}$ and $D_{\max ,2}=15.3{\mathrm{mm}}^{-1}$, while $\delta {\beta }_0=0.7{\mathrm{mm}}^{-1}$ (which corresponds to a soliton duration ${\tau }_0=10\mathrm{fs}$). Of key importance is the fact that the negative-norm branch (shown for $D<0$ by the thick dashed blue line) crosses the positive-norm one (solid thick black line) for $D_{\mathrm{cross}}=-28.4{\mathrm{mm}}^{-1}$, as this phase matching is responsible for the enhancement of the anomalous scattering coefficient.

Figure 3

Squared magnitudes of the scattering coefficients on a soliton of duration ${\tau }_0=10\mathrm{fs}$, with carrier frequency indicated by the black dot in the waveguide dispersion relation of Fig. 1. The vertical dotted line indicates the value of $D_{\min ,1}$ described in Sec. 3a (see also the left panel of Fig. 8 in Appendix pp1). On the left plot, the incident mode lives on the dotted red branch with positive group velocity. In that case, one faces an elastic linear mode mixing involving the probe and the idler (here described by a mode living on the solid black branch). We notice a sharp transition from pure transmission to total reflection occurring very near $D_{\min ,1}$. On the right plot, the incoming mode lives on the solid black branch, whose group velocity is negative and which is crossed by the negative-norm dashed blue branch in Fig. 2. Near $D_{\min ,1}$ we recover the rapid transition from transmission to total reflection, as in the left plot. We also see an enhancement of the transmission coefficient which is caused by the anomalous mode mixing involving the blue modes, whose ${\left|{\beta }_D\right|}^2$ is represented in dashed blue.

Figure 4

Scattering coefficients in the same setup as in Fig. 3, but now shown on a logarithmic scale so that the small coefficients are visible. The top row shows exactly the same scattering coefficients as Fig. 3, i.e., for an incident probe wave on the dotted red branch (upper left plot) and on the solid black branch (upper right plot). On the lower left plot, the incident probe wave lies on the dashed blue branch, while in the lower right plot it is on the light dotted bronze branch. The dotted vertical line indicates the value of $D_{\min ,1}$ in all plots except the lower right, where it shows $D_{\min ,2}$, the extremum of the “local” dispersion relation on top of the pulse at which the light solid green and light dotted bronze branches merge. Note that some scattering coefficients behave in a symmetrical way: in particular, the dashed blue curve on the upper right plot is very similar to the solid black curve on the lower left plot, while the dashed blue curve on the upper left plot is very similar to the dotted red curve on the lower left plot.

Figure 5

Variation of anomalous scattering coefficient with the duration of the pulse. The value of ${\left|{\beta }_D\right|}^2$ shown here describes the scattering of an incident mode on the solid black branch into an outgoing mode on the dashed blue branch. It is given as a function of $D$ for three different pulse durations ${\tau }_0$ [defined by Eq. (10)]: $40\mathrm{fs}$ (dotted curve), $20\mathrm{fs}$ (dashed curve), and $10\mathrm{fs}$ (solid curve). The corresponding total emission rates are found by integrating over $D$, and for the examples here are $1.1\times {10}^{-3}{\mathrm{mm}}^{-1},8.7\times {10}^{-3}{\mathrm{mm}}^{-1}$, and $7.5\times {10}^{-2}{\mathrm{mm}}^{-1}$. We thus see that the total rate scales approximately with ${\tau }_0^{-3}$, as can be rationalized using arguments presented in the text.

Figure 6

Variation of the observed spectrum with the duration of the pulse, assuming an ingoing vacuum state so that the observed photons are generated spontaneously. These photon-number spectra correspond to the same ${\left|{\beta }_D\right|}^2$ shown in Fig. 5, but plotted as a function of $\mathrm{\Delta }\omega$ rather than $D$, and such that their integrals over $\mathrm{\Delta }\omega$ are equal to the integral of ${\left|{\beta }_D\right|}^2$ over $D$. The left plot shows the spectrum observed on the (positive-norm) solid blue branch (with $\mathrm{\Delta }\omega <0$), while the right plot shows that on the (positive-norm) solid black branch (with $\mathrm{\Delta }\omega >0$). Note the significant difference in the scales used here, which is a result of the different slopes of these two branches in Fig. 2.

Figure 7

The real (left panel) and imaginary (right panel) parts of the dispersion relation $D\left(\mathrm{\Delta }\omega \right)$ for linear perturbations on top of a homogeneous background, for different values of $\gamma {\left(A_0^H\right)}^2$: (in order of decreasing line thickness) 0 (black), $2.5{\mathrm{mm}}^{-1}$ (red), and $5{\mathrm{mm}}^{-1}$ (blue). These values correspond to the nonlinear shift $\delta {\beta }_0$ of solitons with durations ${\tau }_0=5.3$ and $3.7\mathrm{fs}$, respectively. When $D$ is real, the solid and dashed curves correspond, respectively, to the positive- and negative-norm branches of the dispersion relation. In the right panel, we recover the standard modulation instability around $\mathrm{\Delta }\omega =0$, and the additional narrow instability where phase matching occurs.

Figure 8

The real (top row) and imaginary (bottom row) parts of the “instantaneous” dispersion relation $D\left(\mathrm{\Delta }\omega \right)$ at the peak of the soliton, as a function of $\mathrm{\Delta }\omega$. In the top row, we show in black (thicker line) the asymptotic dispersion relation, and in red (thinner line) that at the peak of a soliton of duration ${\tau }_0=10\mathrm{fs}$. Solid and dashed curves correspond to positive- and negative-norm modes, respectively, while dotted curves indicate where $D\left(\mathrm{\Delta }\omega \right)$ is complex. In the left panel, the extrema $D_{\max ,1}$ and $D_{\min ,1}$ are indicated, as are the other extrema $D_{\max ,2}$ and $D_{\min ,2}$. The right panel shows a zoom on the region near the crossing (the shaded rectangle on the left panel). On the bottom row are the imaginary parts of $D\left(\mathrm{\Delta }\omega \right)$ in the two regimes where $D\left(\mathrm{\Delta }\omega \right)$ becomes complex, at the peaks of solitons of three durations: (in order of decreasing line thickness) ${\tau }_0=10\mathrm{fs}$ (black), $20\mathrm{fs}$ (red), and $40\mathrm{fs}$ (blue). The left panel shows a descendant of the standard modulation instability on an inhomogeneous background, while the right panel shows the narrow instability that is the origin of the anomalous mode mixing studied in this paper. The width of the former is seen to vary as ${\tau }_0^{-1}$ while that of the latter varies as ${\tau }_0^{-2}$. Their heights instead follow the same scaling law, being both proportional to ${\tau }_0^{-2}$.

Figure 9

Evolution of the full fields $A(\tau ,z)=A_0(\tau ,z)\pm \delta A(\tau ,z)$, where $A_0$ initially contains the soliton field (of duration ${\tau }_0=10\mathrm{fs}$) and $\delta A$ contains a wave packet incident on the soliton. The mean field $A_0$ and the perturbation $\delta A$ are extracted by taking the half-sum and difference, respectively, of the two runs. The left column shows the real part of both $\delta A$ (in light gray) and $A_0$ (in thick black), in units where the peak of the initial soliton field $A_0(\tau ,z=0)=1$. On the right is shown the squared magnitude of their Fourier transforms in $\mathrm{\Delta }\omega$ space, again in units where the peak of ${\tilde{A}}_0(\mathrm{\Delta }\omega ,z=0)=1$. The rows correspond to different values of $z$: $0\mathrm{mm}$ (top row), $20\mathrm{mm}$ (middle row), and $40\mathrm{mm}$ (bottom row).

Figure 10

Conservation laws relevant to the scattering in Fig. 9. In the left column are shown the “instantaneous” squared magnitudes of the scattering coefficients, as described in the text. In the lower left panel is shown the “instantaneous” difference of the total norm from 1 [see Eq. (24)]. The larger deviations at intermediate $z$ arise due to the nonlinearities. In the right column are shown the variations in $N_A$ [see Eq. (7)] applied to the total field, and restricted to different regions of $\mathrm{\Delta }\omega$ space. The upper and middle panels show, respectively, the change in $N_A$ of the probe and of the soliton. The dashed curve in the middle panel takes into account both the effect of the probe on the soliton and the Cherenkov emission, while the solid curve shows the probe contribution only. In the lower right panel is shown the relative variation of the total photon number, integrated over all frequencies. One observes that it is constant up to one part in ${10}^{12}$, which indicates that it stems from numerical errors (given the much larger values of the scattering coefficients). In all panels, the dotted curves are for simulations with exactly the same initial conditions, but including a linear loss $\mathrm{\Gamma }$ such that $10\%$ of the energy is lost after $40\mathrm{mm}$. To compensate for the trivial steady decrease induced by the loss, we have multiplied the squared amplitudes by $e^{2\mathrm{\Gamma }z}$. As explained in the text, the residual effects visible in the plots are essentially accounted for by the reduction of the soliton power caused by $\mathrm{\Gamma }$.

Figure 11

Emission rates of spontaneously created photon pairs. We use the Fourier transform in $\tau$ of the linear perturbation $\delta A_{+}(z,\tau )$ generated in the upper component of the doublet when the lower component is fed by initial vacuum fluctuations [see the initial conditions given in Eqs. (B5)]. The soliton is removed by applying a window filter which vanishes for $\tau$ close to zero and is equal to 1 far from the soliton. We take an average growth rate of the spectrum with respect to $z$, and bin many discrete frequencies together to smooth out the profile (which is highly oscillatory for a single realization). Since the two peaks have very different spreads in frequency, we use two different binning rates: the discrete frequencies used in the simulation have a spacing of about $1.26\times {10}^{-4}\mathrm{PHz}$, and we bin these into groups of 40 for the left plot, and groups of 600 for the right plot. The black curves show the predictions from the linearized treatment at fixed $D$; they are exactly the solid curves already presented in Fig. 6.

Figure 12

Squared magnitude of the Fourier transform in both $z$ and $\tau$ of the linear perturbation $\delta A_{+}(z,\tau )$, for the same simulation used in Fig. 11 (i.e., where only the lower component of the doublet is fed by initial vacuum fluctuations). The soliton is removed by applying a filter which vanishes for $\tau$ close to zero and is equal to 1 far from the soliton. We use a logarithmic scale (with arbitrary units) to make most of the relevant data visible. As expected, the signal is found along the positive-norm dispersion relation, and is maximal near the mode crossing taking place for $|D_{\mathrm{cross}}|=28.4{\mathrm{mm}}^{-1}$.

Figure 13

Correlations between spontaneously produced photons when the initial state is vacuum. On the left is shown the magnitude of the first-order correlation of Eq. (B8) numerically obtained by integrating Eqs. (B2) and by extracting an average by the procedure described in the text. The normalization of the intensity is such that the maximal value is 1. The right panel shows the expected loci of the correlations in the $\left(\mathrm{\Delta }\omega ,\mathrm{\Delta }{\omega }^{\prime }\right)$ plane. The dotted lines show normal correlations involving products of the form ${\tilde{A}}^{★}\left(\mathrm{\Delta }\omega \right)\tilde{A}\left(\mathrm{\Delta }{\omega }^{\prime }\right)$, while the solid (and dashed) lines show anomalous correlations having the form $\tilde{A}\left(\mathrm{\Delta }\omega \right)\tilde{A}\left(\mathrm{\Delta }{\omega }^{\prime }\right)$. The dashed blue lines show where the correlations between the black and blue branches lie, and correspond to the strongest signals seen in the left panel.