# Field renormalization in photonic crystal waveguides

Pierre Colman

► To cite this version:

Pierre Colman. Field renormalization in photonic crystal waveguides. Physical Review A: Atomic, molecular, and optical physics [1990-2015], 2015, 92 (1), pp.013827. 10.1103/physreva.92.013827. hal-03543676

HAL Id: hal-03543676

https://hal.science/hal-03543676

Submitted on 25 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Field renormalization in Photonic Crystal waveguides

Pierre Colman*

DTU Fotonik, ørsted plads, Kgs. Lyngby, Denmark

IEF Institut d'Electronique Fondamentale, Université Paris-Sud, Orsay, France

(Dated: June 4, 2015)

A novel strategy is introduced in order to include variations of the nonlinearity into the nonlinear schrödinger Equation. This technique, which relies on renormalization, is in particular well adapted to nanostructured optical systems where the nonlinearity exhibits large variations up to two orders of magnitude larger than in bulk material. We show that it takes into account in a simple and efficient way the specificity of the nonlinearity in nanostructures that is determined by geometrical parameters like the effective mode area and the group index. The renormalization of the nonlinear Schrödinger Equation is the occasion for physics oriented considerations and unveils the potential of Photonic Crystal waveguides for the study of new nonlinear propagation phenomena.

Keywords: Photonic Crystal, Nonlinear optics, Slow-light, NLSE

# I. INTRODUCTION

Nanostructured optical systems offer a great plasticity because most of their optical properties can be engineered through the design [1–13]. Thus the control by the geometry allows one to get on demand properties, precisely adapted to the phenomena one wants to investigate [7, 14–16]. This is both interesting for applications where it is then possible to optimize all optical functions to a large extent [17] and for fundamental investigations where the experimental conditions could be precisely set to maximize a given effect.

This plasticity comes usually along with large variations of the optical properties over a small bandwidth because effects of geometry depend highly on the actual wavelength size. For instance in Photonic Crystal waveguides, variations of the nonlinear effective Kerr effect up to 100% can occur over a 10nm bandwidth [18]. Consequently the modeling of the nonlinear propagation of optical pulses in such structures is essential for understanding precisely the interplay between the different effects taking place; such a task is challenging [19–21].

The use of the generalized Nonlinear Schrödinger Equation (GNLSE) is usually preferred to direct nonlinear FDTD simulation, though that latter being more accurate [22, 23], because it has very low computation burden and also provides a direct link between the effective coefficient in the GNLSE and the phenomena that are observed. Namely it is straightforward to add phenomena specifically found in semiconductor optics, like for instance the effect of nonlinear absorption and free carriers [24–27]. Moreover discussion and interpretation are made easier as the coefficients in the GNLSE are derived into effective lengths like the dispersion length $L_d = T_0^2/\beta_2$, the nonlinear length $L_{NL} = 1/(\gamma P_0)$ or the shock formation length $L_{shock} = 0.43L_{NL}/|\tau_{NL}|$ [28], etc. The effective lengths [29] provide a rapid overview on the relative strength of the different competing effects. However, because

the GNLSE relies on several approximations, it is important to check that the physics governing the optical system is still correctly described.

Higher order nonlinear effects, for weak perturbations, are taken into account by adding successive corrective terms to the initial NLSE. Here we focus on the inclusion of variations of the nonlinear Kerr response with the angular frequency. Using the transformation $\Delta\omega \rightarrow i\partial_t[30]$, the Taylor expansion of the Kerr coefficient (hence its variation in the frequency domain) can be included into the GNLSE (time domain). Such decomposition -limited to the first order- takes into account self-steepening effects and can be used to model the formation of optical shock front [29, 31].

Despite being easy to implement, such technique suffers from two major limitations.

First it uses the derivation of an aggregated effective parameter computed for a given frequency, whereas by essence nonlinearity involves also the interactions between several waves at different frequencies. Notably, not only the effective mode area, but also the variations of the modes overlap -directly related to cross phase modulation (XPM)- are important in integrated optics. A simple Taylor expansion of the aggregated Kerr coefficient cannot take into account that latter effect. This issue was addressed in Ref. [32] where it was demonstrated that the physics related to variations of the mode area overlap could be in principle retrieved by the addition of extra operators; and that such corrections could actually have a noticeable impact on pulse propagation.

The second limitation is due to the nature of the Taylor expansion itself: it is intended for weakly varying functions and convergence tends to be slow as soon as the functions exhibit large or non-monotonic variations. In case of large variations of the nonlinear effective coefficients, the best option would then be to split the total bandwidth into sub-domains wherein variations of the nonlinear properties are negligible [21, 33]. Thus each sub-domain is described with good accuracy using a GNLSE; then the different equations are connected through ad-hoc parameters that describe with precision

<footnote>* pierre.colman@espci.org</footnote>
effects like cross phase modulation (XPM) or four wave mixing (FWM). This method is especially well suited to describe FWM because for most practical application pumps, signal and idler spectra are well separated; moreover the small number of beams involved allow setting up a set of coupled mode equations with only a few equations. However, if the spectrum is continuous over a large bandwidth, splitting the initial simulation domain into several pieces is complex because there is no obvious choice for the separation frequency between the sub-domains; and numerical artefact might appear at the junction between two sub-domains.

Consequently neither of the two methods gives satisfactory results in case the single pulse bandwidth extends over large variations of the nonlinear effective parameters. This is especially true for dispersion-engineered Photonic Crystal waveguides (PhCWGs) [18] where the nonlinearity arises mainly from the slow light enhancement [19]. Indeed any changes of the group index is immediately reported through the slow light factor $S^2 = (n_g/n_0)^2$ [34] on the Kerr nonlinearity, which is expressed as $\gamma_{eff} = \omega n_{2I}/(cA_{eff}) * S^2$. Variations of the modal area are also important in PhCWGs, but still much weaker than those associated with the slow light variation. Note that while S ranges from 2 to 10 in PhCWGs [24], its effect is less sensitive in Nanowires on silicon which in contrast are much more impacted by variations of the modal area. The slow light effect in Nanowires is roughly $S \approx 1.2$ [35] and is in fact already implicitly taken into account in the expression of the effective area [36].

The method we present here consists in defining and solving the GNLSE for another field than the power flux (cf. Fig. 1). The point is to find the pseudo field that leads to the GNLSE exhibiting minimal variation of the nonlinear coefficient. Especially this technique is well-suited for nanostructured optical systems where the nonlinearity is defined by the effective area and the group index. Previous studies focused on the ab initio definition of the effective parameters of the propagating equation once the eigen modes profile is known. It was demonstrated that the presence of slow light must be taken into account in the normalization in order to guarantee that the correct energy flow through the medium is preserved despite the presence of strong material dispersion or slow light [37]. This first normalization step is crucial. However, the resulting nonlinear propagation equation will contain also a few pre-factors; this aspect is not important for few modes problems like FWM or parametric generation where one can define a set of coupled equations, but becomes problematic for continuous mode problems. For example the propagation equation can be turned back into the form of a classical NLSE, but at the cost of extra approximations which ultimately alter the accuracy of the numerical solution. The second step consists then in the renormalization of the energy flow into a variable whose nonlinear propagating equation has a simpler form. Our finding is that instead of computing

