Modeling of strain-induced Pockels effect in Silicon

We propose a theoretical model to describe the strain-induced linear electro-optic (Pockels) effect in centro-symmetric crystals. The general formulation is presented and the specific case of the strained silicon is investigated in detail because of its attractive properties for integrated optics. The outcome of this analysis is a linear relation between the second order susceptibility tensor and the strain gradient tensor, depending generically on fifteen coefficients. The proposed model greatly simplifies the description of the electro-optic effect in strained silicon waveguides, providing a powerful and effective tool for design and optimization of optical devices.


Introduction
In recent years, strain engineering is emerging as a new frontier in micro and nano-technology. By varying the elastic strain it is possible to turn on physical and chemical properties that are absent in the unstrained material. As a result, electronic, optical, magnetic, phononic and catalytic properties of a material can be tuned by compressive or tensile stress [1,2].
In optics, tensile strained germanium and strained silicon are attracting a great deal of interest. Tensile strained germanium-on-silicon can be used as active material for the shortwavelength infrared light and it can be an efficient solution for manufacturing monolithic lasers and optical amplifiers [3,4]. On the other hand, Pockels effect has been experimentally meas-ured in strained silicon [5], making it a promising candidate material for realizing very fast integrated optical modulators and switches. The modeling of strain-induced electro-optic Pockels effect in silicon is the main object of this work.
The electro-optic effect consists in the change of the refractive index induced by an electric field that varies slowly compared with the frequency of an optical signal [6,7]. In the particular case of the Pockels effect, also called linear electro-optic effect, the change in the refractive index is proportional to the applied electric field, providing an efficient physical mechanism for optical modulation. However, a peculiarity of the Pockels effect is that it arises only in crystalline solids lacking of inversion symmetry [8]. As a consequence, for centro-symmetric crystals (like silicon) the Pockels effect can be observed only when the inversion symmetry is broken, e.g., by the presence of significant surface/interface effects [9][10][11] or by an inhomogeneous mechanical stress.
Since the seminal work [5] there have been considerable progresses in the fabrication of electro-optic modulators based on strained silicon. In 2011 the first fully integrated Mach-Zehnder interferometer (MZI), based on strained silicon rib waveguides, was manufactured and the value of 122 pm/V for the effective electro-optic susceptibility was measured [12]. Two years later, the same authors presented a detailed investigation of the local strain distribution and of the induced optical nonlinearity as a function of the waveguide width, measuring the record value of 190 pm/V for a 300 nm large rib waveguide [13]. The dependence of the second order dielectric susceptibility on both wavelength and waveguide width was later investigated for a channel waveguide, showing that higher values of the effective susceptibility can be reached for narrow waveguides and large wavelength [14].
The more recent analysis performed in [15] and [16]suggests that the phase shift observed in previous MZI experiments (and thus the corresponding effective index variation) can be also related to the presence of free carrier variation inside the waveguide. The inversion of the phase shift observed when switching the applied tension, which was previously interpreted as the smoking-gun signal for the Pockels effects, is now attributed also to the surface charge present in the silicon nitride cladding. The results presented in [17], where the free-carrier plasma dispersion effect in silicon waveguides has been theoretically characterized, further support the possibility that the results obtained in [12,13] and [14] suffer from contamination of free-carrier contribution.
The most natural way to avoid contamination from free-carriers would be to perform highspeed measurements, with temporal resolution smaller than the free-carrier response times. In particular, an unambiguous indication for the presence of a strain-induced second order susceptibility is given by second harmonic generation (SHG) measurements. The experimental value of the second order dielectric susceptibility extracted from SHG measurements is of 40 pm/V [18], which is expected to be of the same order of magnitude of the one related to the Pockels effects. The comparison of these two effects can however be only qualitative, since they are associated with different frequency components of the nonlinear susceptibility.
On the one hand, when the free-carrier concentration are known, the corresponding effective index variation can be predicted by the empirical formula proposed by Soref and Bennett in 1987 [19], and latterly improved by Nedeljkovic et al. [20]. On the other hand several approaches have been proposed to model the strain-induced Pockels effect in silicon, however no one can be effectively used for practical purpose.
In [21] it is shown that a simplified classical model of a 2D centro-symmetric lattice is not able to reproduce the correct order of magnitude of the experimentally measured susceptibilities, moreover, in this approach, the way in which the strain enters the computations does not appear to be completely justified. A more sophisticated model is used in [18], where the linear electro-optic effect in strained silicon is studied by using the time-dependent density-functional theory [22,23]. While this ab initio method is theoretically well founded, it has the obvious drawback of being computationally very expensive.
Some intuitive ideas already present in the literature relate the effective susceptibility to the strain gradient (see, e.g., [13,14,18,24]). Indeed, as explicitly noted in [14]: "Despite the lack of general proof for this claim available in the literature yet, it has been widely accepted that the second order nonlinear effects in strained silicon are caused by the variations of strain i.e. strain gradients inside the crystal." In this work we will show that such a relation can be deduced by using just symmetry arguments and the specific case of the strained silicon will be investigated as a particularly interesting example. The final result will be a simple linear relation between the second order effective susceptibility tensor and the strain gradient tensor (weighted by the electromagnetic modes), depending generically on fifteen independent coefficients that can in principle be obtained from experimental measurements. Once these coefficients are known, the computation of the electro-optic effect is reduced to a standard strain computation and electromagnetic mode analysis, thus providing an easy framework for the optimization of optical devices.

