FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides

A deep insight into the inherent anisotropic optical properties of silicon is required to improve the performance of silicon-waveguide-based photonic devices. It may also lead to novel device concepts and substantially extend the capabilities of silicon photonics in the future. In this paper, for the first time to the best of our knowledge, we present a three-dimensional finite-difference time-domain (FDTD) method for modeling optical phenomena in silicon waveguides, which takes into account fully the anisotropy of the third-order electronic and Raman susceptibilities. We show that, under certain realistic conditions that prevent generation of the longitudinal optical field inside the waveguide, this model is considerably simplified and can be represented by a computationally efficient algorithm, suitable for numerical analysis of complex polarization effects. To demonstrate the versatility of our model, we study polarization dependence for several nonlinear effects, including self-phase modulation, cross-phase modulation, and stimulated Raman scattering. Our FDTD model provides a basis for a full-blown numerical simulator that is restricted neither by the single-mode assumption nor by the slowly varying envelope approximation. © 2010 Optical Society of America OCIS codes:(040.6040) Silicon; (250.4390) Nonlinear optics, integrated optics; (160.1190) Anisotropic optical materials; (050.1755) Computational electromagnetic methods; (230.4320) Nonlinear optical devices. References and links 1. H. Iwai and S. Ohmi, “Silicon integrated circuit technology from past to future,” Microelectron. Reliab. 42, 465–491 (2002). 2. S. E. Thompson and S. Parthasarathy, “Moore’s law: The future of Si microelectronics,” Mater. Today 9, 20–25 (2006). 3. M. Schulz, “The end of the road for silicon?” Nature 399, 729–730 (1999). 4. B. Jalali, “Can silicon change photonics?” Phys. Status Solidi (a) 2, 213–224 (2008). 5. R. A. Soref, “The past, present, and future of silicon photonics,” IEEE J. Sel. Top. Quantum Electron. 12, 1678– 1687 (2006). 6. G. T. Reed and A. P. Knight, Silicon Photonics: An Introduction (Wiley, Hoboken, NJ, 2004). 7. L. Pavesi and D. J. Lockwood, Eds., Silicon Photonics(Springer, New York, 2004). 8. B. Jalali, S. Yegnanarayanan, T. Yoon, T. Yoshimoto, I. Rendina, and F. Coppinger, “Advances in silicon-oninsulator optoelectronics,” IEEE J. Sel. Top. Quantum Electron. 4, 938–947 (1998). 9. B. Jalali and S. Fathpour, “Silicon Photonics,” J. Lightwave Technol. 24, 1460–4615 (2006). 10. R. Claps, D. Dimitropoulos, and B. Jalali, “Stimulated Raman scattering in silicon waveguides,” Electron. Lett. 38, 1352–1354 (2002). #133336 $15.00 USD Received 12 Aug 2010; revised 22 Sep 2010; accepted 22 Sep 2010; published 23 Sep 2010 (C) 2010 OSA 27 September 2010 / Vol. 18, No. 20 / OPTICS EXPRESS 21427 11. M. Dinu, F. Quochi, and H. Garcia, “Third-order nonlinearities in silicon at telecom wavelengths,” Appl. Phys. Lett. 82, 2954–2956 (2003). 12. E. Dulkeith, Y. A. Vlasov, X. Chen, N. C. Panoiu, and R. M. Osgood, Jr., “Self-phase-modulation in submicron silicon-on-insulator photonic wires,” Opt. Express 14, 5524–5534 (2006). 13. I.-W. Hsieh, X. Cheng, J. I. Dadap, N. C. Panoiu, R. M. Osgood, S. J. McNab, and Y. A. Vlasov, “Ultrafast-pulse self-phase modulation and third-order dispersion in Si photonic wire-waveguides,” Opt. Express 14, 12380– 12387 (2006). 14. I.-W. Hsieh, X. Chen, J. I. Dadap, N. C. Panoiu, R. M. Osgood, Jr., S. J. McNab, and Y. A. Vlasov, “Crossphase modulation-induced spectral and temporal effects on co-propagating femtosecond pulses in silicon photonic wires,” Opt. Express 15, 1135–1146 (2007). 15. T. Liang, L. Nunes, T. Sakamoto, K. Sasagawa, T. Kawanishi, M. Tsuchiya, G. Priem, D. V. Thourhout, P. Dumon, R. Baets, and H. Tsang, “Ultrafast all-optical switching by cross-absorption modulation in silicon wire waveguides,” Opt. Express 13, 7298–7303 (2005). 16. Q. Lin, J. Zhang, P. M. Fauchet, and G. P. Agrawal, “Ultrabroadband parametric generation and wavelength conversion in silicon waveguides,” Opt. Express 14, 4786–4799 (2006). 17. Y.-H Kuo, H. Rong, V. Sih, S. Xu, M. Paniccia, and O. Cohen, “Demonstration of wavelength conversion at 40 Gb/s data rate in silicon waveguides,” Opt. Express 14, 11721–11726 (2006). 18. H. Fukuda, K. Yamada, T. Shoji, M. Takahashi, T. Tsuchizawa, T. Watanabe, J.-I. Takahashi, and S.-I. Itabashi, “Four-wave mixing in silicon wire waveguides,” Opt. Express 13, 4629–4637 (2005). 19. R. Claps, V. Raghunathan, D. Dimitropoulos, and B. Jalali, “Anti-Stokes Raman conversion in silicon waveguides,” Opt. Express 11, 2862–2872 (2003). 20. V. Raghunathan, R. Claps, D. Dimitropoulos, and B. Jalali, “Parametric Raman wavelength conversion in scaled silicon waveguides,” J. Lightwave Technol. 23, 2094–2102 (2005). 21. Ö. Boyraz and B. Jalali, “Demonstration of a silicon Raman laser,” Opt. Express 12, 5269–5273 (2004). 22. H. Rong, A. Liu, R. Jones, O. Cohen, D. Hak, R. Nicolaescu, A. Fang, and M. Paniccia, “An all-silicon Raman laser,” Nature433, 292–294 (2005). 23. R. Claps, D. Dimitropoulos, V. Raghunathan, Y. Han, and B. Jalali, “Observation of stimulated Raman amplification in silicon waveguides,” Opt. Express 11, 1731–1739 (2003). 24. Q. Xu, V. R. Almeida, and M. Lipson, “Demonstration of high Raman gain in a submicrometer-size silicon-oninsulator waveguide,” Opt. Lett. 30, 35–37 (2005). 25. T. K. Liang, L. R. Nunes, M. Tsuchiya, K. S. Abedin, T. Miyazaki, D. Van Thourhout, W. Bogaerts, P. Dumon, R. Baets, and H. K. Tsang, “High speed logic gate using two-photon absorption in silicon waveguides,” Opt. Commun.265, 171–174 (2006). 26. C. Manolatou and M. Lipson, “All-optical silicon modulators based on carrier injection by two-photon absorption,” J. Lightwave Technol. 24, 1433–1439 (2006). 27. D. J. Moss, L. Fu, I. Littler, and B. J. Eggleton, “Ultrafast all-optical modulation via two-photon absorption in silicon-on-insulator waveguides,” Electron. Lett. 41, 320–321 (2005). 28. R. Dekker, A. Driessen, T. Wahlbrink, C. Moormann, J. Niehusmann, and M. F örst, “Ultrafast Kerr-induced all-optical wavelength conversion in silicon waveguides using 1.55 μm femtosecond pulses,” Opt. Express 14, 8336–8346 (2006). 29. R. Espinola, J. Dadap, R. Osgood, Jr., S. McNab, and Y. Vlasov, “C-band wavelength conversion in silicon photonic wire waveguides,” Opt. Express 13, 8336–8346 (2005). 30. R. Salem, G. E. Tudury, T. U. Horton, G. M. Carter, and T. E. Murphy, “Polarization-insensitive optical clock recovery at 80 Gb/s using a silicon photodiode,” IEEE Photon. Technol. Lett. 17, 1968–1970 (2005). 31. Ö. Boyraz, P. Koonath, V. Raghunathan, and B. Jalali, “All optical switching and continuum generation in silicon waveguides,” Opt. Express 12, 4094–4102 (2004). 32. I.-W. Hsieh, X. Chen, X. Liu, J. I. Dadap, N. C. Panoiu, C.-Y. Chou, F. Xia, W. M. Green, Y. A. Vlasov, and R. M. Osgood, “Supercontinuum generation in silicon photonic wires,” Opt. Express 15, 15242–15249 (2007). 33. E. K. Tien, N. S. Yuksek, F. Qian, and O. Boyraz, “Pulse compression and modelocking by using TPA in silicon waveguides,” Opt. Express 15, 6500–6506 (2007). 34. L. Yin, J. Zhang, P. M. Fauchet, and G. P. Agrawal, “Optical switching using nonlinear polarization rotation inside silicon waveguides,” Opt. Lett. 34, 476–478 (2009). 35. B. A. Daniel and G. P. Agrawal, “Vectorial nonlinear propagation in silicon nanowire waveguides: Polarization effects,” J. Opt. Soc. Am. B27, 956–965 (2010). 36. I. D. Rukhlenko, I. L. Garanovich, M. Premaratne, A. A. Sukhorukov, G. P. Agrawal, and Yu. S. Kivshar, “Polarization rotation in silicon waveguides: Analytical modeling and applications,” IEEE Photonics Journal 2, 423– 435 (2010). 37. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, Boston, 2005). 38. M. Born and E. Wolf,Principles of Optics , 7th ed. (Cambridge U. Press, Cambridge, 1999). 39. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). #133336 $15.00 USD Received 12 Aug 2010; revised 22 Sep 2010; accepted 22 Sep 2010; published 23 Sep 2010 (C) 2010 OSA 27 September 2010 / Vol. 18, No. 20 / OPTICS EXPRESS 21428 40. R. M. Joseph and A. Taflove, “FDTD Maxwell’s equations models for nonlinear electrodynamics and optics,” IEEETrans. Antennas Propag. 45, 364–374 (1997). 41. P. M. Goorjian, A. Taflove, R. M. Joseph, and S. C. Hagness, “Computational modeling of femtosecond optical solitons from Maxwell’s equations,” IEEE J. Quantum Electron. 28, 2416–2422 (1992). 42. F. L. Teixeira, “Time-domain finite-difference and finite-element methods for Maxwell equations in complex media,” IEEE Trans. Antennas Propag. 56, 2150–2166 (2008). 43. J. Schneider and S. Hudson, “A finite-difference time-domain method applied to anisotropic material,” IEEE Trans. Antennas Propag. 41, 994–999 (1993). 44. R. Holland, “Finite-difference solution of Maxwell’s equations in generalized nonorthogonal coordinated,” IEEE Trans. Nucl. Sci. NS-30, 4589–4591 (1983). 45. L. X. Dou and A. R. Sebak, “3D FDTD method for arbitrary anisotropic materials,” Microwave Opt. Technol. Lett. 48, 2083–2090 (2006). 46. D.-X. Xu, P. Cheben, D. Dalacu, A. Del âge, S. Janz, B. Lamontagne, M.-J. Picard, and W. N. Ye, “Eliminating the birefringence in silicon-on-insulator ridge waveguides by use of cladding stress,” Opt. Lett. 29, 2384–2386 (2004). 47. J. Pastrnak and K. Vedam, “Optical anisotropy of silicon single crystals,” Phys. Rev. B 3, 2567–2571 (1971). 48. T. Chu, M. Yamada, J. Donecker, M. Rossberg, V. Alex, and H. Riemann, “Optical anisotropy in dislocation-free silicon single crystals,” Microelectron. Eng. 66, 327–332 (2003). 49. R. W. Boyd,Nonlinear Optics, 2nd ed. (Academic Press, Boston, 2003). 50. E. B. Graham and R. E. Raab, “Light propagation in cubic and other anisotropic crystals,” Proc