the propagation of the power flux, it is better to weight it by the photonic density of states. The group index plays then an important role in this second normalization step because the group velocity corresponds in fact to the first derivative of the band diagram that is directly linked to the optical density of state. To be more precise, the optical density of state is proportional to inverse of the group velocity and to the inverse of the volume occupied by the photons [38, 39]. Considering that the effective area in PhCWGs does not vary much compared to the group index, then a variation of the group index could be interpreted directly as a variation of the optical density of state.

First, we will first review how the slow light effect is introduced into the nonlinear propagation equations; and thus insist on the two contributions of slow light on the nonlinear enhancement. We will then introduce our renormalization technique and apply it to the computation of nonlinear pulse propagation in dispersion-engineered PhCWGs.

## II. NLSE

We start from the Maxwell's equations. As we are dealing with unidirectional guided propagation (i.e. waveguides), we have split the spatial coordinates into transverse coordinates $\vec{r}_{\perp}$ and $z$ along the waveguide direction.

$$ \nabla \times \vec{E}(\vec{r}_{\perp}, z, t) = -\mu_0 \partial_t \vec{H}(\vec{r}_{\perp}, z, t) \quad (1) $$

$$ \nabla \times \vec{H}(\vec{r}_{\perp}, z, t) = \epsilon(\vec{r}_{\perp}, z) \partial_t \vec{E}(\vec{r}_{\perp}, z, t) + \partial_t \vec{P}_{NL}(\vec{E}) \quad (2) $$

When $\vec{P}_{NL} = 0$, the above equations accept a set of eigen solutions in the frequency domain:

$$ \vec{e}(\vec{r}_{\perp}, z, \omega) = \vec{e}_{bloch}(\vec{r}_{\perp}, z - n.a., \omega)e^{-i\omega t + iK(\omega)z} \quad (3) $$

$$ \vec{h}(\vec{r}_{\perp}, z, \omega) = \vec{h}_{bloch}(\vec{r}_{\perp}, z - n.a., \omega)e^{-i\omega t + iK(\omega)z} \quad (4) $$

with $n \in Z$ and a being the PhC lattice parameters. $K(\omega)$ corresponds to the propagation constant at the frequency $\omega$. For periodic waveguides like PhCWGs, the electric and magnetic field distribution is determined by a Bloch mode which is a periodic function with periodicity a. Hereafter, the dependence on $(\vec{r}_{\perp}, z)$ of the eigen and Bloch modes is implicitly assumed and we will use the short notation $\vec{e}(\omega)$, $\vec{h}(\omega)$, $\vec{e}_{Bloch}(\omega)$, $\vec{h}_{Bloch}(\omega)$. Similarly as detailed in [40], the forward propagating equation is obtained by multiplying eqs. (1-2) by the solution of the unperturbed ($P_{NL} = 0$) system (i.e. $e(\omega)^*$ and $h(\omega)^*$). Subtracting these two new equations from each other leads, after some algebra including the identity $\vec{A} \cdot \nabla \times \vec{B} = \vec{\nabla} \cdot (\vec{A} \times \vec{B}) + \vec{B} \cdot \nabla \times \vec{A}$ to [41]:

$$ \vec{\nabla} \cdot (\vec{E} \times \vec{h}(\omega)^* - \vec{H} \times \vec{e}(\omega)^*) = -\vec{e}(\omega) \partial_t \vec{P}_{NL}(\vec{E}) - \partial_t (\epsilon_0 \vec{E} \cdot \vec{e}(\omega) + \mu_0 \vec{H} \cdot \vec{h}(\omega)) \quad (5) $$
We will now look for a perturbative solution, in the frequency domain, of the form

$$ \vec{E}(\vec{r}_{\perp}, z, \Delta\omega) = A(z, \Delta\omega)e^{-i\Delta\omega t}\vec{e}_{\text{Bloch}}(\omega)e^{-i\omega_{\text{ref}}t+iK(\omega_{\text{ref}})z} \quad (6) $$

$$ \vec{H}(\vec{r}_{\perp}, z, \Delta\omega) = A(z, \Delta\omega)e^{-i\Delta\omega t}\vec{h}_{\text{Bloch}}(\omega)e^{-i\omega_{\text{ref}}t+iK(\omega_{\text{ref}})z} \quad (7) $$

Where $\Delta K(\omega) = K(\omega) - K(\omega_{\text{ref}})$ and $\Delta\omega = \omega - \omega_{\text{ref}}$. Note that the phasor term in eqs. (6-7) is different from the one found in eqs. (3-4), hence we would like to compute directly the evolution of the pulse envelope

centered at frequency $\omega_{\text{ref}}$. Also we assume the envelope $\partial_z A/A \ll 1/a$ varies slowly, so $\partial_z A(z, \omega)$ is constant over a single PhC lattice. We inject eqs. (6-7) into eq. (5) and integrate over a unit PhC cell; this leads to:

$$ \begin{aligned} & (\partial_z A(z, \Delta\omega) - i\Delta K(\omega)A(z, \Delta\omega)) \iiint_{\text{cell}} \vec{z} \cdot (\vec{e}_{\text{Bloch}}^*(\omega) \times \vec{h}_{\text{Bloch}}(\omega) - \vec{h}_{\text{Bloch}}^*(\omega) \times \vec{e}_{\text{Bloch}}(\omega)) dr^3 \\ & = i\omega \iint_{\text{cell}} \vec{e}_{\text{Bloch}}(\omega)^* \cdot \vec{P}_{\text{NL}}(\Delta\omega, \vec{E}) dr^3 \end{aligned} \quad (8) $$

$A(z, \Delta\omega)$ is defined with regard to a central frequency reference $\omega_{\text{ref}}$. Thus eq. (8) gives direct access to the envelope of the pulse, without the fast oscillations in space and time $\{\omega_{\text{ref}}, K(\omega_{\text{ref}})\}$. The integral in the left hand side of eq. (8) corresponds to twice the integral of the Poynting vector $\vec{\Pi}_{\omega}$ over the PhC cell. The link between the Poynting vector and the Bloch mode energy is set through the relationship:

$W_{\omega}$ is a normalization factor, which corresponds to the Bloch mode energy $1/2 \iint_{\text{cell}} \epsilon |e_{\text{Bloch}}(\omega)|^2 + \mu |h_{\text{Bloch}}(\omega)|^2 dr^3$. Finally we have

$$ \begin{aligned} \iint_{\text{cell}} \vec{z} \cdot \Re(\vec{e}_{\omega}^* \times \vec{h}_{\omega}) dr^3 &= \iint_{\text{Surface}} \vec{z} \cdot \vec{\Pi}_{\omega} dr^2 * a \\ &= v_g(\omega)W_{\omega} \end{aligned} \quad (9) $$

$$ \partial_z A(z, \Delta\omega) = i\Delta K(\omega)A(z, \Delta\omega) + i \frac{\omega}{2v_g(\omega)W_{\omega}} \iint_{\text{cell}} \vec{e}_{\text{Bloch}}(\omega) * \vec{P}_{\text{NL}}(\Delta\omega, \vec{E}) dr^3 \quad (10) $$

