An accurate envelope equation for light propagation in photonic nanowires: new nonlinear effects

We derive a new unidirectional evolution equation for photonic nanowires made of silica. Contrary to previous approaches, our formulation simultaneously takes into account both the vector nature of the electromagnetic field and the full variations of the effective modal profiles with wavelength. This leads to the discovery of new, previously unexplored nonlinear effects which have the potential to affect soliton propagation considerably. We specialize our theoretical considerations to the case of perfectly circular silica strands in air, and we support our analysis with detailed numerical simulations.


Introduction
Evolution equations that are able to accurately describe the propagation of ultra-short pulses in strongly nonlinear media are extremely valuable tools in modern computational nonlinear optics, and especially in fiber optics [1]. On one hand, it is always possible to directly solve Maxwell's equations in presence of linear dispersion and nonlinear polarization [2], capturing the full complexity of the nonlinear propagation, but in the majority of cases very large computational resources are needed for this, which means that one can only effectively simulate short propagation distances along an optical fiber. On the other hand, it is possible to drastically reduce the complexity of Maxwell's equations by using the envelopes of the electromagnetic (EM) fields, while still maintaining an accurate description of the important dynamical degrees of freedom. This can lead to different approximation schemes based on various different versions of the nonlinear Schrödinger equation (NLSE), that have been used quite successfully in the past to describe a variety of linear and nonlinear effects in large-core optical fibers. Amongst them, we can mention phenomena such as soliton propagation [3], modulational instabilities and four-wave mixing [4], and third and higher-order dispersive effects [1]. A specific extended version of the NLSE (also called the generalized NLSE, or GNLSE for brevity), which also includes terms describing the Raman effect [5], self-steepening [1] and the full complexity of the group velocity dispersion (GVD) of the waveguide, has been used in recent times to describe the important phenomenon of supercontinuum generation (SCG), i.e. the generation of a broad coherent spectrum from a short input pulse, with spectacular success [6]. This equation has led to a number of important advances in the theoretical understanding of SCG in ultra-small core fibers such as tapered fibers (TFs) [7] and photonic crystal fibers (PCFs) [8,9]. In particular the phenomena of emission of solitonic dispersive radiation (also called resonant radiation or Cherenkov radiation [10]), and the stabilization of the Raman self-frequency shift (RSFS) of solitons by means of such radiation [11,12], are recent important advances that have been possible only by means of a thorough theoretical understanding of the structure of GNLSE.
Equations like the GNLSE are based on several assumptions and approximations. The most important approximation used in all NLSE-based evolution equations is the slowly-varying envelope approximation (SVEA) in both space and time [13]. This amounts to decoupling the fast oscillating spatial and temporal frequencies from the envelope function, which is the only complex function used to describe the dynamics. This reduction is only possible if the envelope profile contains many oscillations, which is typically the case in the visible and mid infrared regime when the pulse duration is larger than a few tens of femtoseconds. SVEA allows to reduce the second-order wave equation for the electric field into a first-order equation (which only contains the forwardpropagating modes and neglects the backward-propagating ones) in which the evolution coordinate is z, the longitudinal spatial coordinate along the fiber. In this way, the computational complexity of numerical simulations is drastically reduced. The most crucial assumption on which the GNLSE is based is that the z-component of the electric field is very small with respect to its transverse components (weak guidance regime). This allows to approximately decouple the transverse dynamics from the longitudinal dynamics. The weak guidance regime represents a very good assumption for large-core fibers (when the core size is much larger than the wavelength of light) and for fibers with a low contrast between the refractive indices of the core and the cladding.
Recently, however, new types of waveguides with a sub-wavelength size of the core (broadly referred as photonic nanowires) have become accessible to fabrication, which both possess a complex cladding structure, that allows one to support solid cores with a sub-wavelength diameter, and have a strong contrast between the refractive indices of the core and the cladding. To this class of waveguides belong some specific examples of PCFs, such as extruded PCFs [14]; the extremely small interstitial features of Kagome hollow-core PCFs [15]; and TFs, i.e. silica rods with sub-micron diameters surrounded by air or vacuum [7]. As mentioned, the common feature of all these structures is the nonnegligible nature of the longitudinal component of the electric field (strong guidance regime) [16]. The importance of such component in the dynamics of light has not been clearly recognized in the past, until very recently [16,17]. In these latter works the authors demonstrate that the calculation of the nonlinear coefficient of photonic nanowires performed in the scalar theory [18,19,20] underestimates its real magnitude of approximately a factor 2, while the correct result can be found by using the more complete vector theory [16]. Such interesting conclusion is also supported by recent experiments performed on ultra-small core, highly nonlinear tellurite PCFs [21].
Inspired by these recent works, in the present paper we wish to further advance the theoretical understanding of the evolution equations in sub-wavelength core structures. Here we derive a new, logically self-consistent forward-evolution equation for photonic nanowires, based on SVEA, that simultaneously takes into account (i) the vector nature of the EM field (polarization and z-component included), (ii) the strong dispersion of the nonlinearity inside the spectral body of the pulse (i.e. the variation of the effective modal area as a function of frequency), and (iii) the full variations of the mode profiles with frequency. In all relevant previous works on the subject [13,16,17,22,23,24], all these properties were not included in the derivation at the same time. Considering (i), (ii) and (iii) altogether and with a minimal amount of assumptions, leads to new qualitative results that have not been possible to investigate with previous formulations. In particular, we demonstrate through a series of accurate simulations that (iii) leads to the emergence of new nonlinear terms that are also indirectly responsible for appreciable modifications in the dynamics of RSFS, in particular a suppression of the latter. Thus, even though the vector nature of the EM field in photonic nanowires may lead to an enhancement of the nonlinear coefficient in nanowires [16], our results show that this is actually counterbalanced by terms and effects that partially suppress such enhancement. However, our findings are not in contradiction with previous works on the subject, for reasons that we shall explain in the text.
The structure of the paper is the following. In section 2 we detail the derivation of our master equation, Eq. (19). In section 3 we describe the symmetries of the master equation, and we demonstrate that our model coincides with previously published models under various limiting cases. In section 4 we approximate Eq. (19) for long pulses, and find that new nonlinear terms arise, which were not considered in any previously published formulation. Such terms become increasingly important when decreasing the size of the core, and for a given core diameter their impact is maximized around a certain optimal wavelength. Finally, in section 5 we show results of numerical simulations, which underline the difference between our new theory and other known models. Conclusions are given in section 6.