Introduction
The microelectronics industry has thrived on the unique electrical and mechanical properties of silicon, which have enabled significant miniaturization and mass-scale integration of electronic components on a single silicon chip [1,2].Unfortunately, this miniaturization process cannot be continued much beyond the current level because of the inherent heat generation and the metal-interconnects bottleneck effect [3].A promising way to overcome these fundamental limitations is to use photons, instead of electrons, as information carriers.This would lead to a gradual replacement of electronic integrated circuits (ICs) by photonic ICs.Silicon is well suited for realizing photonic ICs, owing to the recent advances in the fabrication technologies and the ease of integrating electronics with silicon photonics [4][5][6][7][8][9].A relatively high refractiveindex contrast provided by the silicon-on-insulator (SOI) technology allows one to confine the optical mode to subwavelength dimensions and to realize nanoscale optical waveguides (often called photonic nanowires), for on-chip optoelectronic applications [8].
To study the impact of anisotropic nonlinear effects on light propagation in silicon, a complex numerical simulator that takes into account the tensorial nature of third-order susceptibility is required.The finite-difference time-domain (FDTD) scheme is the best choice for sophisticated modeling, since it allows Maxwell's equations to be exactly solved in three dimensions with few assumptions [37].This scheme is suitable for analyzing complex waveguide geometries with few-cycle pulses, which is beyond the capability of other popular modeling techniques.In this paper, we present an FDTD model of optical phenomena in silicon waveguides, which incorporates the anisotropy of silicon nonlinearities and the vectorial nature of the electromag- (C) 2010 OSA netic field.In Section 2, we describe the structure of the general FDTD algorithm for nonlinear anisotropic media.Starting with Maxwell's equations, we outline the steps required for obtaining the update equations for the electromagnetic field and set up notations employed in the rest of the paper.In Section 3, we consider the relation between the material polarization and the electric field inside a silicon waveguide and derive the contributions of different optical effects with emphasis on the underlying physics.The full material polarization is used in Section 4 to derive the update equations for the electric field in the three-dimensional case.We then simplify these equations to obtain a computationally efficient algorithm that can be used to simulate SOI-based photonic devices.In Section 5, we apply the simplified algorithm to investigate the impact of optical anisotropy on pulse propagation in silicon waveguides.We summarize our study and conclude the paper in Section 6.