Before going further into details, let us focus an instant on the peculiar normalization choice for both $A(z, \omega)$ and $W_{\omega}$. In the linear regime, the choice is to preserve the power flux, so the computed field directly corresponds to the energy flowing through the medium. Using eq. (9), one finds that the total energy transmitting through the waveguide at frequency $\omega$ is $P(z, \omega) = |A(z, \omega)|^2 v_g(\omega) W_{\omega}/a$. As a result the Bloch modes are normalized such as $v_g(\omega) W_{\omega}/a = 1$; so one gets directly from the NLSE $P = |A|^2$. However such normalization is done independently of the nonlinear problem which is considered; consequently it may not necessarily be the

best choice from a pure numerical point of view.

Regarding the nonlinear polarization $P_{\text{NL}}$, we assume that it is due to the $\chi^{(3)}$ response expressed as: $P_{\text{NL}}(r, \omega_0) = 3/2 \epsilon_0 \chi_{111}^{(3)} \int (\vec{E}^*(r)_{\omega_1} \cdot \vec{E}(r)_{\omega_2}) \cdot \vec{E}(r)_{\omega_3} \delta(\omega_2 + \omega_3 - \omega_1 - \omega_0) d\omega_3 \quad \text{, where } \delta \text{ stands for the Dirac delta function. The exact form of } P_{\text{NL}} \text{ is directly related to the nonlinear tensor [42]; and consequently will have a different formulation depending on the material that is actually considered. However the conclusions presented here are general and could be easily extended to any peculiar form of } P_{\text{NL}}. \text{ Finally, we}$
obtain

$$ \partial_z A(z, \omega_0) = \imath K(\omega_0) A(\omega_0) + \frac{\imath \omega_0 n_{2I}}{c A_{eff}(\omega_0)} \frac{\sqrt{n_g \omega_0}}{n_0^2} \int E(z, \omega_1)^* E(z, \omega_2) E(z, \omega_3) g h_{\omega_0}(\omega_1, \omega_2, \omega_3) \delta(\omega_2 + \omega_3 - \omega_1 - \omega_0) d\omega_i^3 \quad (11) $$

We have defined $n_{2I} = 3\chi^{(3)}/(4c\epsilon_0 n_0^2)$. We also introduced the error function $gh_\omega$ which takes into account the fact that the mode overlap integrals may deviate from $A_{eff}(\omega_0)$ for frequency $\omega_i$ different than $\omega_0$. Hence

$$ A_{eff}(\omega_0) = \frac{\left( \iiint_{Cell} \epsilon(r) |e_{\omega_0}|^2 dr^3 \right)^2}{an_0^4 \iiint_{NLmat.} |e_{\omega_0}|^4 dr^3} \quad (12) $$

$$ gh_{\omega_0}(\omega_1, \omega_2, \omega_3) = \frac{\iiint_{NLmat.} (\vec{e}_{\omega_0}^* \cdot \vec{e}_{\omega_2})(\vec{e}_{\omega_1}^* \cdot \vec{e}_{\omega_3}) dr^3}{\iiint_{NLmat.} |e_{\omega_0}|^4 dr^3} \quad (13) $$

NL\_mat indicates that the integral is performed over the nonlinear material in the PhC unit cell. At this point, one notes that the nonlinearity depends on the intensity of electric field E while in eq. (11) the nonlinear pulse evolution is set for the power flux $P = |A|^2$ [19]. Eq. (9) is then used once more and one finally gets:

$$ \partial_z A(z, \omega_0) = \imath K(\omega_0) A(\omega_0) + \imath \frac{\omega_0 n_{2I}}{c A_{eff}(\omega_0)} \frac{\sqrt{n_g \omega_0}}{n_0^2} \int A(z, \omega_1)^* A(z, \omega_2) A(z, \omega_3) g h_{\omega_0}(\omega_1, \omega_2, \omega_3) \sqrt{n_g(\omega_1)n_g(\omega_2)n_g(\omega_3)} \delta(\omega_2 + \omega_3 - \omega_1 - \omega_0) d\omega_i^3 \quad (14) $$

Unfortunately, the $\sqrt{n_g(\omega_{i=1,2,3})}$ and $gh_{\omega_0}(\omega_1, \omega_2, \omega_3)$ terms in the right hand side depend on $\omega_{i=1,2,3}$, hence these prefactors cannot be set aside the integral. How-

ever, the formulation of a standard NLSE equation - expressed in the frequency domain- can be retrieved assuming that $gh_\omega = 1$ and that $\sqrt{n_g(\omega_i)}$ do no vary:

$$ \partial_z A(z, \omega_0) = \imath K(\omega_0) A(z, \omega_0) + \imath \frac{\omega_0 n_{2I}}{c A_{eff}(\omega_0)} \frac{n_{g\omega_0}^2}{n_0^2} \int A(z, \omega_1)^* A(z, \omega_2) A(z, \omega_3) \delta(\omega_2 + \omega_3 - \omega_1 - \omega_0) d\omega_i^3 \quad (15) $$

Note that eq. (15) could have also been directly obtained through the derivation of eq. (14) for a monochromatic wave propagation. In fact this is most of the time how the effective nonlinear coefficient of the NLSE are derived. But the purpose here is precisely to point out explicitly the frequency dependence of the different effective parameters; and the approximations done in regard with it. The propagating equation corresponding to eq. (15) in the time domain is:

$$ \partial_z A(z, t) = \imath D(\imath \partial_t) A(z, t) + \imath \gamma_0 |A(z, t)|^2 A(z, t) \quad (16) $$

Here $\gamma_0 = \gamma(\omega_{ref})$ is the effective Kerr nonlinearity with $\gamma(\omega) = n_{2I}\omega/(cA_{eff}(\omega)).(n_g/n_0)^2$. The dispersion operator $D(\imath\partial_t) = \sum_{n \ge 2} (\partial_\omega^n k)(\imath\partial_t)^n/n$ accounts for dispersion at all orders, with t being the retarded time in the moving frame at velocity $c/n_{g-\omega_0}$. This NLSE does not include any variations of the effective nonlinearity over the simulation domain. Effects related to the dispersive nonlinearity are then re-introduced using the first
order perturbation $\gamma_1$ which is accordingly defined as
$\gamma_1 = \partial_{\omega}\gamma(\omega_{ref})$.

$$
\begin{equation}
\begin{split}
\partial_z A(z,t) = {}& iD(i\partial_t)A(z,t) \\
& + i(\gamma_0 + \gamma_1 i\partial_t)|A(z,t)|^2 A(z,t)
\end{split}
\tag{17}
\end{equation}
$$

Eq. (17) is derived thanks to the fact that the different $\sqrt{n_g(\omega_i=\{1,2,3\})}$ factors in eq. (14) are set outside the integral sign. Such an operation is only valid if the group index $n_g$ does not depend on the different frequencies $\omega_i$ that contribute to the nonlinear effect. This (mathematical) approximation is equivalent to the hypothesis according which the strength of nonlinearity only depends on the angular frequency which the nonlinear effect takes place at; hence not on the different frequencies contributing to the nonlinear effect. There is no evidence that this assumption is valid in general. Especially when it comes to nanostructures, one may think that the photons that are the most confined contribute the most to the nonlinearity. Although the equi-contribution hypothesis does not hold in general, it can happen that the $\sqrt{n_g(\omega_1)n_g(\omega_2)n_g(\omega_3)}$ prefactor depends only on $\omega_0$. Indeed the different frequencies $\omega_i$ are not independent, but must satisfy the energy conservation condition $\omega_0 + \omega_1 = \omega_2 + \omega_3$. In such a case the equi-contribution approximation still describes the nonlinear photon dynamics with accuracy; and we will refer to such photons as being dispersive photons.

