A full vectorial model for pulse propagation in emerging waveguides with subwavelength structures part I: Kerr nonlinearity

The propagation of pulses through waveguides with subwavelength features, inhomogeneous transverse structure, and high index contrast cannot be described accurately using existing models in the presence of nonlinear effects. Here we report the development of a generalised full vectorial model of nonlinear pulse propagation and demonstrate that, unlike the standard pulse propagation formulation, the z-component of guided modes plays a key role for these new structures, and results in generalised definitions of the nonlinear coefficient $\gamma,$ $A_{eff}$, and mode orthognality. While new definitions reduce to standard definitions in some limits, significant differences are predicted, including a factor of $\sim2$ higher value for $\gamma$, for emerging waveguides and microstructured fibers.


Introduction
Nonlinear optical processes in optical fibers and waveguides have attracted significant interest because of the unique environment that they provide for nonlinear interactions, including tight confinement (high intensity), long interaction lengths, and control of propagation constants (see [1][2][3] and references therein). Recent and rapid progress in design and manufacturing of complex structured microstructured fibers and planar waveguides with subwavelength features (including both subwavelength inclusions and voids) has further extended the opportunities for guided-wave nonlinear optics and nonlinear devices by enabling extreme nonlinearity to combine with tailorable chromatic dispersion [1][2][3][4].
The nonlinear optical phenomena that occur in waveguides are determined through two main factors; the linear and nonlinear properties of the constituent bulk materials, and the optical properties of the waveguide. Two recent advances, as indicated below, have provided great potential to accelerate the field of guided-wave nonlinear optics: 1) the design and fabrication of complex structured waveguides with high contrast linear refractive indices and inhomogeneous cross sections, especially through postprocessing techniques. 2) the design and fabrication of waveguides with subwavelength features have opened up extensive opportunities for tailoring the nonlinear processes in waveguides.
While the characterisation of the nonlinear properties of bulk materials is indeed a rich and established field, the possibility of using postprocessing techniques to fill or coat complex structured waveguides with highly nonlinear materials has provided extended flexibilities in tailoring the nonlinear effects in waveguides and hence has opened new horizons for applications of nonlinear optical phenomena. The propagation of guided modes in waveguides with inhomogeneous cross sections are affected by their structure both directly, through the linear part of the refractive index which determines the modal characteristics, and indirectly, since the waveguide modes, under propagation through the structure, experience different losses, nonlinearity etc. Examples of such structures include liquid-filled [5][6][7][8] , gas-filled [9][10][11][12][13] , silicon nanocrystalsfilled [14], atomic vapor-filled [15][16][17][18] , or surface-functionalised recently attracted significant interest, both in planar waveguides, [3,[20][21][22][23][24][25][26][27][28][29] and fibers [30][31][32][33]. For such waveguides, it has been shown that a high intensity layer forms at the low refractive index side of the interface between two dielectric media due to the discontinuity of the normal component of the electric field. The intensity enhancement is proportional to the refractive index contrast of the two media. For the case of subwavelength voids in dielectrics, the enhanced intensity region forms at the surface of the void in the air side and extends over the whole void region since negligible evanescent decay of the field can occur within a subwavelength void. This key characteristic of the subwavelength features, when deployed within high intensity regions within the mode cross-section, can be used to achieve arbitrary distributions of high-intensity regions within waveguides, see for example [33]. We refer to this new class of optical waveguide, with high index contrast, inhomogeneous and complex structure, or subwavelength features, as 'Emerging waveguides' throughout this document.
Despite the importance and growing interest, applications, and publications in the field of nonlinear processes in these emerging waveguides, linear and nonlinear pulse propagation models for these structures still mainly rely on the well-known scalar Helmholtz equation [1,25,27,28,[34][35][36][37][38][39][40][41][42]. This equation is based on the weak guidance approximation, and assumes that the waveguide cross-section is homogeneous structure. However, even a cursory inspection of the key characteristics of these emerging waveguides (i.e. inhomogeneous, high index contrast transverse structure incorporating subwavelength features) reveals that these waveguides operate far from the weak guidance regime. Indeed, one can argue that the emerging waveguides considered here exhibit strong guidance. For example, it is observed that no Helmholtz wave equation can be obtained for the emerging waveguides since ∇.D = 0 in Maxwell's equation does not result in ∇.E = 0 because of the inhomogeneous nature of susceptibility tensor ε(x, y) [4]. Also, it has been pointed out that for subwavelength structures, such as optical nanowires, pulse propagation based on scalar theory does not give a good approximation [37].
There have been some reports of new models of nonlinear pulse propagation that take into account the inhomogeneous and vectorial solutions of Maxwell's equations [4,[43][44][45][46][47]. However, Refs. [4,43,44,45,47] ignored the potential for coupling between different modes including different polarisations, which as we will show later, is a key characteristic of linear and nonlinear pulse propagation when the full vectorial solutions of Maxwell's equations are considered. In fact, we demonstrate here for the first time that there are parameter regimes for which this modal coupling makes a significant contribution to both the predicted nonlinear and dispersive effects. Also, in some reports, [43,44,46], no consideration is given to the possibility of using an inhomogeneous cross-section which, as mentioned above, is a key feature of emerging waveguides. In addition, one vital aspect of the vectorial formulation of nonlinear pulse propagation, which has not been considered to date to the best of our knowledge, is the impact of the longitudinal component of the modal fields [the component along the propagation direction, (z)] on the dispersion, nonlinear, and modal/polarisation coupling behaviour of a pulse propagating through an emerging waveguides. One reason why the z-component of the fields has not been considered before, especially in studies that report high nonlinearity in waveguides e.g., [4,48], could be the fact that for slot waveguides TE modes, for which the z-component of the electric field is zero, have higher nonlinearity than those of TM modes.
Here a general vectorially-based Nonlinear Schrödinger Equation (VNSE), is derived for pulse propagation through waveguides with complex transverse structure including inhomogeneous refractive index profiles and subwavelength features. We demonstrate that in the strong guidance regime, the propagating modes have significant components along the direction of propagation, which causes the propagating modes to be non-transverse. As a result, this formalism predicts that a range of new tempo-spatial effects should be observable within emerging waveguides, such as dispersion-induced depolarisation.
Based on this VNSE, we derive a new and generalized equation for the A e f f , the parameter that defines the effective mode area, and γ, the parameter commonly used to describe the effective nonlinearity of an optical fiber [42]. These new definitions take into account both an inhomogeneous refractive index profile and subwavelength features. We apply these definitions to nanowires, and show that in some regimes, the value of γ can be a factor of two higher than that obtained using the standard definition. We provide an analysis of the value of γ predicted by this new generalised model. The new model also predicts new coupling terms between different modes or polarisations of propagating modes, which are due to non-transverse nature of the modes. Some preliminary results of the concept of the new model was presented in Ref. [49]. Here, we present the extension of our theoretical model and detailed results. We develop the theory of vectorially-based Nonlinear Schrödinger Equation (VNSE) in Sections 2, in which we derive new definitions for γ, A e f f , and mode orthogonality. We apply the new model to step index cylindrical waveguides and analyse and compare the results of the new model with those of the standard model in section 3. Concluding remarks are given in Section 4.
Next we consider Eqs. (4) and (5) for two sets of fields; unperturbed fieldsẼ 0 (r, ω 0 ) and H 0 (r, ω 0 ), which represent the electromagnetic fields of narrowband pulses at frequency ω 0 for which the dispersion, loss, and nonlinearity terms are zero, and perturbed fieldsẼ(r, ω) and H(r, ω), representing electromagnetic fields of frequency ω associated with wideband pulses centred at ω 0 , where the dispersion, loss and nonlinearity terms are nonzero. Vectorial solutions of Maxwell's equation for the unperturbed fields results in a complete orthonormal set of forward, backward and radiation propagating modes (labeled μ) with propagating constants of β μ and forward modal fields of; [50] where e μ (x, y, ω) × h * ν (x, y, ω).ẑdA = N μ δ μν (8) The propagation constant β ν and modal field distributions, e ν (x, y, ω 0 ) and h ν (x, y, ω 0 ), of propagating modes of a waveguide, in general, can be obtained through various numerical methods including Finite Element methods [51], the Multipole Method [52,53], etc. By constructing a function F C defined as F C =Ẽ 0 ×H * +Ẽ * ×H 0 and using Eqs. (4) and (5) for perturbed and unperturbed fields, we find [50] ∇ Next we expand the perturbed fieldsẼ andH according to the orthonormal and complete modal set of forward, backward, and radiation modes of the unperturbed field as [50]: Here index −μ refers to backward propagating modes and both forward and backward modes are orthogonal to radiation modes. Here, we only consider unidirectional pulse propagation for which we neglect the back scattering of a forward propagating laser beam and the nonlinearity associated with it [43,54]. This is not strictly true, especially for nonlinear and coupling processes where counter-propagating fields exist in the fiber. It can be shown that the backscattered field affects the overall nonlinearity for the forward modes but we leave the full investigation of this effect to future publications. Therefore, we only consider the first term in Eqs. (10) and (11) for expanding perturbed fields and hence ignore the coupling between the unperturbed field with the backward and radiation modes of the perturbed field. This will be discussed further later in this section. Unlike other reports that consider the modal expansion only for the nonlinear term [4,43,45], we consider the modal expansion in Eqs. (10) and (11) for both dispersion and nonlinear effects. It should also be noted that the frequency dependence of perturbed fields is totally contained within the coefficientsã μ (z, ω). Assuming that the unperturbed fieldsẼ 0 (r, ω) andH 0 (r, ω) are one of the propagating modes (e.g., mode ν) of unperturbed case i.e.,Ẽ and using the reciprocal theorem [50] ∂ ∂ z results in Here, Eq. (15) is a general first order differential equation that describes the propagation of amplitudes of the coupling coefficient of the perturbed field based on unperturbed one. This equation is similar to those reported in Refs. [43,46], except that the dispersion terms in Eq. (15) include the coupling between different modes.
Although,ã ν s are the coefficients of the perturbed fields, Eq. (15) is exact in the sense that no perturbation has been considered for dispersion and nonlinearity and hence this equation can be applied to describe, in general, any nonlinear or dispersion-based processes in an optical waveguide. The first and the second terms on the right hand side of the equation represent the dispersion and nonlinearity, respectively. Next, we perform a Taylor series expansion for the dispersion term in Eq. (15), around ω 0 . Depending on the bandwidth of the pulse around ω 0 , higher orders in the Taylor series can be considered to achieve better approximation for dispersion terms. For some cases, such as supercontinuum generation where extra wideband pulses propagating along the waveguide, however, it may be more appropriate to work with Eq. (15) directly. We separate the sum over the modes in Eq. (15) into self and cross terms to find; Here, to avoid confusion, superscripts (1) and (n) correspond to the first and higher order dispersion of the propagating modes, respectively, and subscripts μ label the mode number. It is straightforward to show that β is the group velocity as given in Ref. [50]. One important aspect of Eq. (18), which has not been reported before, is the existence of cross dispersive terms β Equations (21) and (22) result in a new process when the μ and ν refer to the two polarisations of 1 and 2 of one mode. It is well known that waveguides with three or higher-fold symmetries are not birefringent [55] i.e., for these waveguides β 1 = β 2 . In this case, the phase terms in Eqs. (21) and (22) are equal to unity and hence do not average to a negligible value as in the non phase matched case. The cross dispersive terms, β (1) 12 and β (n) 12 in Eqs. (21) and (22) are basically modifications to the group velocity β (1) 1 and higher order dispersion terms β (n) 1 of the polarisation 1. They have non-zero values which, as it will be shown later, is mainly due to the fact that in the strong guidance regime the dot product of the two polarisations of one mode i.e., e 1 .e * 2 and h 1 .h * 2 are non-zero because of strong z-component of the fields which results in non-transversality of the modes. This key finding is discussed in more detail later in this section. The physical consequence of this is dispersion-induced depolarisation of the guided mode, i.e., a polarised guided mode depolarises even if the incident beam is initially coupled perfectly to one of the polarisation axes of the waveguide. For instance, assuming that the incident beam is perfectly launched along the polarisation 1, then it can be deduced from Eq. (18) that the amplitude of the field along the polarisation 2, i.e.,ã 2 inside the fibre grows through ∑ n (Δω) n n! β (n) 21ã 1 . Next, we develop the time domain equivalent of Eq. (18). We multiply both sides of Eq. (18) by e −i(ω−ω 0 )t , integrate with respect to ω, consider the following definitions where a ν (z,t) and P NL (r,t) are the inverse Fourier transforms ofã ν (z, ω) andP NL (r, ω), respectively, to find (23) is a general equation that describes the nonlinear pulse propagation in the time domain. The first two terms on the right hand side of this equation describe the dispersion of a pulse propagating through a waveguide. The last term includes all the nonlinear effects and considers a shock term of (1 + τ shock ∂ /∂t) which is responsible for self phase modulation and self steeping of the pulse. This term, in various forms, has been considered in many publications [34,38,40,[55][56][57][58][59][60][61].
A few points should be noted here about the shock term: 1) similar to [43], it is naturally derived through the (∂ /∂t)P NL in the Maxwell's equations without any approximation of a second order time derivative in the scalar wave equation (Helmholtz equation) which is usually used to describe the nonlinear pulse propagation [34,[55][56][57].
2) The whole frequency dependence of the perturbed field is inherently included through the expansions Eq. (10) and (11) and is contained completely withinã 2 coefficients. Thus there is no need to include the frequency dependence of the propagating modes in the shock term, through A e f f , as has been done in Refs. [34,38,42,56,57]. 3) There is no dispersive term associated with χ (3) , i.e., (∂ /∂ ω)χ (3) = 0 since, as we show in Appendix, assuming a delta function form for the response function of Kerr nonlinearity results in frequency independence of χ (3) . For the nonlinear term in Eq. (23), since χ (2) = 0 for isotropic medium such as glasses, we only consider the third order Kerr nonlinearity for which we approximate P NL (r,t) ≈ P (3) (r,t) and assume that the nonlinear response function can be expressed in terms of delta functions, see Appendix, and hence P (3) (r,t), for Kerr nonlinearity, can be written as [63] where χ (3) is a rank four tensor and | indicates tensorial multiplication. Other third order nonlinear effects such as Raman scattering, for which the response function is not an instantaneous function of time will be a subject of future publications. Also it is assumed that P NL (r,t) ≈ P (3) (r,t) is a small perturbation compared to the linear induced polarisataion P L = ε 0 χ (1) (−ω; ω)Ẽ(r, ω), and higher order nonlinear effects are negligible. This is usually a valid approximation at low intensity fields and typical optical glasses due to their relatively low nonlinear properties. However, for some materials such as semiconductor-doped glasses [63][64][65][66][67] and some organic materials such as paratoloune sulphonate (PTS) [65,68] optical processes based on higher order nonlinearity can be observed, due to their higher order nonlinear susceptibilities, at moderate pulse intensity. For such materials, higher order terms must be considered in the nonlinear polarisation fieldP NL (r, ω) = ∑ ∞ n=2P (n) (r, ω). The components of χ (3) depends on the class symmetry of the crystal. Silica glasses have isotropic crystal structure [42] and silicon crystal, which is usually used in waveguides, have m3m point-group symmetry [1]. For isotropic materials, it can be shown that among 81 elements of χ (3) i jkl (i, j, k, l = x, y, z) only 21 are nonzero, which depend on only three independent quantities [69] i.e.; χ where Considering Eqs. (25), Eq. (24) can be written as; where i and j refer to x, y, z. For Kerr nonlinearity with the choice of frequencies in Eq. (24), i.e., χ (3) (−ω 0 ; ω 0 , ω 0 , −ω 0 ), the condition of permutation symmetry requires that χ xyxy . The magnitude of the terms in the right hand side of Eq. (26) depends on the origin of the nonlinear term. In the case of silica and other glasses they are mainly nonresonant electronic origins for which χ [42,69] and hence Eq. (27) can be simplified to For silicon, however, the third order nonlinearity can be expressed based on four independent values as χ where χ d ≡ χ xxxx − χ xxyy − χ xyxy − χ xyyx represent the nonlinearity isotropy. Similar to Silica, for the choice of frequencies in the third order susceptibility and photon energieshω well above E g χ xxyy (−ω; ω, −ω, ω) = χ xyyx (−ω; ω, −ω, ω) ≈ χ xyxy (−ω; ω, −ω, ω) [1]. As a result, χ i jkl becomes χ where ρ ≡ 3χ xxyy /χ xxxx characterizes the nonlinear anisotropy and its value in the telecom band is real and close to 1.27 [1]. Using this, Eq. (24) can be written for silicon as . Within the rest of this paper we ignore the last term of Eq. (30), which in fact affects the polarization dependence of nonlinear phenomena inside silicon waveguides, and thus Eq. (30) is the same as Eq. (28) except the factor ρ. Considering the expansion in Eq. (10) and using Eq. (28) we can evaluate the integrand in the last term of Eq. (23), i.e., (1/ √ N ν )e −iβ ν z e * ν .P NL (r,t) as; where, Greek indices μ, ν, η, ζ represent different modes of the waveguide. The terms on the right hand side of this equation, once integrated over the waveguide cross section, are overlap integrals representing how different propagating modes of the fiber couple to each other through the nonlinearity. Equation (31) can be expanded as sum of terms with and without phase terms as: other phase terms}.
The first two terms on the right hand side of Eq. (32) are automatically phase matched while the rest of the terms require phase matching in order to make significant contributions. The phase terms are responsible for nonlinear-induced depolarisation or four-wave-mixing [42]. They can be phased match, depending on Δβ νμ = β ν − β μ . This can be achieved by employing the flexibility in controlling the dispersion properties of MOfs through structure design and glass choice. By substituting Eq. (32) into Eq. (23), a first order differential equation is obtained which describes the nonlinear pulse propagation in a multimode waveguide. An important aspect of our formalism is related to the orthogonality of the waveguide propagating modes. Contrary to the standard formalism [42,69,70], for which e ν s are approximated to be transverse modes and e * ν .e μ dA = e * νt .e μt dA = 0, or in the case of different polarizations e ν .e μ = 0, in our formalism e * μ .e ν dA = 0 (or e ν .e μ = 0 if μ and ν are the two polarizations of the same mode) because the modes are non-transverse, i.e., they have non-zero z-component. The generalized orthogonality condition, which is valid even in the strong guidance regime is ( e ν × h * μ ).ẑdA = δ νμ inherently includes the z−component of the fields. Considering that and [50] one can show that which considering the general orthogonality relation ( e ν × h * μ ).ẑdA = δ νμ results in where subscript t refers to transverse component of fields and operators. Eq. (36) clearly shows that e νt · e * μt dA = 0 in the parameter regime where z−component of electromagnetic fields are non zero. It should also be noted that Eq. (32) has been obtained by ignoring the backward and radiation terms in Eqs. (10) and (11). Considering these two terms in the expansion Eqs. (10) and (11) results in coupling between forward-backward and forward-radiation modes, which are represented by dot products of forward modes with backward and radiation modes in Eq. (32). These terms describe the power coupling between a forward propagating mode and backward and radiation modes due to nonlinearity. We leave the full investigation of these coupling effects to future publications.
In the case of single mode fibers, where two independent polarizations exist in the waveguide, μ and ν refer to the two polarizations 1 and 2 and Δβ νμ is the linear birefringence of the waveguide. For waveguides with strong birefringence, the beat length L B = 2π/Δβ is short, and hence for fiber lengths L >> L B the phase terms oscillate very fast and hence have negligible contributions. However, for waveguides with weak birefringence, for which L < L B the phase terms are not negligible and should be taken into account. In the following sections we develop a model for nonlinear pulse propagation in single mode waveguides for both weak and strong birefringence.