Three-dimensional FDTD scheme for nonlinear anisotropic media
Classical electromagnetic theory is based on Maxwell's equations, which admit exact analytical solutions only for a limited number of relatively simple problems [38].In modern nonlinear optics, these equations are usually solved numerically using different computational schemes.The FDTD method is one of the most widely used methods for this purpose, due to the simplicity of its implementation, its applicability to arbitrary scattering media (e.g., media exhibiting anisotropic and/or nonlinear responses), and the high accuracy of the generated results [37].
In order to establish the notations used in the following sections, we sketch the main steps involved in the FDTD algorithm that solves the source-free Maxwell's equations in their differential form [38] ∂ B(r,t) where B(r,t) and D(r,t) are the magnetic and electric inductions, and E(r,t) and H(r,t) are the electric and magnetic fields.For nonmagnetic media, the constitutive relations between the four vector fields in Eq. ( 1) are given by where µ 0 and ε 0 are, respectively, the permeability and permittivity of vacuum.The response of an inhomogeneous anisotropic medium to an external electric field is taken into account by the material polarization vector P(r,t).
An arbitrary vector field V = {B, D, H, E, P} can be split into its Cartesian components V ξ along the axis ξ = {α, β , γ} of the FDTD coordinate system as V(r,t) = αV α (r,t) + βV β (r,t) + γV γ (r,t), where α, β , and γ are the unit vectors.In the FDTD scheme, each of these components are discretized in time and space domains according to the recipe of Yee [37,39].
The discretization scheme proposed by Yee can be visualized using the cubic unit cell and the time line shown in Fig. 1.The electric and magnetic fields are defined on different points of the unit cell and at alternative time steps.Mathematically, each component V ξ (r,t) becomes a function of four discrete arguments, where {i, j, k, n} ∈ N, and ∆ξ and ∆t are the space and time increments.The values of d µ (µ = 1, 2, 3, 4), corresponding to the arrangement adopted in Fig. 1, are given in Table 1.It should be noted that not all components shown in Fig. 1 belong to the same unit cell.In order to obtain the update equations for the discrete electromagnetic field, one needs to replace the time and space derivatives in Eq. ( 1) by finite differences.With this approach, the update equations for the vectors B, H, and D can be readily derived from Eqs. ( 1) and (2).For example, for α components the update equations are: where S = c∆t/∆ξ is the Courant stability number and c is the speed of light in vacuum.For brevity, from here on we omit space indexes on the right side of update equations that coincide with those on the left side.With these and similar equations for other field components, we can calculate the values of B ξ , H ξ , and D ξ at any time step, if their values and the values of E ξ are known at the previous time step.
To update the electric field component E ξ (ξ = α, β , γ) at the time step n + 1, we need to Positions of electric field components in the FDTD grid.The values of components E β and E γ at the node (i, j, k) are calculated by averaging their four values over the nodes know the components D ξ and P ξ at the time steps n + 1 and n.According to Eq. ( 2), the update equation takes the form In turn, the material polarization can be calculated if the electric field inside the medium is known.This polarization includes contributions from different physical phenomena, which may be either isotropic (i) or anisotropic (a), i.e., P(r,t) = ∑ i P (i) + ∑ a P (a) .
The relations between the electric field and the induced polarizations P (i) and P (a) are available either in the time or the frequency domain, depending on the specific models of the dielectric response [40][41][42].The frequency-domain relations can be transformed into the time domain either by using the standard replacement iω → −∂ /∂t or by employing the higher-order Fourier transforms.The resulting time-domain relations can be converted into the update equations for material polarization using finite differences for the derivatives.
If the material response is isotropic, then the ξ component of the polarization vector depends only on the ξ component of the electric field, and we can write symbolically where f (i) ξ (. ..) denotes the functional dependence.In the case of the anisotropic material response, the value of P (a) ξ is determined by all three components of the electric field such that This equation assumes that the material response is local and that the values of three components of the electric field are to be taken at the same grid node as the polarization component.ξ [43][44][45].For example, the components E β and E γ are projected onto the position of the component P (a) α using the relations (see Fig. 2) Equations ( 3)-( 5) are solved iteratively for a fixed value of D until the desired accuracy is achieved.