### III. RENORMALIZATION

As depicted by Fig 1 we shall try to circumvent the problem caused by fluctuations of the nonlinear effective parameters in eq. (14) by finding a proper referential wherein the nonlinear response is flat. Looking back at eq. (14) tells us that this is indeed the key in order to obtain a time-domain propagation equation like eq. (16) : the different frequencies should have an **equi-contribution** to the nonlinear process; hence no frequency dependent pre-factor must appear inside the integral sign of eq. (14). To do so, we decided to weight the frequencies by a ad-hoc $m(\omega)$ contribution.
Note that we neglect at first any variation of the mode field distribution over the pulse bandwidth ($gh_{\omega_0}(\omega_i) = 1 \forall \omega_0$). The reason for this is that variations of the group index in PhCWGs account for about 75% of the total variation of the effective nonlinear coefficients [18, 21]. Consequently, dealing with the variations of the slow-light factor would be the first step. Moreover this simplified case is a good test case to present the renormalization technique. We will show in next section how the variations of the mode field distribution can be taken as well into account by this technique.

It appears that a natural choice would be to solve the propagation equation for the field $\Psi(\omega) = \sqrt{n_g(\omega)/n_0}A(\omega)$ so eq. (14) becomes

$$
\begin{equation}
\partial_z \Psi(z, \omega_0) = i K(\omega_0) \Psi(z, \omega_0) + i \frac{\omega_0 n_{2I}}{c A_{eff}(\omega_0)} \frac{n_g \omega_0}{n_0} \int \Psi(z, \omega_1)^* \Psi(z, \omega_2) \Psi(z, \omega_3) \delta(\omega_2 + \omega_3 - \omega_1 - \omega_0) d\omega_i^3 \tag{18}
\end{equation}
$$

<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Schematic of the two methods used here to compute the nonlinear pulse propagation in PhC waveguide. a) Standard method: the GNLSE is solved for the energy flux. b) Input conditions (and the GNLSE) are renormed, and the nonlinear propagation is solved for a pseudo-electric $\Psi$ field.</figcaption></figure>

Instead of solving the nonlinear propagation of the power flux $P = |A|^2$, we are now solving the propagation of a pseudo-field $\Psi$. Looking in details at the differences between eq. (18) and eq. (15), we see that the slow light enhancement factor enhancement $S^2 = (n_g/n_0)^2$ has been replaced by $S$. Interestingly, if we had intro-

duced higher order nonlinear effects such as three photons absorption (ThPA) which are associated to a slow light enhancement of $S^3$, then the slow light pre-factor would have been turned into $S$ as well. A first feature is that the strength of the nonlinearity appears to be much weaker. Besides if we consider the variations of the Kerr effect $\partial_{\omega}\gamma(\omega)$, the contribution of the slow-light to the characteristic self-steepening time [18] $\tau_{NL} = \partial_{\omega}\gamma(\omega)/\gamma(\omega)$ has been halved! Consequently eq. (18) exhibits weaker nonlinearity and even weaker relative nonlinear dispersion than eq. (15). Weaker relative variations of the nonlinearity means that the actual nonlinear variations could be taken into account with more accuracy by eq. (18)- which could also include second order perturbative corrections-. Note that although the nonlinearity appears much weaker, the input field has been renormalized as well and is now much stronger. Consequently global parameters like the soliton number are preserved.

Usually, such a strategy is not convenient because it requires additional transformation back and forth between computed quantities ($\Psi \neq A$) and measured ones
(power flux $P = |A|^2$) -cf. Fig. 1. Also one could wonder whether such a renormalization would not simply lead to unphysical solutions, hence simpler equations but not describing correctly the physics. Any renormalization could be applied to eq. (14), as long as it is done in consistency with the math -e.g. the integral sign in eq. (14)-; just that not any renormalization would actually lead to a simpler formulation. Regarding the specific choice of $m(\omega) = \sqrt{n(\omega)}/n_0$, eq. (9) tells that because $n_g|A|^2 \propto |E|^2$, then the pseudo-field $\Psi$ is actually directly proportional to the Bloch mode electric field. This is consistent with the fact that the electric field density is the physical quantity that matters for nonlinearity, not the power flux. Because the group velocity governs most of the physics (slow-light enhancement) in PhCWGs, it is then not surprising that it plays a role both in the initial normalization of the Bloch field and also here in the renormalization of the computed power flux.

In brief, the renormalization technique poses the question whether the natural choice that is usually made to compute the evolution of the power flux is right. We think it is not for nanostructured systems.

## A. Implementation and comparison

Nonlinear pulse propagation can now be dealt with in two different ways. On one hand, one can use eq. (17), expanding eventually the Taylor series of the nonlinearity beyond the first order in order to get a better match; on the other hand eq. (18) will also provide an accurate result, and might be easier to implement. The outcome of both equations should be about the same, given that they are indeed describing the same physical system.

To investigate the differences between eq. (17) and eq. (18), we take as a test case the nonlinear pulse propagation experiment performed in conditions similar to those in Ref. [16]. A 2.3 ps Fourier limited pulse is send close to the zero group velocity wavelength (ZDW) of an 1.5mm-long dispersion-engineered PhC waveguide. The pulse peak power is 8W and corresponds to a soliton number of $N = 2.1$. For consistency with real systems, the equations have been adapted to include the specificity of optical semiconductor nanostructures, mainly the effect of nonlinear absorption and the presence of free carriers [24–27]. The power flux propagation is computed by mean of the following generalized Schrödinger equation[43].

$$ \frac{\partial A}{\partial z} = -\frac{\alpha}{2}A - \alpha_3|A|^4A - D(\iota\partial_t)A + \iota(\gamma_0 + G(\iota\partial_t))|A|^2A - (\sigma + ik_0\delta)NA \quad (19) $$

$\alpha_0 = 2dB/mm$ and $\alpha_3 = 25/(W^2.mm)$ [26] stand for the linear propagation loss and three photons absorption (ThPA). The dispersion operator $D(\iota\partial_t) = \sum_{n \ge 2}(\partial_n^k)(\iota\partial_t)^n/n$ accounts for dispersion at all orders, with t being the retarded time in the moving frame at velocity $c/n_g$ (calculated at the input wavelength).

Besides, we introduce the Kerr operator $G(\iota\partial_t) = \sum_{n \ge 1}(\partial_n^\gamma)(\iota\partial_t)^n/n!$ to take into account the dispersion of the Kerr coefficient with the angular pulsation (higher order shock terms). Such expansions are intended to provide accurate numerical results, though it limits the insight into the physical parameters that govern the pulse propagation. $\sigma$ and $\delta$ account for the free carriers absorption and dispersion respectively. Owing to the fact that we are considering here a high bandgap material like GaInP that exhibits solely three photons (ThPA) and no two photons absorption, the self generated plasma does not impact much the overall dynamics. Practical details related to the way the free carriers effects are computed and added to the NLSE are found in Ref. [43–45].