Nonlinear susceptibility and Pockels effect
In this section we recall, for the benefit of the reader and to fix the notation, the basic properties of the quadratic nonlinear susceptibility that will be needed in the following.
The quadratic nonlinear susceptibility is conventionally defined, for a local and causal medium, by the following relation between the (quadratic component of the) polarization vector P P P (2) (t) and the electric field E E E(t): where ε 0 is the free-space dielectric permittivity, χ i jk is the susceptibility tensor. Sum over repeated indices is always understood in Eq. (1) and followings unless otherwise explicitly stated. The dependence on the position x x x is omitted for the sake of the simplicity but will be important in the following sections. The previous relation can be rewritten in the frequency domain as following where the first argument of the susceptibility is just the sum of the two other frequencies, a notation conventionally adopted in the literature. The symmetry properties of the tensor χ (2) i jk (ω 1 + ω 2 ; ω 1 , ω 2 ) will be particularly important in our analysis. From the definition and the reality of the fields it easily follows that and χ (2) For a lossless and weakly dispersive medium, when no external static magnetic fields are present, it can be shown that χ (2) i jk (ω 1 + ω 2 ; ω 1 , ω 2 ) is real and it is invariant under any permutation of the indices, provided the arguments are similarly permuted (see e.g. [25] for an explicit computation in perturbation theory or [26] for an indirect argument). In particular, in the limit of zero frequencies it is symmetric under a generic permutation of the indices.
The Pockels effect consists in the variation of the index of refraction at frequency ω when a static electric field is applied. To investigate the Pockels effect we thus have to study P P P (2) (ω) when the electric field is E E E(t) = E E E dc + Re[E E E opt (ω)e −iωt ], where we denoted by E E E dc the static (real) electric field and by E E E opt (ω) the component of frequency ω of the optical field. The relevant frequency components of the quadratic susceptibility are thus χ (2) i jk (ω; ω, 0) and the general symmetries previously discussed become now thus χ (2) i jk (ω; ω, 0) is invariant under permutation of the first two indices. Using these properties it is simple to obtain the following form for the polarization in the case of the Pockels effect For the sake of simplicity, we will suppress in the following the superscript of the optical field E E E opt that will be simply denoted by E E E. The phenomenon of second harmonic generation can also be described using similar relations, the main difference being that the relevant frequency components of the second order susceptibility are in that case χ (2) i jk (2ω; ω, ω), wich are symmetric under the exchange of the last two indices.
Note that, in addition to the previously discussed symmetries of the quadratic susceptibility, the lattice symmetry must also be added in the case of crystals: the tensor χ χ χ (2) is invariant under the symmetry group of the lattice. In particular, for centro-symmetric lattices (that are invariant under the inversion symmetry x x x → −x x x) the tensor χ χ χ (2) has to vanish, like all the invariant tensors with an odd number of indices.
In the following, when discussing lattice symmetries, it will obviously be convenient to work in the crystallographic frame, however attention has to be paid to the fact that this frame typically does not coincide with the device coordinate one, so that a change of reference frame is needed to obtain expressions of direct physical application [27].