Derivation of master equation
In this section we present a detailed derivation of our master equation, Eq. (19). Our starting point are Maxwell's equations written for a non-magnetic medium in the Lorentz-Heaviside units system [25]: where E(r, t) and H(r, t) are the electric and magnetic field respectively, t is time, r is the coordinate vector, c is the speed of light in vacuum and P NL (r, t) is the nonlinear polarization field that we shall specify later. The displacement field is given by D(r, t) = ǫ ⊗ E, where ǫ(r ⊥ , t) is the dielectric function that models the linear response of the medium, i.e. its dispersion. Its Fourier transform ǫ(r ⊥ , ω) gives the square of the refractive index along the transverse direction r ⊥ . The symbol ⊗ indicates a convolution: In the frequency domain, one has D(r, ω) = ǫ(r ⊥ , ω)E(r, ω).
In absence of nonlinear polarization (P NL = 0), the fiber has linear eigenmodes in the form of plane waves, solutions of the linearized Maxwell equations at a specified frequency ω [26]. These are given by: and where m is a collective index that indicates both the polarization state and the mode number, and β m (ω) is the propagation constant. N mω is a mode-and frequency-dependent normalization constant that will be defined shortly. Functions e m,ω (r ⊥ ) and h m,ω (r ⊥ ) give the profile of the fiber eigenmodes [26]. Note that e, h have all three components, including a z-component which becomes comparable with the other two transverse components for a sufficiently small core and for a relatively large contrast between the core and the cladding refractive indices, as pointed out in Ref. [16]. respectively, for a pump wavelength λ = 1 µm. Note the gradual expulsion of the total (transverse and longitudinal) field from the core towards the low-index medium (air in this case). For PCFs and for fibers with a highly non-trivial cladding structure, the vector profiles of the electric and magnetic eigenmodes must be calculated numerically by means of a mode solver, but for the specific case of optical fibers with a perfectly circular core, e and h are known analytically, provided that the propagation constant is known, see Ref. [26]. The full fields are expanded in terms of the eigenmodes as follows: where we have used the convention f (t) = [2π] −1 dωf (ω)e −iωt for the direct Fourier transform and f (ω) = dtf (t)e iωt for its inverse. It is customary to define the normalized fieldê mω ≡ e mω / √ N mω (and similarly for the magnetic field), in order to have an orthonormal basis on which to expand any physical field, analogously to what is done in quantum mechanics. We further note that in Ref. [16] the fields are expanded by assuming that the modesê mω andĥ mω are unchanged around a reference frequency ω 0 , which is arbitrary and taken to be equal to the central pulse frequency. Thus in Ref. [16] it is assumed that all spectral components inside the pulse will have the same mode profile, and there is no integration over the frequencies analogous to Eqs. (6-7), which is instead considered in Refs. [23,24,22]. In the present work we wish to be able to describe the full variations of the mode profiles with frequency because, as it was already demonstrated in Ref. [27], this actually leads to a large dispersion of the nonlinearity that strongly affects the dynamics of the resulting equations, and, most importantly, new nonlinear effects that we shall analyze in detain in section 4.
The correct orthogonality relations between the linear eigenmodes are given by the following cross product integral [26]: However, what is ultimately important is not the absolute normalization, but the relative change of the normalization constant (and thus of the modal effective area) with wavelength. That is why we decide to use in the following, without loss of generality, the orthonormal quantitiesê,ĥ, which satisfy the simpler conditions It is customary [16,23] to define a generalized Poynting vector as follows: By using Maxwell's equations (1-3) and the vector theorem , one can find the divergence of S: After integrating both members of Eq. (11) in the transverse dimensions, performing a time integration as in [23], and by using the reciprocity theorem [26] one obtains the following forward-propagating equation (first derived by Kolesik et al. in 2004, see Refs. [23]): Equation (12) is the basic equation for the common Fourier component of the EM fields, on which any forward-propagating theory must be based. Note that Eq. (12) is written for a generic frequency ω, thus no arbitrary reference frequency has been selected up to this point. Attempts to directly solve Eq. (12) have been given, in various forms, in Refs. [23] and [24]. However, this approach, although exact in absence of backward-propagating modes in the waveguide, does not provide a transparent understanding of the physical processes involved, and it does not lead to any serious reduction of the computational costs with respect to solving directly Maxwell's equations.
To further progress from Eq. (12), one needs to write explicitly the nonlinear polarization field for silica in terms of the electric field components. This can be done in several ways, depending on which particular medium or model one is interested in. For instance, in Ref. [16] the authors consider the vector Kerr nonlinearity of silica in absence of Raman effect. In Ref. [17] the Raman effect has been included, but emphasis was given to the analysis of stimulated Raman scattering. For our purposes we choose the following form for the nonlinear polarization of silica (see also Ref. [24]): where χ xxxx (r ⊥ ) is the third-order susceptibility coefficient, which is the only approximately independent component of the susceptibility tensor that survives the symmetries [28], and is a function of the transverse coordinates. R(t) is the Raman-Kerr convolution function that describes both the instantaneous (Kerr) and non-instantaneous (Raman) nonlinearities in the medium. In the case of silica one has: where is the Raman response function set by the vibrations of silica molecules induced by the EM field, with τ 1 = 12.2 fs and τ 2 = 32 fs. Coefficient θ parameterizes the relative importance between Raman and Kerr effect in Eq. (14), and the experimental value for silica is θ ≃ 0.18 [1]. In Eq. (14), δ(t) is the Dirac function and Θ(t) is a Heaviside (step) function that is introduced to preserve causality. In addition, the response function is normalized in such a way that +∞ −∞ R(t)dt = 1. Other, more precise forms of the Raman response functions for silica are known [1]. In a more accurate formulation, the isotropic and the anisotropic parts of the nuclear response should be included [29]. However, the following discussions are largely independent of the precise form of R(t), and are valid in general for any functional form of the response function.
After substituting expression (13) into Eq. (12), we explicitly obtain: Only at this point we choose an arbitrary reference frequency ω 0 ; furthermore, in order to capture the variations in the mode profile around ω 0 , we expand the linear modesê in Taylor series: where ∆ω ≡ ω − ω 0 is the frequency detuning from the reference frequency ω 0 , and the quantity f (j) m,ω0 is given by the term inside the square brackets in the last member of Eq. (16), and is proportional to the j-th derivative of the mode profile. A Taylor expansion equivalent to our Eq. (16) was also proposed in Eqs. (28) and (37) of Ref. [24]. However, in this latter work the author (being mostly interested in large core, low-contrast PCFs) gets rid of the longitudinal components of the EM fields [see his Eq. (13)], analyzes the scalar case, and does not elaborate further on the physical meaning of the corrections induced by the higher-order terms. Expansion (16) is also not considered in Refs. [16,30]. In our derivation, the expansion (16) is a crucial step. In fact, by taking into account the modifications of the quasi-linear mode profiles, we are led to a correct identification of new nonlinear effects that, under certain circumstances, strongly compete with the RSFS. We shall describe the latter effect in section 4.
By inserting expression (16) into Eq. (15), and after integrating over the transverse coordinates of the fiber, we can single out the following constant nonlinear coefficients: We shall analyze this complicated but very important object and its large symmetries in the next section. Γ jhpv mkld [in the following, the upper indices jhpv shall always indicate the order of the derivatives in the Taylor expansion in (16), while the lower indices mkld shall always denote the polarization and mode states] contains information about the nonlinearity experienced by the different modes or polarizations in the fiber. Note that the structure of Γ jhpv mkld closely resembles the one of the nonlinear polarization P NL inside Eq. (15). In addition, all components of Γ jhpv mkld have the same physical dimensions. In particular, it is convenient to use the usual units for nonlinear coefficients of W −1 m −1 for the quantities k 0 Γ jhpv mkld . The final step in our derivation is to find an equation in the time domain for the envelope of the electric field. In order to do this, we decouple the fast oscillating phase by introducing the envelope function Q, therefore using SVEA: In this paper we shall deal with perfectly circular photonic nanowires, i.e. silica strands surrounded by air or vacuum. In order to simplify our treatment we assume that the waveguide is operated in a single-mode regime, so that the index m = 1, 2 only represents the polarization states of the fundamental mode, which in perfectly circular fibers are degenerate [β 1 (ω) = β 2 (ω)] and orthogonal in the sense of the cross product, Eq. (9). This assumption can of course be relaxed without problems in our formulation, but here we wish to focus on the novel properties of our equation. After some regrouping and a Fourier transformation over the new independent variable ∆ω, we finally obtain our master equation: where is the dispersion operator that encodes all the complexity of the fiber GVD around the reference frequency [23], and is a function that naturally contains the dynamics of the shock term [through the term (1 + ∆ω/ω 0 )]. The convoluted nonlinear fields in Eq. (19) are given by: that also depends on polarization state indices and derivative orders in the expansion of Eq. (16). Equation (19) is the central result of the present work. It differs from other formulations in that it has been derived by making no assumptions or approximations on the frequency dependence of the field profiles of the linear modes of the waveguide. Higher-order derivatives of the mode profiles embedded into Eq. (19) are responsible for the appearance of additional nonlinear terms in the equations. These new terms obviously become increasingly small with increasing orders of the derivatives, since for every derivation there is a factor 1/ω 0 . However, as we demonstrate in section 4, the new first-order term has a strength comparable with or sometimes even larger than the Raman and Kerr terms, and there are circumstances in which it must have observable consequences.