<figure><img src="image_3.png" /><figcaption>FIG. 2. Parameters used for the simulation. a) group index $n_g$ as function of the wavelength. b) Effective Kerr coefficient. c) Effective three photons absorption (ThPA) coefficient. The square mark indicates the position of the input wavelength in Fig 3</figcaption></figure>

As a guideline we have $\beta_2 = -6.7ps^2/mm$, $\beta_3 = -1.7ps^3/mm$, $\gamma_0 = 2200/(W.m)$ and $\gamma_1/\gamma_0 = -170fs$, the waveguide dispersion and the variation of the nonlinear coefficients (Kerr and ThPA) with the angular frequency are shown in Fig 2.

Briefly, the renormed NLSE equation is obtained from eq. (17) by (i) expressing it in the frequency domain, (ii) dividing the nonlinear coefficient by $m(\omega)^n$ where $n = 2$ for $\chi^{(3)}$ effect (e.g. Kerr) and $n = 4$ for $\chi^{(5)}$ (e.g. Three photons absorption), (iii) the input field $\Psi$ is then obtained by multiplying the initial input field $A(\omega)$ by $m(\omega)$.

Comparison between the results of eq. (17) (grey) and its counterpart solved for the pseudo electric field $\Psi$ (thick black) are shown in Fig 3. The input soliton has undergone a blue Soliton-Self Frequency Shift (SSFS) of a few nm and dispersive waves are generated in the normal dispersion region [46, 47]. In the temporal domain the complex interplay [18] of the dispersive nonlinearity [48], free carriers effects [45], and the SSFS, results in a pulse advance of a few picoseconds. We see the spectral
<figure><img src="image_4.png" /><figcaption>FIG. 3. Comparison between the two numerical models. a) Spectra. b) Temporal trace. Thin black: input pulse; thick black: with renormalization $m = \sqrt{n_g(\omega)/n_0}$. Thick grey: without renormalization. $\lambda_0$ indicates the position of the zero group velocity dispersion point (ZVD).</figcaption></figure>

position of the dispersive wave (in the normal dispersion region) and the time advance (10ps vs. 15ps) differ between the two models. Indeed after the initial generation, the dispersive wave interacts with the soliton through cross-phase modulation and cascaded four wave mixing [49]: therefore its amplitude and position depends greatly on the exact form of the Kerr nonlinearity. The amplitude of the dispersive waves grows much stronger when the renormalization procedure is employed; and because the spectral recoil appears in reaction to the emission of the dispersive wave, the SSFS is much stronger for the renormalized case. Divergences between the two numerical models are clearly visible.

## B. Comparison with an analytic solution

Now that we have shown that the two approaches (the nominal GNLSE and its renormalized counterpart) lead to different results, we must determine which model is (the most) correct. The case that we just discussed corresponds to a realistic case where both the medium (dispersion, nonlinearity, absorption) and the input parameters (power, duration) are within reach of current experiments. So it would be in principle possible to perform such an experiment and compare the two models with the experimental results. However such measurements are not available yet. Another way would be to confront directly the results of the two models with an analytic case.

Recent work by Erkintalo et al. [49] demonstrated that nonlinear pulse propagation (continuous spectrum) and cascaded FWM (discrete spectrum) are closely related and that the nonlinear pulse dynamics can be described as a cascade of FWM events. This means that the capacity of an equation to accurately reproduce the reality

is intrinsic to its capacity to deal correctly with FWM. Although no analytic solution exists for the exact case we just studied we can still simplify our problem to the propagation in a lossless and $L = 100\mu m$ short PhCWG of a $P_0 = 8W$ single continuous wave beam of same central frequency ($\omega_{ref}$) as previously. In such a situation the short length of the waveguide and the absence of propagation loss render the undepleted pump approximation valid. The FWM conversion efficiency [7] depending on the pump signal detuning ($\delta\omega$) is then expressed through the analytic formulation:

$$ \eta(\delta\omega) = (\gamma_{FWM}(\delta\omega)P_0L)^2 \left( \frac{\sinh(g(\delta\omega)L)}{g(\delta\omega)L} \right)^2 \quad (20) $$

$$ g^2(\delta\omega) = (\gamma_{FWM}(\delta\omega)P_0)^2 - (\Delta K_{L\omega Ref}(\delta\omega) + \Delta K_{NL\omega Ref}(\delta\omega))^2/4 \quad (21) $$

$$ \Delta K_{L\omega Ref}(\delta\omega) = 2K(\omega_{ref}) - K(\omega_{ref} + \delta\omega) - K(\omega_{ref} - \delta\omega) \quad (22) $$

$$ \Delta K_{NL\omega Ref}(\delta\omega) = 2P_0(\gamma_{XPM}(\delta\omega) + \gamma_{XPM}(-\delta\omega) - \gamma_{SPM}(\omega_{ref})) \quad (23) $$

$\gamma_{FWM}$, $\gamma_{XPM}$, $\gamma_{SPM}$ account respectively for the FWM nonlinear coupling coefficient, the cross phase modulation (XPM) between the strong pump and the weak signal/idler and the self-phase modulation (SPM) of the pump [21]. Especially it takes into account the different modes fields overlap. According to the notation used previously, these effective nonlinear coefficients are expressed as:

$$ \gamma_{SPM}(\omega_{ref}) = \frac{\omega_{ref} n_{2I}}{c A_{eff}(\omega_{ref})} g h_{\omega_{ref}}(\omega_{ref}, \omega_{ref}, \omega_{ref}) \left( \frac{n_g(\omega_{ref})}{n_0} \right)^2 \quad (24) $$

$$ \gamma_{XPM}(\omega) = \frac{\omega n_{2I}}{c A_{eff}(\omega)} g h_{\omega_{ref}}(\omega_{ref}, \omega, \omega_{ref}) \left( \frac{n_g(\omega_{ref}) n_g(\omega)}{n_0^2} \right) \quad (25) $$

$$ \gamma_{FWM}(\delta\omega) = \frac{(\omega_{ref} + \delta\omega)n_{2I}}{c A_{eff}(\omega + \delta\omega)} g h_{\omega_{ref}}(\omega_{ref} - \delta\omega, \omega_{ref}, \omega_{ref}) \left( \frac{n_g(\omega_{ref})\sqrt{n_g(\omega_{ref} - \delta\omega)n_g(\omega_{ref} + \delta\omega)}}{n_0^2} \right) \quad (26) $$

These coefficients are computed for each pump-signal detuning $\delta\omega$; as a result we get the FWM gain curve as shown in thick plain black in Fig 4. We compare now this analytic curve to the results given by the different models.

First we see that the GNLSE (thick dashed green) does not appear to converge any better than the standard NLSE (light dashed blue). The NLSE does not take
<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) $\eta = P_{\text{idler}}(L)/P_{\text{signal}}(0)$ for various pump-signal detuning after a propagation of $L = 100\mu m$. Thick Black: Analytic model. Thin dashed blue: results of the NLSE. Thick dashed green: GNLSE. Plain grey: GNLSE renormalized wherein only the group index variations are taken into account. cyan dots: GNLSE wherein both slow light and effective area variations are included into the renormalization.</figcaption></figure>