The strain-induced Pockels effect
In the linear theory of elasticity, a small deformation x x x → x x x+u u u(x x x) is described by the symmetric strain tensor ε ε ε, defined by where u u u(x x x) represents the displacement of a material point [28]. In order to determine the relation between χ χ χ (2) and ε ε ε we will follow the same philosophy of effective field theories in theoretical physics: the most general expression compatible with the symmetries of the problem is considered, then the various possible terms are classified according to their "strength" retaining only the most relevant ones. The final relation will depend on a number of unknown constants, which are usually called "low energy constants" and have to be fixed by comparing with experimental data. As a starting point we assume that χ χ χ (2) is a local functional of ε ε ε, i.e. that χ χ χ (2) at point x x x depends only on the values of ε ε ε and its derivatives at x x x (the possible dependence of observable quantities on the strain gradient tensor is well known in the litature, see e.g. [29]). We then assume that this dependence is analytic and we develop everything in Taylor series, thus arriving to the following expression where we introduced the shorthand and where dots stand for terms involving higher derivatives of the strain tensor and higher orders of the Taylor expansion. The subscript | ε=0 means that the derivatives have to be computed at vanishing deformation.
If we now specialize to the case of centro-symmetric crystals, all the coefficients with an odd number of indices identically vanish, thus the previous expression becomes In particular all the terms depending on ε ε ε but not on its derivatives disappear. This is a direct consequence of the fact that a uniform strain does not break the inversion symmetry and thus cannot induce a non-vanishing quadratic susceptibility. The various term in the right hand side of Eq. (10) can be classified according to their power of the strain and their number of derivatives. In the limit of small deformation and (by writing explicitly the dependence on the frequencies and on the position) the leading contribution is and we thus expect a linear relation between the tensors χ χ χ (2) and ζ ζ ζ . The tensor T T T inherits some symmetries from χ χ χ (2) and ε ε ε: it is symmetric for α ↔ β and, as far as the Pockels effect is concerned, for i ↔ j (see Eqs. (5)). It is thus useful to adopt the contracted index notationT {i j}k{αβ }γ (see App. A), that points out the fact that only 324 of the 3 6 = 729 components of T T T are linearly independent. Since T T T is an invariant tensor for the lattice symmetry, this number can be further largely reduced. In the case of the lattice octahedral symmetry typical of the silicon crystal [8], the number of independent components is in fact only 15. The general procedure to identify the independent elements is reported in App. B, while the explicit form of Eq. (11) in the crystal reference frame is given in App. C.