Optical phenomena in silicon waveguides
As we have seen, both the isotropic and anisotropic contributions of the optical response of a dielectric medium can be easily handled with the FDTD scheme.However, its implementation requires establishing a functional relation between the material polarization induced in a silicon waveguide and the electric field of the propagating mode inside the waveguide.
A number of linear and nonlinear optical phenomena should be taken into account for accurate modeling of silicon-based photonic devices.The major linear effects are dispersion, absorption, and scattering losses.Linear absorption in silicon is isotropic and much stronger than that in fibers; however, its effect on optical propagation is usually weak as opposed to the impact from nonlinear losses.Similar to optical fibers, linear dispersion plays a significant role in silicon waveguides when the spectrum of a propagating pulse is sufficiently broad.Due to a high refractive-index contrast, the dispersive properties of a SOI waveguide depend strongly on its geometry; this dependence leads to a relatively large waveguide dispersion.In particular, large modal birefringence (∼ 10 −3 ) can arise in SOI waveguides with asymmetric cross sections [46], even though the intrinsic optical anisotropy of silicon is relatively low (∼ 10 −6 ) [47,48].Modal birefringence and waveguide dispersion are automatically included in the FDTD scheme when the linear optical response of silicon is characterized by a scalar susceptibility, χ(1) (ω).
If there is no inhomogeneous mechanical stress breaking the inversion symmetry of the crystal field, the second-order electro-optic effect is absent in silicon [49][50][51].In this paper, we only consider silicon waveguides that are free from mechanical stresses and ignore all second-order nonlinear effects.
The nonlinear interaction of an optical field with bound electrons of silicon atoms results in the Kerr effect, two-photon absorption (TPA), and Raman scattering [52][53][54].These effects are of the third order with respect to the applied field.They have anisotropic character and are described by two fourth-rank tensors, χe κλ µν and χR κλ µν [49,53,55].The process of TPA may generate a considerable number of free carriers, which absorb light and reduce the effective refractive index of a silicon waveguide.In addition, TPA generates heat and raises the waveguide temperature leading to an increase in the refractive index; this phenomenon is referred to as the thermo-optic effect (TOE) in the literature [56][57][58][59].The impact of TOE on pulse propagation is typically much stronger than that due to free-carrier dispersion (FCD) [56].Free-carrier effects and TOE are isotropic, but strongly nonlinear; their strength grows in proportion to the fifth power of the electric field and can be approximately described by an intensity-dependent susceptibility, χFC (|E| 4 , ω).With the preceding nonlinear processes in mind, the polarization Pκ (r, ω) induced inside a silicon waveguide by the electric field Ẽκ (r, ω) can be represented in the form [49,53] where χ(3) κλ µν = χe κλ µν + χR κλ µν , ω 3 = ω − ω 1 − ω 2 , and we used the following convention for the Fourier transform: Equation ( 6) is written in the FDTD coordinate system, whose axes do not generally coincide with the crystallographic axes of the waveguide (see Fig. 3).Therefore before using this equation, it is necessary to transform the third-order susceptibility tensor from the crystallographic coordinates (x, y, and z) into the FDTD coordinates (α, β , and γ).The general transformation for a fourth-rank tensor is given by [60]: where the Latin indices denote the crystallographic axes and a is the transformation matrix.
It is convenient to choose the FDTD axes along the edges of the silicon waveguide.Without loss of generality, we select the α axis along the waveguide length and take the other two axes to form a right-handed coordinate system.Since silicon waveguides are often fabricated on the (001) plane, we assume that the γ axis coincides with the [001] direction (see Fig. 3).In this case, the FDTD axes α and β are obtained by rotating the crystallographic axes x and y about the z axis.If the rotation angle is ϑ and the rotation is performed as shown in Fig. 3, then the transformation matrix is of the form As an example, ϑ = π/4 corresponds to the widely used waveguides fabricated along the [110] direction [53,61].

Linear absorption and linear dispersion
Linear optical properties of silicon waveguide are described by the complex-valued dielectric function ε(ω The real part of this function accounts for dispersion of the linear refractive index n(ω) ≈ ε′ (ω), while the imaginary part is responsible for linear losses, ε′′ (ω) ≈ cn 0 α L /(2ω), where α L is the linear-loss coefficient and n 0 ≡ n(ω 0 ) is the refractive index at the absolute maximum of the input field spectrum located at frequency ω 0 .In practice, α L includes scattering losses occurring at waveguide interfaces.
Dispersion of the refractive index has been measured for silicon crystals over a broad spectral range.The experimental data is often approximated with an appropriately chosen analytical formula.One of the simplest formulae for n 2 (ω) is given by the 2s-pole Sellmeier equation [62][63][64] where a j and ω j are the fitting parameters.In Section 5, we use this equation in the specific case s = 2 with a 1 = 9.733, a 2 = 0.936, ω 1 = 1032.49THz, and ω 2 = 817.28THz.These values provide an accurate fit of the refractive-index data for silicon in the wavelength range from 1.2 to 2 µm [65].

The Kerr effect and two-photon absorption
At high optical intensities, electrons residing at the outer shells of silicon atoms start moving anharmonically.This leads to an increase by ∆n NL of the refractive index in proportion to the local intensity I of light.This phenomena, referred to as the Kerr effect, is characterized by the nonlinear coefficient n 2 = ∆n NL /I ≈ 6 × 10 −5 cm 2 /GW.The Kerr effect is responsible for self-phase modulation (SPM), soliton formation, supercontinuum generation, and modulation instability, among other things [12,64,66,67].
In the telecommunication band around 1.55 µm, energy of infrared photons (about 0.8 eV) exceeds one half of the indirect band gap in silicon (0.56 eV).As a result, silicon exhibits twophoton absorption (TPA) accompanied with the emission of transverse optical (TO) phonons [68,69] and generation of free carriers.Despite the fact that TPA is a nonlinear-loss mechanism and is undesirable in majority of photonic devises, it has been successfully employed to create ultrafast optical modulators and switches [25][26][27]33].
In the crystallographic coordinate system, the anisotropy tensor E klmn has the following wellknown form [49,53]: where δ kl is the Kronecker delta function and the real number ρ characterizes anisotropy of the electronic susceptibility (ρ ≈ 1.27 near λ = 1.55 µm).Transforming E klmn according to Eqs. ( 7) and ( 8), we obtain the following nonzero elements of the anisotropy tensor in the FDTD coordinates: where A ϑ = 1 2 (ρ − 1) sin 2 (2ϑ ) and B ϑ = 1 4 (ρ − 1) sin(4ϑ ).Using the preceding results and Eq. ( 6), we obtain the Kerr-induced material polarization in the time domain in the form where the auxiliary functions S K κ are given by In the approximation of instantaneous electronic response, the anisotropy of TPA is identical to that of the Kerr effect.Owing to this fact, the TPA-induced polarization can be written similar to Eq. ( 9), In Section 4, we convert this equation into the time domain and use it together with Eqs. ( 9) and (10)