into account any variations of the nonlinearity, hence $\gamma_{SPM} = \gamma_{XPM} = \gamma_{FWM}$. While the NLSE tends to overestimate the FWM bandwidth, the GNLSE underestimates the FWM gain. Thus the inclusion of self-steepening (i.e. the dispersion of the nonlinearity) does not appear as a great improvement. One must note that the GNLSE is still essential to model effect specifically related to self-steepening like the formation of a shock front or to explain the energy dependent time-advance of the nonlinear pulse [18, 48]. We see that the renormalization of the slow light enhancement factor (plain grey) improves the convergence of the GNLSE: for small detuning ($\delta\lambda < 2nm$) the GNLSE and the analytic model converge; for larger detuning some discrepancies still remain, but the overall error is still weaker than for the un-renormalized GNLSE or the simple NLSE. The remaining error is due to the fact that the effect of the variation of the effective modal area is neglected at first. We see that if we then renormalized the GNLSE to take as well into account that later effect (cyan dots), we obtain a very good agreement between the analytic model and the renormalized GNLSE. We present in the last section how the renormalization technique is general and can take into account variation of the modal area.

This demonstrates that the renormalization technique that we present here really constitutes an improvement compared to previous formulations. More generally our technique could also be seen as a generalization of Ref. [21] that deals in a seamless way with the generation of multiple signal and idler orders, hence not only limited to a discrete set of a few (usually four) beams.

### C. Discussion

The differences between the two models (with or without the renormalization) indicate that the physics governing the nonlinear pulse is necessarily different. Especially, one of the largest changes between eq. (17) and eq. (18) lies in the ratio of Self-Phase (SPM) to Cross-Phase Modulation (XPM) intensity. Usually the XPM is twice as many times stronger as the SPM for a given frequency. After the renormalization of the slow light variations, the photons are weighted by the quantity $m(\omega)^2 = n_g(\omega)/n_0$ and they do not have the same contribution to the nonlinear index change. The weight factor corresponds to the increase of the photonic density of states in the PhCWGs compared to the bulk material. As said earlier, the first derivative of the band diagram could be both interpreted as the group velocity or the density of states: photons with a higher density of states will have a higher probability to interact with the nonlinear medium. By analogy, the inclusion of the variations of the effective modal area would correspond to dividing the photon eigen contribution by an effective volume. Thus the enhancement of the nonlinearity due to the slow light and the tight confinement of light can be seen as an enhancement caused by the increase of the optical density of states. The photons with a high density of states are also less influenced by the other photons with a lower density of states ($SPM < 2XPM$).

Another interesting effect lies in the magnitude of the shock term: the stronger this term is, the faster will the pulse form a shock front. The formation of an optical shock -more precisely the presence of the Kerr shock term- might play a very important role in the generation of dispersive linear waves (DSW) that could then be generated without strong dispersion requirements, for instance the presence of ZVD point is not mandatory [50]. The intensity of the XPM and its dependence with the angular frequency play then a predominant role. After the renormalization, the dispersive nonlinearity ($\gamma_1/\gamma_0$) in PhCWGs is only about one third of its nominal value in eq. (17); while at the same time the weight of the photons has been strongly modified. Consequently the behavior of PhCWGs with regards to this new DSW generation scenario is different to what is found in other systems. More generally, the renormalization redefines in a non-trivial way the interaction, mediated by the material nonlinearity, between the photons. Consequently, in the laboratory reference frame (i.e. considering only the power flux and no renormalization) the PhCWG behavior would be different to what is a priori expected.

### IV. IMPACT OF MODE AREA

Thus far, we only presented how to include the variations of the slow light factor in the renormalization process. The main reason is that the slow-light is responsible for most of nonlinearity variations in PhCWGs. However
the change in the effective modal area still accounts for about 25% of the dispersive nonlinearity; and we have shown through the comparison with an analytic test case that it has a noticeable impact on the overall nonlinear dynamics. Consequently we now show how the renormalization method is also able to include variations of the effective modal area.

By analogy with what has been done in the previous section, the renormalization procedure is applicable as well for the modal area subject that $gh_{\omega_0}$ can be decomposed as $gh_{\omega_0}(\omega_1, \omega_2, \omega_3) = g(\omega_0)h(\omega_1)h(\omega_2)h(\omega_3)$. If such is the case, then we define the renormalization function $m(\omega) = h(\omega)\sqrt{n_g(\omega)}$, and the propagation equation to solve becomes:

$$ \partial_z \Psi(z, \omega_0) = iK(\omega_0)\Psi(z, \omega_0) + i \frac{\omega_0 n_{2I}}{c A_{eff}(\omega_0)} \frac{n_g \omega_0}{n_0} g(\omega_0) h(\omega_0) \int \Psi(z, \omega_1)^* \Psi(z, \omega_2) \Psi(z, \omega_3) \delta(\omega_2 + \omega_3 - \omega_1 - \omega_0) d\omega_i^3 \quad (27) $$

Consequently $h(\omega)$ weights the individual photons' contribution to the Kerr nonlinearity: it includes implicitly most of the dispersive nonlinearity (in PhCWGs). On the contrary $g(\omega_0)h(\omega_0)/A_{eff}(\omega_0)$ stands for the dispersive part of the modal area. The problem of the factor decomposition of $gh_{\omega_0}$ is directly related to the question whether the nonlinearity in PhCWGs can be modelled accurately using an analytic function [20]. It has been demonstrated that, for dispersion engineered PhCWGs like the one considered in the present paper, the nonlinearity can be fitted by a Morse type potential function with four adjustable parameters. However the decomposition of such function in factor decomposition will only gives an approximate. More generally the best way to decompose the variations of the effective modal area is still an open question.

In any case, it is always possible to choose a decomposition that preserves the self-phase modulation ($SPM(\omega) \propto g(\omega)h(\omega)^3$) and cross-phase modulation ($XPM(\omega_{ref}, \omega) \propto g(\omega)h(\omega_{ref})^2h(\omega)$. $\omega_{ref}$ is defined as the center of the frequency domain.

$$ h(\omega) = \sqrt{\frac{SPM(\omega)}{XPM(\omega_{ref}, \omega)}} \quad (28) $$

$$ g(\omega) = \sqrt{\frac{XPM(\omega_{ref}, \omega)^3}{SPM(\omega)^3}} \quad (29) $$

If the decomposition is consistent -i.e. if $gh\omega_i$ can indeed be decomposed in factors-, $h(\omega)$ and $g(\omega)$ do not depend on the central frequency $\omega_{ref}$ that is chosen. Otherwise, only the XPM created by a pulse centered at $\omega_{ref}$, as well as the SPM for any frequency, are included correctly. Such approximation could still be sufficient if the propagation is dominated by a single strong pulse. Checking how $h(\omega_{ref})$ and $g(\omega_{ref})$ depends on $\omega_{ref}$ is essential to assess the validity of the renormalization for taking into account effects related to the effective modal area. In Fig 5-a), we show the value of $h(\omega)$ computed according to eq. (28) depending on the central frequency $\omega_{ref}$ (Y-axis).