The effective susceptibility
We have shown in the previous section that a linear relation between the χ χ χ (2) (x x x) tensor and the strain gradient ζ ζ ζ (x x x) has to be expected, however this local relation is not easily accessible by experiments. What is typically measured (see [12][13][14]) is the variation of the effective refractive index n eff for a propagating waveguide mode induced by the switching-on of the static electric field. The change of the effective refractive index due to the Pockels effect can easily be obtained (see App. D) and it is given by where the spatial dependence of the electric fields and of χ χ χ (2) is implied, A is the cross-section of the silicon waveguide and the normalization factor N is the active power of the optical mode that propagates along the waveguide, given by where i i i z denotes the unit vector parallel to the direction of propagation and A ∞ is the whole plane orthogonal to the waveguide. In a typical experimental setup, the metal contacts are much larger that the optical waveguide and the electrostatic field can be assumed constant in Eq. (12) [12][13][14].
In order to compare different experimental results, obtained using different devices, it is convenient to introduce an effective susceptibility that relates ∆n eff to E E E dc . It is worth noting that, since ∆n eff is a scalar quantity, an effective susceptibility defined in such a way will not be a 3−index tensor like χ χ χ (2) , but simply a vector. Another aspect that is usually not fully appreciated is the fact that the effective susceptibility not only provides a quantitative estimate of the non-vanishing properties of χ χ χ (2) , but it is also related to the effectiveness of the device in maximizing ∆n eff .
There are several definition of the effective susceptibility that are used in the literature and indeed, for this reason, some care is required in comparing results from different experimental groups (like, e.g., [12,13] and [14]). The simplest definition is just given by where the index k is the direction of the static field, like in Eq. (12). When only one component of the static field is different from zero, it is straightforward to derive the corresponding entry of the effective susceptibility vector. The definition in Eq. (14) appears to be the natural generalization to anisotropic optical waveguides of the simple result n∆n = χ (2) E dc , which is valid for the Pockels effect in homogeneous and isotropic materials, and it follows from n 2 = ε/ε 0 and ε = ε 0 (1 + χ (1) + 2χ (2) E dc ).
Using in this definition of effective susceptibility the relation of Eq. (11) leads to where we defined the weighted strain gradient as It can be noted that the frequency dependence of χ eff k (ω; ω, 0) originates from several sources: the possible explicit dependence of T T T (ω; ω, 0) on ω, the presence of n eff in the definition of the weighted strain gradients and, implicitly, through the form of the optical waveguide mode E E E(x x x) in Eq. (16).
Since the experimentally observed dependence of χ eff on the frequency is very mild (see [14]), it is reasonable to assume that the tensor T T T is frequency independent, i.e. that all the dependence of χ eff k on ω can be explained by the dependence of the weighted strain gradient on ω. It should be clear that this is not expected to be always true and this assumption can only be justified/disproved by its ability/inability to reproduce experimental data. When new, more precise, data will appear, it is conceivable that this hypothesis will have to be relaxed.

Investigation of strain-induced susceptibility in fabricated devices
In this section we investigate two waveguide fabricated and characterized in [12,13] and [14], respectively. The 15 independent entries of tensor T T T can be potentially computed by fitting the experimental results. However, because the results in the literature might be contaminated by the free carrier effect [15,17], we limit our analysis to some general properties of the straininduced χ χ χ (2) , showing the role played by the geometry to provide different effective second order susceptibility.  Figure 1. Cross-section of the strained silicon based MZI studied in [12,13] and [14], respectively. In (a) the slab waveguide cross-section described in [12,13]. The cases of waveguide width (w Si ) equal to 300 nm, 350 nm, 400 nm, 450 nm, and 500 nm have been investigated. In (b) the channel waveguide cross-section described in [14]. The cases of waveguide width (w Si ) equal to 385 nm, 435 nm, and 468 nm have been investigated. Pictures are not to scale.

Devices under investigation
The waveguide cross-section used in [12,13] is schematically shown in Fig. 1(a). A silicon rib waveguide was manufactured on a silicon-on-insulator (SOI) substrate with a 220 nm thick (100)-oriented top silicon layer over a 3 µm thick buried oxide. Silicon waveguides are fabricated by an etching process that leaves a slab thickness of 45 nm. A 350 nm thick Si 3 N 4 layer is deposited using remote plasma enhanced chemical vapor deposition and after the annealing process, a protective SiO 2 cladding layer is deposited on the top (850 nm thick).
The geometry adopted in [14] is slightly different and it is schematically shown in Fig. 1(b). A fully etched channel waveguide is fabricated on a 260 nm thick (100)-oriented top silicon layer with 2 µm buried oxide. The silicon waveguide is eventually covered by a single 700 nm layer of Si 3 N 4 . This choice was motivated by the fact that the protective SiO 2 cladding layer on Si 3 N 4 was observed to reduce the overall stress in the silicon wafer and the induced nonlinearity.