Symmetries of the master equation and reduction to known models
The nonlinear coefficient Γ jhpv mkld given in Eq. (17) possesses very large symmetries that we wish to analyze in this section, following similar arguments as in Ref. [30]. Some of these symmetries are intrinsic, i.e. valid for any fiber geometry, while some other are strongly geometry-dependent. As we shall see, the symmetries of Γ jhpv mkld can be effectively used to reduce enormously the complexity of Eq. (19). As a practical important example, here we specialize the description of such symmetry relations for the case of a silica nanowire with a perfectly circular core embedded in air. In this case, we consider only the fundamental mode of the fiber, which has two polarizations m = 1, 2 with the same propagation constant [β 1 (ω) = β 2 (ω)]. In the following we always assume that when index m is fixed, then n = 3 − m denotes the opposite polarization. As before, indices j, h, p, v ≥ 0 are used to denote the orders of the derivatives in the definition (17) of Γ, and they can take any positive integer value.
Directly from the definition (17), one immediately gets the relation Γ jhpv mkld = Γ jphv mlkd , so the two lower an upper central indices can be simultaneously interchanged without changing the value of the nonlinear coefficient. Other simple relations (that can be directly verified by the reader) are the following: Γ jhpv mklm = Γ vphj mlkm , Γ jhhh mmnn = Γ hjhh mmnn , Γ jhhh mmmm = Γ hjhh mmmm . All the above equalities are always valid for any geometry of the waveguide.
In the special case of perfectly circular fiber, from direct computation of Eq. (17) by using the explicit formulas for the two (even and odd) degenerate polarizations of the fundamental mode given in Ref. [26], one obtains the useful equalities Γ jhpv , Γ jhpv mkld = Γ vhpj dklm , and Γ jhhh mnnm = Γ hjhh mnnm . Most importantly, many of these coefficients vanish identically: Γ jhpv mnmm = Γ jhpv mmnm = Γ jhpv mmmn = Γ jhpv nmmm = Γ jhpv nmnn = Γ jhpv nnmn = Γ jhpv nnnm = Γ jhpv mnnn = 0. We now demonstrate that our propagation equation (19) can be reduced to the previously published model of Ref. [16] by using the above symmetries, and by taking into account only the lowest order in the modal profile expansion Eq. (16). For this purpose, let us write the zero-th order approximation (j = h = p = v = 0) of Eq. (19): We now explicit the summations in Eq. (22) for the nonlinear functions Φ 000 kld , obtaining: According to the symmetries reported above, the fourth, fifth, sixth and seventh term of Eq. (23) vanish identically. Moreover, the eighth and ninth term have an identical nonlinear coefficient. A careful comparison between all terms of Eq. (23) shows that, indeed, the structure of this equation is identical to the one reported in Eq. (45) of Ref. [16], in the limiting case of a single mode regime, and without the Raman effect. However, our definition of the envelope Q as given by Eq. (18) automatically incorporates into the equations the frequency variations of the normalization constant (which is proportional to the effective mode area, see Ref. [27]), something that has not been done in the approach of Ref. [16], where the constancy of N mω around ω 0 was assumed.