<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) a) Color map: $h(\omega)$ computed for different $\omega_{ref}$ (y-axis). b) $h(\omega)$ averaged on the different $\omega_{ref}$ for the normal dispersion region (blue, $\lambda_{ref} > 1560$nm); and the anomalous dispersion region (red, $\lambda_{ref} < 1560$nm). Shaded area indicates the standard deviation of the average. c) Same as b), but showing the dispersive part of the nonlinearity $g(\omega)h(\omega)/A_{eff}(\omega)$.</figcaption></figure>

We observe two main zones: one ranging 1525 – 1560nm (anomalous dispersion) and the other one ranging 1560 – 1610nm (ranging from the first zero GVD to the band edge). The two zones are separated by the first zero group velocity dispersion (GVD) point (cf. Fig 2-a)). Inside each zone, $g(\omega)h(\omega)$ and $h(\omega)$ do not depend on the central reference frequency as seen in Fig 5-b,c). This brings forth two major conclusions. First it appears that inside each zone, the decomposition of $gh(\omega_i)$ holds; and it is therefore possible to describe correctly the pulse evolution using the renormalization technique (of course providing that the simulation spectral domain being confined within one zone). Second, the presence of two distinct and well defined zones indicates that there is actually a change in the physics governing the photon evolution.

For high frequencies (small wavelengths), the photons
have a nonlinear dispersive behavior corresponding to the fact that $h(\omega)g(\omega)/A_{eff}(\omega)$ varies while $h(\omega)$ is almost constant (Fig 5-b), blue). This indicates that the photons have a quasi equi-contribution to the nonlinearity. On the contrary, for low frequencies (long wavelengths), the individual contribution of photons is more pronounced (Fig 5-b), red) and $h(\omega)$ exhibits variations up to 60% while $h(\omega)g(\omega)$ remains flat. Weighting the eigen photons contribution to the nonlinearity is essential.

In dispersion-engineered PhCWGs, the slow light at long wavelength (low frequency) is caused by the presence of a complete Photonic Bandgap (PBG). Close to the band edge, the physics is then similar to what is found in Bragg Gratings [51]. In contrast, the slow-light obtained at higher frequency through dispersion-engineering has another nature and arises thanks to the complex interferences occurring inside the Bloch mode. Although the decomposition of $gh(\omega)$ into a product of functions is not mathematically exact, we see here that it could nevertheless be a useful metric to sort-out slow light [52] into categories depending on the equi-contribution to the nonlinearity that the photons have -or do not have-.

In the introduction, we disregarded the coupled set of equations as a possible solution for modeling systems with large variation of the nonlinear parameters. The fact that the photons are split here into two well defined domains tends to rehabilitate a posteriori such strategy. Still one must be careful on the way photons are taken into account at the limit between the two domains, especially considering that this point is precisely the zero group velocity dispersion (ZVD) point.

Finally we focused on SPM, XPM in the factor decomposition of $gh(\omega_i)$. The photon weight $h(\omega)$ depends on this peculiar choice. Therefore a material exhibiting a different nonlinear response, like the presence of a $\chi^{(2)}$ nonlinearity or a different form for the $\chi^{(3)}$ tensor, would lead to a different $h(\omega)$.

V. CONCLUSION

We have presented a new method to incorporate in an efficient way into the GNLSE the variations of the nonlinearity that exist in systems with structured slow-light: effects linked to variation of the slow-light enhancement factor could be taken into account through renormalization. Fluctuations of the mode effective area can be dealt with by this technique as well. This would be of importance especially for system like nanowires, where variations of the effective mode area are important. As a result this paper gives practical hints regarding the way nonlinear pulse propagation in nanostructured systems could be computed; and what could be the limitation of current models.

This is crucial as more studies are precisely focusing on higher order nonlinear effects and their mutual interplay [16, 18, 45, 48, 50]. Our model is consistent with an analytic set of equations derived to model discrete FWM events [21]; and could be considered as a generalization of that article. Besides, it is worth to note that our technique does not increase the computational burden compared to the resolution of a standard GNLSE.

This study was also the occasion for more fundamental considerations. In particular we found relevant not to compute directly the propagation of the power flux but to weight first the photon contribution by their optical density of states. Although the modal area depends on the considered nonlinear effect and its associated tensor, we found two classes of slow light in dispersion-engineered PhCWGs: engineered slow light with a normal nonlinear dispersive behavior, and a region close to the photonic band gap where the weight factor of the photons contributes greatly. Such study might change how we perceive and understand slow-light effects which are in fact more related to high density of states light effect.

Usually, the improvement of the models which deal with nonlinear pulse propagation comes along with the addition of extra operators in order to describe the new effects. Here the renormalized equation keeps the formulation with no extra terms added. Indeed, the new phenomena that are observed lie inside the renormalization function $m(\omega)$, not in the GNLSE itself. Namely the imbalance between SPM and XPM could be interpreted as inertial forces that appear because of the non-trivial relationship between the laboratory referential and the PhCWGs one, where some photons appear more immune to perturbation or prone to perturb others. Within this new reference frame, all the semi-analytic method/models developed so far, like the momentum method, remain valid [53–58].

Finally the present discussion only focused on PhCWGs, a system where the nonlinear variations are extreme. Part of our conclusions would anyway also apply to other nano-structured systems like nanowires where the relatively weaker nonlinear variations must be seen in regards to the very large optical bandwidth these systems support. Especially we have shown that variations of the effective area -which are dominant over slow-light in nanowires- can be taken into account by the renormalization method; and that this method can describe with accuracy both self-phase and cross phase modulation effects.

Acknowledgements

The author would like to thank A. Degiron for advice and direction, and acknowledges fruitful discussion with C. Husko, S Trillo and A De Rossi. This work has been supported through the Villum Fonden funded centre of excellence NATEC.
[1] J. I. Dadap, N. C. Panoiu, X. Chen, I.-W. Hsieh, X. Liu, C.-Y. Chou, E. Dulkeith, S. J. McNab, F. Xia, W. M. J. Green, L. Sekaric, Y. A. Vlasov, and R. M. Osgood, Opt. Express **16**, 1280 (2008).

[2] L. K. van Vugt, B. Zhang, B. Piccione, A. A. Spector, and R. Agarwal, Nano Letters **9**, 1684 (2009), pMID: 19265428, http://dx.doi.org/10.1021/nl900371r.

[3] N. Gutman, W. H. Dupree, Y. Sun, A. A. Sukhorukov, and C. M. de Sterke, Opt. Express **20**, 3519 (2012).

[4] D. Mori and T. Baba, Opt. Express **13**, 9398 (2005).

[5] J. Li, T. P. White, L. O'Faolain, A. Gomez-Iglesias, and T. F. Krauss, Opt. Express **16**, 6227 (2008).

[6] L. O'Faolain, S. A. Schulz, D. M. Beggs, T. P. White, M. Spasenović, L. Kuipers, F. Morichetti, A. Melloni, S. Mazoyer, J. P. Hugonin, P. Lalanne, and T. F. Krauss, Opt. Express **18**, 27627 (2010).