Strain simulation details
In order to investigate the consequences of Eqs. (15) and (16), we need to evaluate the strain gradient for the structures described in the previous section. Since silicon elastic properties significantly depend on the orientation of the crystalline structure, they have been taken into account in the mechanical simulation. To describe the deformation of silicon, one possibility could be to make use of the Hooke's law, i.e. the linear relation between strain and stress, which for materials with cubic symmetry involves only three independent components [30]. A more convenient description, that avoids tensorial transformation, is however the one that makes use of the orthotropic model [30]. A material is said to be orthotropic when it has at least two orthogonal planes of symmetry. Its elasticity can be described by a matrix that takes into account the fundamental elasticity quantities in the axes of interest: the Young's module (Y ), the Poisson's ratio (ν) and the shear modulus (G). In this work we use the letter Y for the Young's module instead of E, traditionally used, to avoid confusion with the electric field.
The most common use of orthotropic expressions for silicon is to provide the elasticity values in the frame of a standard (100) Similar relation can be derived for the device in [12,13], where x and z are switched in Eq. (17). The deformation of the silicon waveguides has been computed using COMSOL multiphysics tools [31]. Assuming 1 GPa compressive stress as the initial condition for the silicon nitride layer, the elastic strain has been computed for the structures shown in Fig.1(a) and Fig.1(b) as a function of the waveguide width and the wavelength. For the silicon waveguide and the buried silicon we used the values of the elastic modulus, the shear modulus and the Young's modulus in Eqs. (17). Since for the solid analysis a 2D model has been considered, the components ε xz and ε yz of the strain identically vanish; for the same reason, the derivative of the strain coefficients with respect to z are assumed equal to zero.

Results
In Fig. 2 we show the behavior of the strain components from which we can evaluate the strain derivatives, in Fig. 3 the results of the electromagnetic mode analysis, performed using COMSOL multi-physics tools [31]. The geometry taken into account is the waveguide cross section for the device described in [13]. The same results for the device described in [14] are respectively shown in Figs. 4 and 5. Focusing on the mode analysis, it is worth noting that in both cases the waveguides show a single mode behavior, however the E y and E z components are not negligible compared to the E x component and thus the mode is not purely transverse electric. The high value of E z is due to the high index step of the waveguides which requires in fact a full vectorial mode solver to be accurately computed [32]. Working in the device coordinate frame reported in Figs. 1(a) and 1(b) and with E E E dc directed along the y-axis, the relation between the effective susceptibility and the strain gradient can be written in the form where the coefficients c i (i = 1, . . . , 15) are the independent entries of the tensor T T T and the terms o i are linear combinations of the weighted strain gradients that will be called overlap functions in the following. Their explicit form is reported in App. C, Eq. (31) and it is the same for both the crystallographic axis orientation considered in [12,13] and [14]. For these orientations some of the overlap factors turn out to coincide: so that only 12 independent coefficients are needed in this case, instead of the 15 ones required for a generic device orientation. It is reasonable to expect that, for practical purposes, the number of numerical constants to be fixed can be further reduced. Indeed we numerically estimated the overlap functions corresponding to the experimental setups used in [13] and [14] and  Fig. 7 these most significants overlap are displayed for various waveguide width and wavelength. Since all the symmetries of the problem have been taken into account when determining the independent components of the tensor T T T , it is natural to expect all the c i to be about the same order of magnitude. As a consequence one expects that the coefficients that will be more important to reliably describe experimental results are the ones multiplying the dominant overlaps. This reduces the number of independent constants to be determined to 7 in our more conservative estimate.
From Figs. 6 and 7 we see that the dominant contributions are the ones denoted by o 9 , o 5p6 and o 2 . It is thus interesting to look at the specific form of these overlaps in Eq. (31). The largest field component is E x and it is natural to expect that these large overlaps are weighted by E x , which is indeed the case. Regarding the strain-gradient components, o 5p6 and o 9 involve ∂ ε xx /∂ y, which has been often assumed in the literature as the only relevant component of the strain-gradient; however o 2 is independent of this component and is related to ∂ ε yy /∂ y, whose importance has been overlooked so far.
As explained in the introduction, we cannot at present estimate the numerical values of the constants c i in Eq. (18) for lack of reliable (i.e. for which we can safely exclude a contamination from free carriers effects) experimental data to compare with. It is however interesting to note (a) (b) Figure 6. Behavior of the overlap functions for the waveguides under investigation with respect to waveguide width for (a) the device used by Chmielak et al. [13] and (b) the device used by Damas et al. [14]. In both cases only the most significant overlaps are plotted at λ = 1550nm. that, also without any knowledge of these constants, we can reach an important conclusion on a point that attracted some attention in the literature, i.e. the relevance of fabrication defects on the strain-induced Pockels effect. In particular in [13] was noted that the strain profile of the studied structure was quite sensitive to the presence of lateral cracks in the Si 3 N 4 overlayer. When such defects are present, we observe a strong increase of both strain components and strain gradients in proximity of the lateral cracks, as shown in Fig. 8. However, as can be seen from Fig. 9, the defects do not induce any sizable modification of the overlap functions. The physical explanation for this result is that the device regions in which the strain is significantly affected by the defect are also the regions where the electromagnetic field is weak, so that the net effect of the crack on the quantities in Eq. (16) is small. This simple fact was previously overlooked, probably due a confusion that is sometimes present in the literature between the strain-gradient dependence of the local non-linear susceptibility in Eq. (11) and that of the effective susceptibility, in which the strain-gradient tensor enters weighted by the electromagnetic field. Having explicitly checked that the overlaps does not appreciably changes when the Figure 8. Strain components profile ε xx (a), ε yy (b), ε zz (c) and ε xy (d) in the strained silicon waveguide [13] in presence of defect fabrication in the Si 3 N 4 slab. lateral crack is present, we can thus safely conclude (also without any knowledge of the numerical coefficient c i ) that the dependence of the effective susceptibility on this type of fabrication defects is very weak, as far as the strain-induced Pockels effect is concerned.