New nonlinear effects in photonic nanowires
The higher-order derivatives in the Taylor expansion performed in Eq. (16) are expected to produce additional nonlinear terms, previously unknown, that become increasingly important when the core diameter is reduced and/or the refractive index contrast between cladding and core is increased. Here we attempt to obtain a physical understanding of these new terms, at least at the first-order in the Taylor expansion. In the following we just want to isolate and understand as clearly as possible the new nonlinear terms that come out of Eq. (19). In order to do this, we just take one polarization state (m = 1) for simplicity, we expand Eq. (19) and collect all the zero-th and first-order terms coming from the products between the G, Γ and Φ functions, obtaining (24) By using the relation (i/ω 0 )∂ t Φ 000 = Φ 100 + Φ 010 − Φ 001 (that can be derived just by explicitly calculating the derivative of Φ 000 ), one can slightly simplify Eq. (24) into where we have defined the constants γ 0 ≡ k 0 Γ 0000 and γ 1 ≡ k 0 Γ 1000 . A plot of the ratio r ≡ Γ 1000 /Γ 0000 = γ 1 /γ 0 is given in Fig. 2(a) for silica strands of various diameters. Interestingly, for a given diameter, r is predominantly positive, and becomes negative only in regions of very short wavelengths, not easily accessible experimentally (typically below 200 nm). Also, r is significantly larger than unity for wavelengths λ > Cd, where C ≃ 1.6 for perfectly circular silica strands in air. In order to gain some insight into the physical meaning of all nonlinear terms in Eq. (25), we consider relatively long pulses (of width t 0 ≫ 100 fs) for which the envelope evolves slowly along the fiber. Under this assumption, one can approximate the convolutions with nonlinear combinations of the field Q and its derivatives [1]. In particular, one has where The above simple Lorentzian model for the Raman response function gives T R ≃ 1.5 fs and T ′ R ≃ −23 fs 2 for silica fibers. Substituting Eqs. (26)(27) into Eq. (25), and neglecting the second order derivatives and all terms of order of ω −2 0 or T R /ω 0 , after some algebra we obtain our final expression (in the first-order of the perturbation theory) for the approximate evolution equation for long pulses: (28) In Eq. (28), the term indicated by T0 represents the well-known shock operator [1]. Terms T1 and T2 in Eq. (28) correspond to the usual and easily recognizable Kerr effect and RSFS respectively [1,31]. The additional term T3 is connected to a nonlinear change of the group velocity, proportional to the field intensity |Q| 2 . Such nonlinear term arises quite unexpectedly, and it is strictly related to the way in which the mode profile changes with frequency. Slightly different frequencies have different linear mode profiles, which in turn will have slightly different group velocities. This change becomes nonlinear, because the time derivative inside φ 001 in Eq. (21) is embedded into the full nonlinear convolution of Eq. (21), thus giving rise to term T3 in Eq. (28). It is now easy to see that ratio r of Fig. 2(a) quantifies the relative importance between the new nonlinearity T3 and the Kerr nonlinearity T1 in Eq. (28). Note that r can be as large as 6 in silica strands inside the frequency window accessible experimentally.
The relative importance between the T3 and T2 terms is regulated by a dimensionless parameter α S ≡ 4γ 1 /(γ 0 ω 0 T R ), which is shown in Fig. 2(b) as a function of wavelength for silica strands of various diameters. Again, this function is positive above 200 nm. Note that α S can be as large as 15 inside the frequency window accessible experimentally.
The presence of the extra first-order term T3 has not been previously reported in the context of unidirectional pulse propagation equations, and it constitutes an important result of the present work. It can be demonstrated [32] that this term may lead to either a suppression or an enhancement of the RSFS, depending on whether the coefficient α S is respectively positive or negative. In addition, the discussions above demonstrate that there are broad regions of wavelengths in which T3 can be several times larger than both the Kerr and the Raman terms. Fig. 3(a) shows the calculated value of the new first-order nonlinear coefficient α 1 ≡ 4γ 1 /(ω 0 t 0 ) as a function of wavelength, for silica strands in air of various core sizes and for t 0 = 100 fs. It is interesting to see that the optimal wavelength for which |α 1 | is largest approximately corresponds to λ ≃ Cd, again with C ≃ 1.6. Moreover, it can be seen in Fig. 3(a) that such a linear law seems to be a universal behavior, once that geometry and core and cladding materials are specified. Fig. 3(b) shows the calculated GVD of silica strands in air of various core sizes. By comparing Fig. 3(b) with Fig. 3(a) one can see that the condition λ ≃ Cd is always satisfied in correspondence of normal GVD for this kind of circular geometry. In our extensive numerical simulations we have seen that the new first-order nonlinear terms make a substantial difference in the pulse propagation only when the input pulse is launched in the anomalous dispersion, i.e. in presence of soliton formation, where it can enter in competition with the RSFS of solitons, while in the normal dispersion regime the contribution of the term proportional to φ 001 on the dynamics is largely negligible. For this reason, it seems that silica circular nanowires are not optimal for evidencing the effects of the new nonlinear term T3 regulated by nonlinear coefficient α 1 . Possibly this is why such new effects have never previously been noticed experimentally, although techniques to fabricate photonic nanowires have already been available since a few years [18]. A detailed study of the optimization of the new nonlinear terms in sub-wavelength-core PCFs and high-refractive index glasses will be carried out in a separate publication [32].
As a final note, we would like to stress again that the expansion of Eq. (27) is valid only for relatively long pulses. For ultra-short pulses, such as t 0 < 50 fs, all the physical effects and terms that are clearly visible in the formulation of Eq. (27) will be naturally embedded into the full convolution contained in Φ 001 , see Eq. (25).