Single mode highly birefringent waveguides
In the case of single mode waveguides with high birefringence, where only the first two terms in Eq. (32) are significant, substituting Eq. (32) into Eq. (23), considering [1,42] where μ,ν = 1, 2 and μ = ν refer to the two polarisations of the fundamental mode. This equation can finally be written in a simple form; Eqs. (38) is the final form of nonlinear pulse propagation inside a single mode birefringent waveguide, which in form is similar to the commonly used equation (see [42]) but with the effective nonlinear coefficients of the waveguide are now given by the generalized forms as in Eq. (39) and (40). By generalizing the definition of the A e f f as the nonlinear coefficient γ ν can be rewritten as where n 2 can be viewed as nonlinear refractive index averaged over an inhomogeneous cross section weighted with respect to field distribution. The advantages of writing γ as in Eq. (42) over the other reported form [4] is that it allows the analysis of γ to be separated into parts describing linear (geometry and n(x, y) =⇒ A e f f ) and nonlinear (mode profile and n 2 (x, y) =⇒ n 2 ) characteristics, providing a more intuitive analysis of γ. In our formalism, A e f f has its standard interpretation as the effective area of the propagating modes, which can be determined purely based on the geometry and the linear refractive index of the waveguide n(x, y), and does not require to be considered as the "effective nonlinear interaction area" as in Ref. [4]. Considering Eqs. (33) and (35), we find that Eq. (41) can be written as This equation in the limit of small z-component of electric field, where ∇ t e z can be ignored in comparison with β e t , simplifies as which is the standard definition of A e f f [42]. Comparing Eqs. (38)-(42) with the expression commonly used for nonlinear birefringent terms [42], shows 1) the vectorial-based γ developed here (referred to hereafter as γ V ) in Eq. (42) includes both the inhomogeneous waveguide structure and vectorial nature of electromagnetic fields and 2) an extra term that induces the depolarisation of an initially polarised beam through a non-zero coupling term γ (1) μν ∝ n 2 (x, y)n 2 (x, y)( e ν .e * μ 2 + e ν .e μ 2 )dA = 0. The nonzero nature of this term, especially in the strong guidance regime as will be shown in Sec. 3, is the direct result of two facts; 1) the two different polarizations are not perpendicular to each other in the common sense, i.e., e μ .e ν = 0 or e μ .e * ν dA = 0 because of the strong z−component of the fields and 2) the transverse integral, due to transverse dependence of n 2 (x, y) and n 2 (x, y), can be evaluated over different regions with different n and n 2 .