Conclusions
We have proposed an effective model that describes the strain-induced dielectric susceptibility in centro-symmetric crystals starting from symmetry arguments. The specific case of the Pockels effect in strained silicon has been investigated as a particularly interesting example, with potential applications in silicon photonics for implementing modulation and switching functions, however the idea of the approach is extemely general and can be applied to any centro-symmetric crystals (in fact also to non centro-symmetric ones) and to other physical phenomena, like e.g. second harmonic generation. The main result of our analysis consists in a simple relation between the second order dielectric susceptibility and the strain gradient tensor. We have shown that the effective second order susceptibility of strained silicon cannot be accurately described by considering only mechanical deformations but the combination of optical modes and mechanical stress analysis is required. As a result the effective susceptibility can be written as a linear combination of the weighted strain gradient components (defined in Eq. (16)) with fifteen independent coefficients in the general case. This approach put on firm ground some intuitive ideas present in the recent literature, that suggested a relation between the strain-induced effective susceptibility and the strain gradient.
In the experimental settings studied in [12,13] and [14] the number of independent coefficients to be determined to completely parametrize the strain-induced Pockels effect reduces to 12, moreover we gave indications that, for all practical purposes, this number can be effectively reduced to 7.
When new experimental measures of the effective strain-induced susceptibility (free from spurious effects like the ones related to free carriers) will be available, it will be possible to estimate these coefficients.This will have a considerable impact on the design and optimization of electro-optic modulators based on strained silicon, since it reduces the computation of electro-optic effect to standard mechanical and optical computations, thus allowing device optimization in terms of silicon waveguide geometries, crystallographic axes orientation and electrode position.