Numerical analysis and comparison with previous models
In this section we show results of the extensive numerical simulations that was performed by using Eq. (19). We shall compare the output spectra obtained by pulse propagation in different dispersive regimes, and for the situation in which Eq. (19) is truncated at the zero-th order, as opposed as when instead all terms of first-order are included. As our representative case, we take a silica strand of core diameter d = 0.6 µm in air, the dispersion of which is displayed in Fig. 3(b). For this diameter, the two zero GVD points are approximately located at λ 1 = 0.468 µm and λ 2 = 0.718 µm. The input pulse duration is always taken to have a temporal width equal to t 0 = 100 fs, with the shape of a hyperbolic secant, and with an input power P such that the 'number of solitons' parameter N ≡ [γ 0 P t 2 0 /|β 2 (λ 0 )|] 1/2 is always equal to N = 7, where β 2 (λ 0 ) is the GVD coefficient calculated at the pump wavelength λ 0 . Propagation distances are given in dimensionless units: | is the second-order dispersion length at the pump wavelength λ 0 [1]. Also, for simplicity we assume for the moment that only one polarization state is excited (m = 1), and we neglect the dynamics of the other polarization state. In the final part of this section we shall also show some result obtained when taking into account the full polarization structure of Eq. (19).  Fig. 4(a) shows the output spectrum of an input pulse (displayed in black solid line) propagating in the normal dispersion regime of the nanowire, when pumping at λ 0 = 1.2 µm, and after a dimensionless propagation distance of ξ = 1 (L D2 ≃ 6 mm). At this wavelength we have r = 1.351 and α S = 2.295, see Fig. 2(a,b). The blue solid line in Fig. 4(a) shows the result of the simulation of equation (19), when considering the expansion of the Taylor series of the mode profiles up to the first order. The green dashed line in Fig. 4(a) shows the results of the simulation of Eq. (19), when considering only the zero-th order in the expansion. It is possible to observe that in the case of normal dispersion there is indeed only a small difference between the green dashed line and the blue solid line, which means that the new first-order nonlinear terms contained in Eq. (19) give only a relatively small contribution in this regime, as anticipated before.
When launching the input pulse at λ 0 = 0.8 µm (L D2 ≃ 10 cm), still in the normal dispersion regime but closer to one of the zero GVD points, all three curves nearly coincide, see Fig. 4(b). In fact, the new nonlinear coefficient is smaller than in the previous case, r = 0.3973 and α S = 0.45. Thus we can conclude that in absence of soliton dynamics, the two models give qualitatively and quantitatively very similar results.
A completely different scenario is observed in presence of solitons, see Fig.  4(c). This figure shows the output spectrum of a pulse launched in the anomalous dispersion regime at λ 0 = 0.55 µm, for a propagation of ξ = 1 (L D2 ≃ 30 cm). At this wavelength, r = 0.1537 and α S = 0.1196. Fig. 4(c) shows that in the initial stages of the propagation, the stabilization of RSFS by means of resonant radiation (at around 800 nm in the figure) emitted near the zero GVD point [11], occurs earlier when the equations do not take into account the new first-order nonlinear term T3 in Eq. (28), as the strongest soliton red-shifts at a faster rate in this latter (less precise) model. It was mentioned in section 4 that the novel nonlinear terms associated to the convolution φ 001 in Eq. (19) can give rise either to a suppression or an enhancement of the RSFS of solitons, depending on whether the sign of coefficient α S is positive or negative, respectively. This is due to a peculiar interplay between the Raman effect and the nonlinear change of group velocity induced by the new terms in the solitonic regime [32]. Despite the smallness of coefficients r and α S at this wavelength, a quite strong suppression scenario is already visible in Fig. 4(c). In fact, one can see that in the full model that takes into account the additional first-order nonlinear terms (blue solid line) solitons are pushed towards longer wavelengths at a much slower rate than in the model without first-order nonlinear effects (green dashed line). If the pump wavelength is located in a frequency region where α S is negative, exactly the opposite scenario would be observed, and the blue solid line would be shifted towards longer wavelengths with respect to the green dashed line, which means that the RSFS rate is increased. However, such a region of positive α S is located in the UV in circular silica strands [see Fig.  2(b)] and it is not reachable experimentally. Thus, for the geometry considered in this paper, the suppression of RSFS is the only physically meaningful scenario. Due to the relative smallness of |α S | inside the anomalous dispersion regime of the circular silica waveguides under consideration, the corrections to RSFS induced by the new nonlinear terms are also not as large as they could be. However, as we shall show in detail in a future publication [32], by using properly designed PCFs and high-refractive index materials, one is able to bring the wavelength for which |α 1 | has a maximum [see Fig. 3(a)] inside the anomalous dispersion regime. In such a case, the new nonlinear terms will completely change the dynamics of solitons, leading to qualitatively new regimes of nonlinear propagation in which the RSFS of solitons can be completely stopped for many dispersion lengths, without the use of stabilization techniques based on solitonic resonant dispersive wave emission [11].
In the same conditions as in Fig. 4(c), but for a slightly longer propagation distance ξ = 1.4 [see Fig. 4(d)], the solid blue and the dashed green lines become qualitatively very similar, apart from the visible differences in the dynamics of generation of resonant radiation, around 800 nm in Fig. 4(d), the bandwidth of which is somewhat reduced in the model truncated at the first-order (blue solid line) than in the model truncated at the zeroth-order (green dashed line).
In all the above calculations we have simulated only one of the two polarization states, in order to clearly display the differences between several models and approximations. However, Eq. (19) has a vector nature and actually represents a pair of coupled equations for the two polarization fields Q 1 and Q 2 . In Fig.  5(a) the final spectrum of a linearly polarized input pulse (95% on the m = 1 axis and 5% along the m = 2 axis) is shown, simulated by using the full Eq. (19). All other parameters for the pulse and the fiber are the same as in Fig.  4(c). Fig. 5(b) shows the time domain picture which corresponds to Fig 5(a). As expected, the formation of several vector solitons of different amplitudes is observed, each of them subject to RSFS according to their individual intensities.