Single mode non-birefringent waveguides
For waveguides with perfect three or higher fold symmetry, the fundamental mode of the waveguide is degenerate or non-birefringent [55], i.e., β 1 = β 2 where 1 and 2 refer to the two polarisations. Therefore, the phase factors in Eq. (32) are equal to 1 for the pair of fundamental modes of these waveguides and hence Eq. (32) and Eq. (23), result in : where μ,ν = 1, 2 and μ = ν and refer to the two polarisations of the fundamental mode. Eq. (45) is the vectorial generalisation of the common nonlinear pulse propagation for the two polarisations of single mode fiber [42]. The main differences are the extra contributions from the different combinations of e μ .e ν . These terms, as indicated in the previous section, have nonzero values in the regime of high index and subwavelength core diameters.

Results and Discussion
The formalism developed in the previous section is general, and can be applied to an arbitrary waveguide. However, here we apply the above formalism to a simple step-index rod waveguide (i.e. a nanowire) and demonstrate that its nonlinear behaviour is predicted to be significantly different in the regime of high index contrast and subwavelength features than the predictions made using the standard formalism [42]. It should be pointed out that, throughout this paper we use the unit of W −1 m −1 for γ instead of the commonly used unit W −1 km −1 [42]. This is justified considering the fact that small core structures in high index silicon or glasses can now provide access to waveguides with extremely large γ values [4,27,39,71]. A comparison between γ V developed here and the standard definition given by Agrawal [42]; where F(x, y) = e t (x, y) is the scalar transverse electric field, indicates that γ V accounts for both inhomogeneous waveguide cross section and full vectorial nature of the propagating modes of the waveguide, especially the z-component of the modes. Figure 1 shows the effective nonlinearity for the two definitions of γ as a function of core diameter for three step-index rods with different host materials of silica [42] ( n = 1.45, n 2 = 2.6 × 10 −20 m 2 /W ), bismuth [72] (n = 2.05, n 2 = 3.2 × 10 −19 m 2 /W ), and silicon [73] (n = 3.45, n 2 = 4.5 × 10 −18 m 2 /W ).
While the γ V based on VNSE approaches γ A in the limit of large core diameter, it is significantly higher for small core diameters. For example, in Fig. 1c γ V is a factor of 2 higher than the γ A for silicon at the core diameter of D = 0.19 μm. Figure 1 also indicates that the difference between the γ values increases as the index contrast of the core and cladding increases. Here, we have also considered another definition of γ, given by Foster et. al. [20], as; Foster et. al. [20] refer this equation to Agrawal [42]. It seems that Foster et. al. have just simply replaced |F(x, y)| 2 = |e t (x, y)| 2 in Eq. (47) with the S z = (e ν × h * ν ).ẑ , arguing that S z is the intensity of light propagating down the waveguide. However, this replacement is only valid for the transverse mode approximation where |F(x, y)| 2 = |e t (x, y)| 2 is proportional to (e ν × h * ν ).ẑ and interpreted as the intensity of the light, see Eq. (35). In general, for full vectorial formalism, the nonlinear-induced polarisation is expanded in terms of different powers of electric field strength, see theory section, which for the χ (3) nonlinearity has the form given in Eq. (30), which in turn results in the new definition of γ V as in Eq. (39).  In Fig. 1, we have also plotted the values of γ F as a function of core diameter for different glass materials. While for silica, with the lowest index, γ V and γ F curves are on the top of each other for silicon with highest index, there is a maximum difference of 40% between the two curves at D = 0.15 μm. Also there is a difference of about 14% between the maximum values of γ V and γ F .
The differences between γ V developed here and γ A and γ F are attributed to the fact that propagating modes of a waveguide are not transverse in strong guidance regime. In order to demonstrate this, we have defined the transversality [74] of a mode as T ν = 1 − e 2 νz dA/ e 2 ν dA and plotted it as a function of core diameter for different materials, as shown in Fig. 2. It indicates that in the regime of large core diameter, the transversality approaches 100%, i.e., the modes become essentially transverse as expected. However, in the limit of small cores, the transversality is reduced, indicating that a large fraction of electric field power is contained within the z-component of the field. This effect is more profound for the waveguides made from high index glass than the low index ones.
μν , is proportional to n 2 (x, y)n 2 (x, y)( e ν .e * μ 2 + e ν .e μ 2 )dA, which contributes to the overall nonlinearity of mode ν, and is due to the overlap of the two different modes (two different polarizations in the case of single mode waveguide). This term does not appear in formalisms [42] where fully-transverse propagating modes and homogenous cross section are assumed since either e ν .e * μ dA = 0, due to orthogonality of transverse modes or e ν .e μ = 0 if the modes are the two polarizations of a single fundamental mode.
In our formalism, however, this term appears because of the non-zero z-component of the fields which result in e ν .e * μ 2 dA = 0 for any two propagating modes of a waveguide or e ν .e μ = 0 even when the two modes are the polarizations of the fundamental mode of a single mode waveguide. Figure 4 shows the magnitude of γ (1) μν relative to γ ν , i.e., γ (1) μν /γ ν where μ and ν are the two polarizations of a simple step index rod, as a function of core diameter. It demonstrates that indeed in the limit of large core, the relative value of γ (1) μν approaches zero but for small core diameters its value is enhanced significantly. Comparing the ratio γ (1) μν /γ ν for three different materials silica (n = 1.45) and bismuth (n = 2.05) and silicon (n = 3.45) demonstrates that the value of γ (1) μν is more significant in subwavelength regime and for large index contrast host materials.
The behaviour of the nonlinear coefficient γ V as a function of wavelength also shows a sig-   (1) μν and γ ν , see Eq. (40) and (39), as a function of core diameter for step index rods with host materials Silica (n = 1.45), Bismuth (n = 2.05), and Silicon (n = 3.45). The cladding material for glasses is air and for silicon is silica. Signs and the solid lines are the calculated data and lines of best fit, respectively. nificant difference between the usual definitions γ A and the one based on our VNSE. Fig. 5a) shows the behaviour of γ as a function of core diameter for different wavelengths; λ = 532, 633, 800, 1064, 1310 , and 1550 nm. As it is also evident from Fig. 5b, the γ values decreases as the wavelength increases, but the decrease in the maximum value of γ is much faster for the γ V of our model compared to that of common definition γ A . It is also evident from Fig. 5a) that the position of maximum of γ shifts to larger core diameters as the wavelength increases. Very high values of γ at short wavelengths are due to the tighter confinement of the propagating mode of the waveguide and hence their higher intensities. The possibility of achieving γ values of 150 W −1 m −1 , see Fig. 5b, suggest that order of π nonlinear phase shift should be achievable for input powers of order of 20 mW and for effective fibre length of 1 m in the visible spectrum. The dispersion properties of the γ values can be better compared by examining Fig. 6 in which the wavelength behaviour of γ at constant core diameters are shown. While for large core diameter, e.g. D = 1.6 μm, the γ (λ = 1550 nm) increases by a factor 3.5 to γ (λ = 532 nm) and the two definitions of γ are very close, for small core diameter D = 0.5 μm the difference between the γ values at λ = 1550 nm and λ = 532 nm is the order of 20 times and a large difference between the two definitions of γ is observed. This indicates higher dispersion of γ at small core diameters than that of large core diameters.