A. Contracted index notation
When dealing with tensors which are invariant under permutation of some indices it is convenient to introduce a compact notation, in order to simplify expressions. It is customary to define the contracted index notation (denoted by curly brackets) as follow Using this contracted notation a symmetric 2−index tensor V V V (i.e., V i j = V ji ) can be written as a vectorV V V with 6 entries, whose components are V {i j} with i ≤ j, and analogous simplifications occurs also for higher order tensors. It is important to note that, while the index contraction make some manipulations easier, the tensor transformation and multiplication properties are changed by this index replacement and some care is required. In contracted notation, the 6-index tensor T T T tensor introduced in Eq. (11) can be written (for the case of the Pockels effect, see the end of Sec. 3) as the the four index tensor where j-indices take value in 1, 2, 3, while i 2 , i 4 = 1, 2, 3, and i 1 , i 3 = 1, . . . , 6.

B. Symmetry analysis of the tensor T
In this appendix we report the algorithm used to write down the independent components of the 6−index tensor T T T introduced in Eq. (11). It is first of all convenient to rewrite the tensor T T T as a column vector T T T c by introducing a lexicographic labelling of the indices. The component ..i 6 (i 1 , . . . , i 6 = 1, 2, 3) corresponds to the component T c i where In order to impose the invariance properties of T T T , the general transformation law for tensors have to be translated in a form suitable to be applied to the vector T T T c . If R R R is the matrix representing the coordinate transformation, let us define the square matrix A A A as where the relation between i and i 1 , . . . , i 6 is the same as Eq. (22) (1) . . .
where I I I is the identity matrix.
To identify the independent elements of the tensor T T T , symmetries under index permutation must also be taken into account. For this purpose, let us introduce the column vectorT T T c asso- The relation betweenT T T andT T T c can be written in matrix form as where C C C is a 729 × 324 matrix, whose rows have only a single element different from zero and equal to 1. Combining Eq. (24) and Eq. (26) we arrive to where N N N is a rectangular matrix. Because N N N is real, rank(N N N) = rank (N N N t N N N (31)

D. The variation of the refraction index
In this appendix we report some details on the expression for the variation of the effective refraction index induced by a non-vanishing χ χ χ (2) , following an argument analogous to the one used in Refs. [35]. Let us consider the following form of the reciprocity theorem for a waveguide (see e.g. [36] §31-1): where the waveguide is directed along the z direction, A ∞ is the plane orthogonal to it and the vector F F F c is defined by In the following E E E 0 and H H H 0 will be the optical fields propagating in the guide when χ χ χ (2) ≡ 0, while E E E and H H H will be the corresponding fields when the nonlinear susceptibility is nonvanishing. It is simple to show that, when all the field have angular frequency ω, Eq. (32) can be rewritten in the form where P P P (2) is the polarization induced by the nonlinear susceptibility (P P P = P P P 0 + P P P (2) ) and A is the section of the wave guide in which χ χ χ (2) is non-vanishing. We will use for the unperturbed fields the form E E E 0 (r r r,t) = e e e 0 (x, y; ω 0 )e i(k 0 z−ωt) , H H H 0 (r r r,t) = h h h 0 (x, y; ω 0 )e i(k 0 z−ωt) , while we will assume for the perturbed fields the expressions E E E(r r r,t) = u(z)e e e 0 (x, y; ω 0 )e i(kz−ωt) , H H H(r r r,t) = u(z)h h h 0 (x, y; ω 0 )e i(kz−ωt) , thus assuming that the nonlinearity does not significantly affect the transverse modes. Using these forms of the fields in Eq. (34), together with Eq. (6), which can be written as P P P (2) = 2ε 0 χ χ χ (2) (ω; ω, 0) : E E EE E E dc , we obtain a differential equation for the envelope function u(z), namely: where we introduced the notation X = 2ε 0 A e e e * 0 · χ χ χ (2) : e e e 0 E E E dc dA A ∞ (e e e 0 × h h h * 0 + e e e * 0 × h h h 0 ) · i i i z dA .
By solving this equation and inserting the solution in Eq. (36) we finally obtain k = k 0 + ωX, thus corresponding to the variation ∆n eff = cX of the effective refraction index.