Conclusions and future work
In conclusion, we have formulated a new set of vector unidirectional evolution equations [Eq. (19)] for silica photonic nanowires, i.e. for waveguides with a core smaller that the wavelength of light, and with a large step in the refractive indices of the core and the cladding, for which the z-component of the electric field is not negligible. Our model differs from previous formulations in that we simultaneously take into account the vector nature of the EM field and the variations of the linear mode profiles with frequency. Such variations generate new higher-order nonlinear terms, given by the Taylor expansion in Eq. (19). An analysis of the long-pulse limit of Eq. (19) has been carried out, which makes the physical understanding of the new first-order nonlinear terms implicitly contained in Eq. (19) much more transparent. This analysis led to Eq. (28), in which a nonlinear group velocity term appears, the strength of which is regulated by coefficient α 1 , which for silica strands in air of given diameter d takes appreciable values in the normal dispersion for wavelengths λ > Cd, with C ≃ 1.6. Numerical simulations of pulse propagation in silica strands were provided. The difference between the output spectra of our model and previous formulations is significant, especially when the propagation occurs in the anomalous dispersion regime, where soliton dynamics is mostly affected by the new nonlinearities.
Future works include the use of ultrasmall solid-core PCFs and interstitial features in the cladding of hollow-core PCFs, as well as high-refractive index nanowires, for the optimization of the novel nonlinear terms. This optimization, which makes use of the great flexibility in the dispersion tailoring of the above fibers, will bring the wavelength range at which the new nonlinearities are most effective (i.e. where |α 1 | is around its maximum) inside the anomalous dispersion regime of the fiber. This will create a truly qualitatively new and unexplored nonlinear behavior in such waveguides. This work is supported by the German Max Planck Society for the Advancement of Science (MPG). The authors would like to acknowledge Prof. P. St. J. Russell for stimulating discussions, and S. Stark for helping with the numerical codes.