Discussion and conclusion
A new frontier in the field of optical waveguides is the design and fabrication of waveguides, referred to here as "emerging waveguides" with three main characteristics; 1) complex and inhomogeneous structure, 2) high index contrast, and 3) subwavelength features. A direct consequence of these features is that the common nonlinear Schrödinger equation which is based on weak guidance approximation does not provide accurate description of nonlinear processes in these waveguides. Here, we developed a vectorial based Nonlinear Schrödinger equation (VNSE), without relying on weak guidance approximation, which can be applied to these emerging waveguides.
An important feature of these waveguides is the fact that their propagating modes have much bigger z-component (z; direction of propagation) in comparison with those waveguides for which the weak guidance approximation is valid. As a result, the modes of these waveguides are not fully transverse and hence different orthogonality conditions govern them. Nonlinear  and dispersion processes can be associated to this third "direction of propagation polarisation", which can also couples the transverse polarisations through dispersion and nonlinear processes. Our model provides a platform for generalising nonlinear processes such XPM, Modulation Instability, Soliton formation and propagation, Four Wave Mixing, Parametric Processes and Raman and Brillouin Scattering for emerging waveguides. The model also predicts new tempospatial processes such as dispersion-induced depolarisation of the guided modes. Despite the complexity that this new concept brings into guided-wave nonlinear optics, early redevelopment of Kerr and Raman processes indicates the great flexibility for "Engineering" nonlinear processes.
Based on the model developed here, we provided a new vectorially-based definition for the effective nonlinear coefficient of waveguides (γ V ). Although the developed model is general, we have applied it to a simple step index cylindrical waveguide and shown that even for such a simple structure, γ V can be a factor of two higher than the common definition of γ, in the regime of high index contrast and subwavelength dimensions. However the full extent of the model in terms of exploring the rich physics behind the new pulse propagation model and the new definitions of effective nonlinearity γ, effective mode area A e f f , and cross-mode effective nonlinearity γ μν , especially for waveguides with inhomogeneous structures, is yet to be explored.
The pulse propagation model developed here [see Eqs. (38) and (45)] adds some complexity in terms of numerical solutions in comparison with the standard model [42]. It includes calculations of different overlap integrals of the propagating fields and the linear and nonlinear index distribution, see Eqs. (38) and (45). These integrals, however, are numerically easy to take and need to be evaluated only once, for any given fibre, before numerically solving the pulse propagation equations. The model also implies that coupled pulse propagation equations of different modes, either different polarisations of the mode of a single mode waveguide or different modes of a multimode waveguide, must be solved to give an accurate picture of Kerr nonlinear process for a propagating pulse, especially in the parameter regime for which our formalism gives very distinct result in comparison with the standard model, i.e, high index contrast, inhomogeneous structure, and subwavelength features.
Experimental measurement of the effective nonlinearity in the parameter regime where there is a distinct difference between the new model and the standard model, will be crucial to confirm the results of new model. It should also pointed out that within this paper, Kerr nonlinearity has been considered for waveguides with inhomogeneous structure assuming instantaneous nonlinear response (i.e., response time much shorter than the pulse width τ R << τ P ) of the materials everywhere in the waveguide. The case of finite and inhomogeneous nonlinear response time will be considered in future publications.
In general the third order induced polarisation is related to electric fields, in time and frequency domains, through the following equations [63]; where R (3) is the third order response function. Assuming that for Kerr nonlinearity the response function is product of delta functions [42] i.e., with S (3) being a constant, we find using Eq. (49);