[7] M. Ebnali-Heidari, C. Monat, C. Grillet, and M. K. Moravvej-Farshi, Opt. Express **17**, 18340 (2009).

[8] P. Colman, S. Combrié, G. Lehoucq, and A. D. Rossi, Opt. Express **20**, 13108 (2012).

[9] F. Wang, J. S. Jensen, and O. Sigmund, J. Opt. Soc. Am. B **28**, 387 (2011).

[10] F. Wang, J. S. Jensen, J. Mørk, and O. Sigmund, J. Opt. Soc. Am. A **29**, 2657 (2012).

[11] J. Vuckovic, M. Loncar, H. Mabuchi, and A. Scherer, Phys. Rev. E **65**, 016608 (2001).

[12] E. Kuramochi, H. Taniyama, T. Tanabe, A. Shinya, and M. Notomi, Appl. Phys. Lett. **93**, 111112 (2008).

[13] Y. Tanaka, T. Asano, and S. Noda, Journal of Lightwave Technology **26**, 1532 (2008).

[14] L. Yin, Q. Lin, and G. P. Agrawal, Opt. Lett. **32**, 391 (2007).

[15] M. A. Foster, A. C. Turner, J. E. Sharping, B. S. Schmidt, M. Lipson, and A. Gaeta, Nature Letters **441**, 960 (2006).

[16] P. Colman, S. Combri, G. Lehoucq, A. de Rossi, and S. Trillo, Physical Review Letters **109** (2012), 10.1103/PhysRevLett.109.093901.

[17] D. I. Yeom, Opt. Lett. **33**, 660 (2008).

[18] C. Colman, P.; Husko, Submitted (2014).

[19] N. A. R. Bhat and J. E. Sipe, Phys. Rev. E **64**, 056604 (2001).

[20] S. Roy, M. Santagiustina, P. Colman, S. Combrie, and A. De Rossi, Photonics Journal, IEEE **4**, 224 (2012).

[21] M. Santagiustina, C. G. Someda, G. Vadala', S. Combrie', and A. D. Rossi, Opt. Express **18**, 21024 (2010).

[22] D. Li and C. D. Sarris, J. Lightwave Technol. **29**, 1003 (2011).

[23] A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method, 2nd ed.*, edited by M. A. H. Norwood (Norwood, MA: Artech House, 2000).

[24] J. Li, L. O'Faolain, I. H. Rey, and T. F. Krauss, Opt. Express **19**, 4458 (2011).

[25] A. Baron, A. Ryasnyanskiy, N. Dubreuil, P. Delaye, Q. V. Tran, S. Combrié, A. de Rossi, R. Frey, and G. Roosen, Opt. Express **17**, 552 (2009).

[26] C. Husko, S. Combrié, Q. V. Tran, F. Raineri, C. W. Wong, and A. D. Rossi, Opt. Express **17**, 22442 (2009).

[27] C. Monat, B. Corcoran, M. Ebnali-Heidari, C. Grillet, B. J. Eggleton, T. P. White, L. O'Faolain, and T. F. Krauss, Opt. Express **17**, 2944 (2009).

[28] D. Anderson and M. Lisak, Phys. Rev. A **27**, 1393 (1983).

[29] M. Trippenbach and Y. B. Band, Phys. Rev. A **57**, 4791 (1998).

[30] G. P. Agrawal, *Nonlinear Fiber Optics*, edited by A. Press (AP San Diego, 2007).

[31] J. R. de Oliveira, M. A. de Moura, J. M. Hickmann, and A. S. L. Gomes, J. Opt. Soc. Am. B **9**, 2025 (1992).

[32] T. X. Tran and F. Biancalana, Opt. Express **17**, 17934 (2009).

[33] X.-M. Liu, Phys. Rev. A **77**, 043818 (2008).

[34] T. F. Krauss, J. Phys. D: Appl. Phys. **40**, 2666 (2007).

[35] W. Ding, A. V. Gorbach, W. J. Wadsworth, J. C. Knight, D. V. Skryabin, M. J. Strain, M. Sorel, and R. M. D. L. Rue, Opt. Express **18**, 26625 (2010).

[36] C. Koos, L. Jacome, C. Poulton, J. Leuthold, and W. Freude, Opt. Express **15**, 5976 (2007).

[37] J. Sipe, J. Opt. A: Pure Appl. Opt. **11**, 114006 (2009).

[38] S. Hughes, Opt. Lett. **29**, 2659 (2004).

[39] V. S. C. MangaRao and S. Hughes, Phys. Rev. B **75**, 205437 (2007).

[40] M. Kolesik, J. V. Moloney, and M. Mlejnek, Phys. Rev. Lett. **89**, 283902 (2002).

[41] M. Kolesik and J. V. Moloney, Phys. Rev. E **70**, 036604 (2004).

[42] R. W. Boyd, *Nonlinear Optics* (Academic Press, 2003).

[43] P. Colman, C. Husko, S. Combrie, I. Sagnes, C. W. Wong, and A. D. Rossi, Nature Photonics **4**, 862 (2010).

[44] A. Blanco-Redondo, C. Husko, D. Eades, Y. Zhang, J. Li, T. Krauss, and B. Eggleton, Nat Commun **5**, (2014).

[45] C. A. Husko, S. Combrie, P. Colman, J. Zheng, A. D. R. Rossi, and C. W. Wong, Science Report **3**, 1100 (2013).

[46] F. Biancalana, D. V. Skryabin, and A. V. Yulin, Phys. Rev. E **70**, 016615 (2004).

[47] M. F. Saleh, W. Chang, P. Hölzer, A. Nazarkin, J. C. Travers, N. Y. Joly, P. S. J. Russell, and F. Biancalana, Phys. Rev. Lett. **107**, 203902 (2011).

[48] F. Raineri, T. J. Karle, V. Roppo, P. Monnier, and R. Raj, Phys. Rev. A **87**, 041802 (2013).

[49] M. Erkintalo, Y. Q. Xu, S. G. Murdoch, J. M. Dudley, and G. Genty, Phys. Rev. Lett. **109**, 223904 (2012).

[50] T. Roger, M. F. Saleh, S. Roy, F. Biancalana, C. Li, and D. Faccio, Phys. Rev. A **88**, 051801 (2013).

[51] B. J. Eggleton, C. M. de Sterke, and R. E. Slusher, J. Opt. Soc. Am. B **16**, 587 (1999).

[52] R. W. Boyd, J. Opt. Soc. Am. B **28**, A38 (2011).

[53] D. Anderson, Phys. Rev. A **27**, 3135 (1983).

[54] P. K. A. Wai, H. H. Chen, and Y. C. Lee, Phys. Rev. A **41**, 426 (1990).

[55] M. Santagiustina, J. Opt. Soc. Am. B **14**, 1484 (1997).

[56] V. Grigoryan, C. Menyuk, and R.-M. Mu, Lightwave Technology **17**, 1347 (1999).

[57] C. J. McKinstrie, J. Santhanam, and G. P. Agrawal, J. Opt. Soc. Am. B **19**, 640 (2002).

[58] S. Lefrancois, C. Husko, A. Blanco-Redondo, and B. J. Eggleton, J. Opt. Soc. Am. B **32**, 218 (2015).