Raman scattering
If an intense optical field at frequency ω p propagates through a silicon waveguide, it generates another two fields at frequencies ω p ± Ω R through a phenomenon known as the spontaneous Raman scattering [55,70].These fields are induced by the spontaneous emission and absorption of optical phonons of frequency Ω R , the so-called Raman shift.The process of Raman scattering is substantially intensified (stimulated) if a second optical field of frequency ω s = ω p − Ω R is launched into the waveguide.In this case, one says that the energy is transferred between the fields via the process of stimulated Raman scattering (SRS) [53][54][55].
In silicon, Raman scattering is produced by optical phonons near the Brillouin zone center [71].The three phonon modes (almost dispersionless) have the same energy of about 65 meV.This energy determines the population of the phonon modes and relative intensities of red-shifted (Stokes) and blue-shifted (anti-Stokes) lines in the optical spectrum.At room temperature (T ≈ 25 meV), the Stokes line is nearly 13.5 times stronger than the anti-Stokes line.Even though the processes of Stokes and anti-Stokes scattering substantially differ in strength, it is important for the accuracy of numerical data that both are easily accounted for in the FDTD scheme.
The dispersion of Raman scattering is adequately described by treating silicon cores as classical oscillators [49,54].With this model, one can write the Raman susceptibility in the form [53] where the dimensionless tensor R κλ µν describes the anisotropy of the Raman scattering and the function H(ω) represents the Raman gain profile, Here, ξ R = 2ε 0 n 0 c 2 g R (ω)/ω, g R (ω) is the Raman gain coefficient, Ω R = 15.6 THz, and 2Γ R ≈ 100 GHz is the gain bandwidth [72].The ratio g R (ω)/ω is almost independent of ω since g R (ω) ∝ ω [53].In the 1.55-µm region, the values of g R reported by different experimental groups lie in the range from 4.3 to 76 cm/GW [73].The fact that these values are more than 450 and 7500 times larger than the Raman gain coefficient in optical fibers (∼ 0.01 cm/GW), enables net gain from SRS even in the presence of strong nonlinear absorption [21,74,75].
In the crystallographic coordinate system, the anisotropy tensor for silicon has the form This expression shows, for example, that the Raman scattering does not occur between two optical beams polarized along the same crystallographic axis (R k j j j = δ k j + δ k j − 2δ k j = 0).Therefore, if a single pulse polarized along any principal axis propagates through a silicon waveguide, the SRS-induced redistribution of energy within the pulse bandwidth does not occur.
As before, we transform R klmn into the FDTD coordinate system using Eqs.( 7) and (8).It is easy to see that the 24 nonzero elements of the anisotropy tensor are given by where C ϑ = sin 2 (2ϑ and D ϑ = 1  2 sin(4ϑ ).Using Eq. ( 12) and introducing a new function H(t) = F −1 [ H(ω)], we arrive at the following time-domain expression for the Raman susceptibility: With this result, the Raman-induced material polarization can be written in the following form suitable for the FDTD implementation: where

Free-carrier and thermo-optic effects
At high optical powers, TPA produces a large number of free electrons and holes, which give rise to another optical-loss mechanism known as free-carrier absorption (FCA).Besides, free carriers reduce the effective refractive index of silicon waveguides; this effect is referred to as free-carrier dispersion (FCD).The strengths of FCD and FCA depend on the free-carrier density N and are proportional to the real and imaginary parts of the susceptibility as where ∆n FC (N) and ∆α FC (N) stand, respectively, for the free-carrier-induced change in the refractive index and for the loss coefficient.We assume that carriers are neither injected into the waveguide nor moved out of it, and the densities of electrons and holes are equal.In this case, it is common to use the following empirical formulas [53]: The evolution of free-carrier density is governed by the interplay between two processes: carrier generation via TPA and carrier recombination through all possible nonradiative channels.These processes are described by the rate equation where τ c is the free-carrier lifetime and I(r,t) = 1 2 cε 0 0 |E(r,t)| 2 is the optical intensity.Since silicon is nondirect gap semiconductor, TPA is accompanied by the emission of phonons.In addition, phonons are emitted during the electron-hole recombination.As a result of these two processes, the temperature and refractive index of a silicon waveguide increase.In silicon, the effect of thermal changes in the refractive index-the thermo-optic effect (TOE)-is much stronger than FCD [56].This is a consequence of the large value of thermo-optic coefficient κ = 1.86 × 10 −4 K −1 , which relates the nonlinear change in the refractive index to device temperature.The fact that all energy dissipated through TPA is eventually converted to heat, allows one to account for TOE by reducing FCD coefficient ζ by the amount where θ ≈ 1 µs is the thermal dissipation time, C ≈ 0.7 J/(g×K) is the thermal capacity, and ρ = 2.3 g/cm 3 is the density of silicon.Generally, both FCD and TOE should be included in the FDTD simulator because both effects are crucial for phase-sensitive problems.

Update equations for the electric field in silicon waveguides
The general FDTD scheme discussed in Section 2 will now be applied to silicon by employing the various susceptibilities specified in the preceding section.This involves discretization of the second constitutive relation in Eq. ( 2) and derivation of algebraic equations for updating the electric field components.

Three-dimensional implementation
We start by considering the general situation in which dimensions of the silicon waveguide and the polarization state of the input field are arbitrary.In this case, all six components of the electromagnetic field should be accounted for in the FDTD scheme.In the frequency domain, the electric field is related to the material polarization of the silicon waveguide as where κ = α, β , and γ.Multiplying both sides of this equation by iω and using Eqs.( 9) and ( 14), the standard replacement of the factor −iω with the finite-difference operator in the time domain yields: where, for compactness, we have omitted the space subscripts .
Conversion into the time domain of the equation for the jth component, ω 2 j − ω 2 PLD κ j = ε 0 a j ω j Ẽκ , gives us the differential equation This equation is readily discretized with the central difference scheme to obtain In a similar fashion, we use Eq. ( 13) to find that the evolution of nine auxiliary variables Ξ R κλ in Eq. ( 14) is governed by the differential equation whose discretized version takes the form The relations between the functions S R κλ and the electric field components are given in Section 3.3.
Since Eq. ( 15) is the first-order differential equation, the functions ∆α FC (N) and ∆n FC (N) are defined at the non-integer time steps.The carrier density at the time step n + 1 2 can be found using its value at the time step n − 1 2 and the stored components E κ .The updated values of the functions ∆α FC and ∆n FC are thus calculated using the equation For stability and accuracy of the FDTD simulations, time step ∆t should be much less than the duration of a single optical cycle.Therefore, we can simplify the algorithm in Eq. ( 16) with the approximation without noticeably affecting the FDTD results.Rearranging the terms in Eq. ( 16), we finally arrive at the following coupled update equations for the electric field components: Notice that each of the auxiliary functions S K κ depends polynomially on all three components of the electric field through Eq. (10).For this reason, Eq. ( 20) cannot be solved analytically for the unknown component E κ | n+1 (κ = α, β , γ).Rather, values of the functions P LD κ j , Ξ R κλ , and N are updated using Eqs.( 17)-( 19) and the electric field components calculated at the previous time step.With these values, Eq. ( 20) is then solved iteratively to update the electric field.

A simplified FDTD model for silicon waveguides
Tight confinement of optical mode inside SOI waveguides results in strong variations of the magnetic field over the mode profile.If the incident beam is polarized perpendicular to the waveguide axis α, such variations generate a relatively large longitudinal component E α of the electric field as an optical beam propagates through the waveguide [35,76].Generally, E α cannot be ignored even in a single-mode waveguide for which spatial variations of the magnetic field are relatively weak.The reason for this is the SRS-induced coupling between the electric field components, which is expressed by the sum in Eq. (20).
However, the buildup of the longitudinal electric field in the vicinity of the mode center can be ignored when the spectrum of the incident pulse is much narrower than the Raman shift of 15.6 THz.This condition is satisfied for optical pulses of widths > 100 fs.For shorter pulses, the adopted models of the Kerr effect and TPA may need be modified to allow for a finite response time of the bound electrons.Indeed, as can be seen from Eq. ( 14), the Raman polarization P R α is not produced by the transverse components E β and E γ when the waveguide is fabricated along the crystallographic direction determined by the angle ϑ j = π 4 j, where j = 0, 1, 2, 3 (see Fig. 3), resulting in D ϑ j = 0.
In this section we consider a one-dimensional FDTD domain in which the transverse spatial distributions of the TE and TM modes are ignored, and Maxwell's equations are solved only along the waveguide length (the axis α).It amounts to assuming that the incident field is polarized in the β -γ plane and the silicon waveguide is fabricated along one of the directions ϑ j .Further, the waveguide dimensions are assumed to be large enough that the axial field components remain small all along the waveguide length, i.e., the optical field preserves its transverse nature along the mode center.The effects of waveguide dimensions are then included only through the effective mode index n 0 and the effective mode area.In practice, the results of such a simplified FDTD model should be applicable for square-shape waveguides with dimensions larger than the optical wavelength inside silicon (λ /n 0 ).
With these limitations in mind, the update equation ( 20) can be simplified considerably near the mode center and takes the form where κ = β or γ, ν = κ, and the functions B ± κ should be used with E α = 0.These equations show that the preceding two assumptions simplify the general three-dimensional FDTD model substantially and reduce it to a pair of one-dimensional FDTD models [37,43].These two one-dimensional models describe propagation of the β -polarized and γ-polarized transverse electromagnetic modes, which are coupled through the Kerr effect, TPA, and SRS.Free-carrier effects do not couple Eqs. ( 21), as they depend on the electric field at the previous time steps,  n and n − 1.The one-dimensional nature of the coupled FDTD models is also evident by their Yee cells shown in Fig. 4.
It is important to keep in mind that the reduced implementation of the FDTD scheme cannot be used to analyze anisotropy stemming from the asymmetry of the waveguide cross section (e.g., modal birefringence, FCA-induced polarization rotation, etc. [35]).However, it allows one to examine the effects of intrinsic anisotropy of silicon nonlinearities without too much computational effort.Since the maximum intensity is achieved at the center of a single-mode waveguide, the one-dimensional implementation also allows one to estimate the strongest impact of nonlinear optical phenomena on the field propagation.

Numerical examples and discussion
The FDTD algorithm developed here allows one to simulate the propagation of a short optical pulse through a silicon waveguide under quite general conditions.In this section we use the simplified FDTD algorithm to analyze polarization-dependent optical phenomena in silicon waveguides.We consider propagation of optical pulses through waveguides fabricated along the [110] direction on the (001) surface (ϑ = π/4).We include all nonlinear effects but neglect thermal effects.In the case of XPM, we assume that the two pulses at different wavelengths are launched at the same end of the waveguide.The parameter values are chosen to be: ε 2 = 1.72 × 10 −19 m 2 /V 2 , β TPA = 0.9 cm/GW, n 0 = 3.17, and α L = 1 dB/cm.

Input-output characteristics
It is well known that the output power in silicon waveguides saturates with incident power, due to the nonlinear absorption resulting from TPA and FCA [52,77,78].For relatively short optical pulses, TPA dominates FCA, and the output intensity exhibits dependence on input polarization.Figure 5 shows the input-output characteristics for three silicon waveguides, pumped by Gaussian pulses with a full width at half maximum (FWHM) of 1.4 ps polarized along the TE and TM directions, by plotting the output peak intensity as a function of the input peak intensity.It can be seen that the output intensity is higher for TM pulses because TPA is larger for the TE mode than for the TM mode by a factor of 1  2 (1 + ρ) ≈ 1.14.More generally, Eqs.(10b) and (10c) show that the anisotropy of TPA depends on the waveguide orientation specified by the angle ϑ as 1 + 1 2 (ρ − 1) sin 2 (2ϑ ).Since the effects of TPA and FCA accumulate with propagation, the peak output intensity is lower in longer waveguides.

Nonlinear polarization rotation
Aside from different rates of TPA, TE and TM modes in silicon waveguides exhibit different nonlinear phase shifts imposed by SPM and XPM.Both the TPA and XPM modify the state of polarization (SOP) of a sufficiently strong optical field [35].The efficiency of the nonlinear polarization rotation can be characterized by the transmittance of a linearly polarized input pulse through the analyzer that blocks the pulse for low input intensities [36].Figure 6 shows such a transmittance as a function of the input polarization angle ϕ for a 1.4-ps Gaussian pulse; ϕ = 0 and ϕ = π/2 correspond to TE and TM modes, respectively (see Fig. 3).We define the transmittance as a ratio of the energy flux, +∞ −∞ I(t) dt at the output and input ends of the waveguide.
One can see from Fig. 6 that the pulse preserves its initial SOP not only when it propagates as either a pure TE or TM mode but also when the input SOP corresponds to ϕ 0 = tan −1 1/ √ 2 ≈ 35 • .Interestingly, the value of ϕ 0 does not depend on the material parameters or pulse characteristics.For an arbitrary fabrication direction ϑ , it can be found by noting that conservation Thus, ϕ 0 reaches its maximum value of π/4 when the waveguide is fabricated along the principal axis [100] or [010].The relative nonlinear phase shift between the TE and TM polarization components increases monotonically with input pulse intensity I 0 and waveguide length L. However, due to nonlinear losses, there are optimal values of I 0 and L that maximize the pulse transmittance [36].

Polarization variations along the pulse
Since the processes of SPM and TPA depend on optical intensity, different parts of a pulse exhibit different polarization changes.It is useful to plot variations of the SOP along the output pulse on the Poincaré sphere.The instantaneous SOP of the pulse can be characterized by the polarization ellipse inscribed by the tip of the electric field vector [79].The time evolution of SOP at the waveguide output is then represented by a closed trace on the Poincaré sphere.Figure 7 shows how SOP changes for the most intense pulse in Fig. 6(a).Each trajectory on the Poincaré sphere starting from the equator corresponds to a specific linear SOP of the input pulse.Consider a pulse with ϕ = π/10 portrayed by the trace AB.As intensity builds along the leading edge of the pulse, anisotropic nonlinear effects come into force and start increasing the ellipticity and azimuth of the instantaneous SOP (arc AA ′ ).After the ellipticity and azimuth peak at the pulse center with maximum intensity (point B), the efficiencies of the TPA and SPM processes start to decrease (arc BB ′ ), and the SOP returns to the initial linear SOP at the trailing edge of the pulse (point A).
It is interesting to note that the arcs AA ′ B and BB ′ A on the Poincaré sphere nearly coincide.A detailed analysis shows that this feature reflects the fact that pulse polarization predominantly changes in response to cross-intensity modulation of the TM and TE modes caused by the degenerate four-wave mixing (FWM), rather than by the SPM-induced relative phase shift between the two modes.The relative nonlinear phase shift between the two polarization components is negligible for picosecond pulses in waveguides shorter than 1 mm, but it becomes important in the quasi-CW regimes or for longer waveguides.Visualization of the SOP on the Poincaré sphere also reveals the significance of polarization angle ϕ 0 in Eq. (22).When ϕ ϕ 0 , the intensity of the TE mode builds up more slowly than the intensity of the TM mode at the leading edge of the output pulse owing to cross-intensity modulation; this results in the development of the right-handed elliptic polarization.In contrast, when ϕ > ϕ 0 , energy transfer between the two polarization components results in steeper growth of the TM mode and in the development of the left-handed elliptic polarization.The TM and TE modes excited by an optical pulse whose linear SOP is characterized by the angle ϕ 0 do not exhibit degenerate FWM and, therefore, do not interchange any energy.

Nonlinear switching through polarization changes
Nonlinear polarization rotation (NPR) in silicon waveguides can be used for optical switching [34,36].Switching of a linearly polarized CW signal is realized by launching an intense pump into the waveguide, which changes the SOP of the CW signal and results in its partial transmittance through the output polarizer (called the analyzer).If the signal is in the form of a pulse, its intensity must be small enough to avoid nonlinear effects in the absence of the pump.Otherwise, different parts of the signal experience amounts of NPR, and the signal cannot be completely blocked by the analyzer.This restriction on signal intensity can, however, be lifted, if the signal is in the form of a CW beam.
To analyze the switching scenario in detail, we consider switching of a 200-THz CW beam (wavelength 1.5 by 350-fs Gaussian pulses whose spectrum peaks at 210 THz.The intensity of the CW beam is assumed to be 1 GW/cm 2 so that it produces almost no nonlinear effects.Figures 8(a)-8(c) show the analyzer transmittance or switching efficiency defined as the ratio of the output signal energy (per unit area near the mode center) to the energy of a 350-fs Gaussian pulse with a peak intensity of 1 GW/cm 2 .The analyzer is aligned perpendicular to the linear SOP of the input CW beam, since modal birefringence is absent and signal intensity is too low to cause significant NPR in the absence of the pump.It is evident that switching efficiency depends strongly on polarizations of the pump and signal.For example, when the pump excites either the TE or TM mode, the maximum switching efficiency is achieved for the signal that excites both of them equally (i.e., for ϕ = π/4); if the pump is initially polarized along the direction ϕ = π/4 (TE-TM case), maximum transmittance occurs for ϕ = 0 and ϕ = π/2 [see Fig. 8(a)].
Large values of signal transmittance in Fig. 8 are a consequence of using subpicosecond pump pulses and short waveguides.For a given free-carrier lifetime, the transmittance can be maximized by varying the waveguide length L and the pump intensity I p0 [36].The existence of optimal values for L and I p0 is evidenced by the maxima seen in Figs.8(b) and 8(c).Signal profiles at the output of 0.1-, 0.2-, and 0.4-mm-long waveguides are shown in Fig. 8(d).It is seen that the signal is similar in shape to the pump (although slightly asymmetric) in the case of short waveguides.The asymmetry is produced by FCA, which predominantly attenuates the trailing edge of the pump pulse and reduces the NPR that it causes.In longer waveguides, the signal profile changes dramatically.For example, the profile exhibits a local minimum near the point where the pump intensity peaks for L = 0.4 mm.This minimum in signal transmittance is due to excessive NPR leading to partial restoration of the initial signal SOP.

Polarization dependence of Raman amplification
As a final example, we use our FDTD simulator to illustrate how the SOPs of two pulses, separated in their carrier frequencies by the Raman shift, affect the efficiency of the SRS process that is responsible for Raman amplification.Both the signal and pump pulses are assumed to be 1.4-ps Gaussian pulses with peak intensities of 1 and 100 GW/cm 2 , respectively.Figure 9 shows the evolution of the peak intensity of the signal pulse along the waveguide length for different SOPs of the input pulses.As can be seen there, no SRS occurs, if both the pump and signal pulses propagate in the form of TM modes, a well-known result for silicon waveguides fabricated along the [110] direction [53,54,75].
The different evolutions of the signal peak intensity along the waveguide length for three other polarization combinations (the blue, green, and red curves) are due to the anisotropy of the TPA and SRS phenomena.The Raman amplification is considerably stronger when the pump and signal are launched with orthogonal SOPs, a result that does not appear to have been fully appreciated in the literature on silicon Raman amplifiers.It should also be noted that the reduction of peak signal intensity with increasing waveguide length is caused not only by the cumulative nature of FCA, but also by FWM and CARS [19,63].Calculation of the output signal intensity for other possible polarizations of the pump and signal is required to optimize the performance of silicon Raman amplifiers.

Conclusions
We have presented a comprehensive three-dimensional FDTD model for studying nonlinear optical phenomena in silicon waveguides, which allows-for first time to our knowledge-for anisotropy of the Kerr effect, two-photon absorption, and stimulated Raman scattering.Based on our model, we developed a computationally efficient FDTD algorithm suitable for simulating polarization-dependent propagation of optical pulses through silicon waveguides.We applied this algorithm to examine several polarization effects that are most favorable for applications.
In particular, we demonstrated that the developed algorithm can be used for the optimization of Kerr shutters and silicon Raman amplifiers.Wherever possible, we compared our findings with available experimental data and theoretical results obtained within the framework of the slowly varying envelope approximation.The comparison revealed a high accuracy of our FDTD simulator within its applicability domain.Owing to its generality in handling complex waveguide geometries and short optical pulses, our FDTD simulator is very useful for testing the suitability of anisotropic nonlinearities for different silicon-based photonic devices.Also, we envision using it as a testing vehicle for nonlinear differential equations that describe specific anisotropic phenomena in silicon waveguides.

Fig. 1 .
Fig. 1.Staggered space arrangement (left) and leapfrog time ordering (right) of discretized electromagnetic-field components using a cubic Yee cell.

ξ
by averaging their values over the four nearest nodes surrounding P (a)

Fig. 4 .
Fig. 4. Three one-dimensional Yee cells corresponding to the simplified FDTD model.The electric and magnetic fields along the propagation direction α are ignored.

Fig. 6 .Fig. 7 .
Fig.6.Efficiency of SPM-induced polarization rotation for different input SOPs of a 1.4-ps Gaussian pulse.The output polarizer is perpendicular to the input polarization state of the pulse.We choose τ c = 0.8 ns; other parameters are given in the text.

Fig. 8 .
Fig. 8. Switching efficiency of a CW beam for a 350-fs Gaussian pump pulse as a function of (a) input polarization angle ϕ of the CW beam, (b) waveguide length L, and (c) pump's peak intensity I p0 .Panel (d) shows switching windows for three different waveguide lengths; dashed curve shows the Gaussian pulse that is used in the definition of transmittance.In panels (b)-(d), ϕ = π/4.In all panels input CW intensity is 1 GW/cm 2 and τ c = 0.8 ns; other parameters are given in the text.

#Fig. 9 .
Fig. 9. Peak intensity of TM-and TE-polarized signals amplified via SRS by TM-and TE-polarized pumps.Both input signal and input pump are assumed to be 1.4-ps Gaussian pulses with peak intensities of 1 and 100 GW/cm 2 , respectively; carrier frequencies are 200 and 215.6 THz; g R = 76 cm/GW; other parameters are given in the text.

Table 1 .
Values of d µ specifying relative positions of components V ξ inside a fourdimensional unit cell.

Table 1 )
; the other two components of the electric field are defined at the adjacent nodes.It is common to calculate the values of these two components at the position of P to incorporate the Kerr effect and TPA into the FDTD model.
#133336 -$15.00USD Received 12 Aug 2010; revised 22 Sep 2010; accepted 22 Sep 2010; published 23 Sep 2010 (C) 2010 OSA 27 September 2010 / Vol. 18, No. 20 / OPTICS EXPRESS 21437 because they are common to all functions.We continue to use this notational simplification in what follows.To incorporate linear dispersion into the preceding we calculate the value of P LD at the time step n + 1, by independently updating its s components that correspond to different terms in the Sellmeier equation, i.e., #133336 -$15.00USD Received 12 Aug 2010; revised 22 Sep 2010; accepted 22 Sep 2010; published 23 Sep 2010 (C) 2010 OSA 27 September 2010 / Vol. 18, No. 20 / OPTICS EXPRESS 21440