1. 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–9]. A relatively high refractive-index 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].

The strong optical nonlinearities of silicon give rise to a variety of nonlinear effects, such as stimulated Raman scattering (SRS) [10, 11], self-phase modulation (SPM) [12, 13], cross-phase modulation (XPM) [14], cross-absorption modulation (XAM) [15], four-wave mixing (FWM) [16–18], and coherent anti-Stokes Raman scattering (CARS) [9, 19, 20]. These effects have been successfully utilized in a number of chip-scale photonic devices, including lasers [21, 22], amplifiers [23, 24], switchers [15, 25], modulators [26–28], broad-band frequency converters [20, 29], detectors [30], all-optical logic gates [25], continuum generators [31, 32], and pulse compressors [33]. Since the Kerr effect, two-photon absorption (TPA), and Raman scattering are anisotropic in the case of silicon, the state of polarization of the optical field and waveguide orientation could be tailored to take advantage of optical anisotropy for making novel optical devices. Even though there is a strong potential for anisotropic effects to be used in new devices, none of the above-listed devices utilizes optical anisotropy. Meanwhile, recent theoretical works have demonstrated that the anisotropy exhibited by silicon waveguides can be used to realize fast optical switching [34–36] and power equalization [36].

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 electromagnetic 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.

2. 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]

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 μ ₀ and ε ₀ 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

where $\widehat{\alpha }$ , $\widehat{\beta }$, and $\widehat{\gamma }$ 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} ∈ ℕ, 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.

 

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

Download Full Size | PDF

Table 1. Values of d_μ specifying relative positions of components V_ξ inside a four-dimensional unit cell.

View Table

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 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.,

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–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_{\xi }^{(i)}(\dots )$ denotes the functional dependence.

In the case of the anisotropic material response, the value of $P_{\xi }^{(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. However, only the component E_ξ coincides in space with $P_{\xi }^{(a)}$ (see 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_{\xi }^{(a)}$ by averaging their values over the four nearest nodes surrounding $P_{\xi }^{(a)}$ [43–45]. For example, the components E_β and E_γ are projected onto the position of the component $P_{\alpha }^{(a)}$ using the relations (see Fig. 2)

 

Fig. 2 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 (i ± 1/2, j ± 1/2, k) and (i ± 1/2, j, k ± 1/2), respectively.

Download Full Size | PDF

Equations (3)–(5) are solved iteratively for a fixed value of D until the desired accuracy is achieved.

3. 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⁻³) can arise in SOI waveguides with asymmetric cross sections [46], even though the intrinsic optical anisotropy of silicon is relatively low (∼ 10⁻⁶) [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, $\tilde{\chi }$ ⁽¹⁾(ω).

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–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–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, ${\tilde{\chi }}_{\kappa \lambda \mu \nu }^{\text{e}}$ and ${\tilde{\chi }}_{\kappa \lambda \mu \nu }^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–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, $\tilde{\chi }$ ^FC(|E|⁴, ω).

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 ${\tilde{\chi }}_{\kappa \lambda \mu \nu }^{(3)}={\tilde{\chi }}_{\kappa \lambda \mu \nu }^{\text{e}}+{\tilde{\chi }}_{\kappa \lambda \mu \nu }^R$, ω ₃ = ω − ω ₁ − ω ₂, 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.

 

Fig. 3 Relative orientation of FDTD (α, β, γ) and crystallographic (x, y, z) axes adopted in the paper. For $\vartheta =\frac{\pi }{4}$, the FDTD axes α, β, and γ coincide, respectively, with the [110], [1̄10], and [001] crystallographic directions. The inset shows the TM and TE polarizations and an arbitrary linear polarization determined by the angle φ.

Download Full Size | PDF

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].

3.1. 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(\omega )\approx \sqrt{\tilde{\varepsilon }\prime (\omega )}$, while the imaginary part is responsible for linear losses, $\tilde{\epsilon }$″(ω) ≈ cn ₀ α_L/(2ω), where α_L is the linear-loss coefficient and n ₀ ≡ n(ω ₀) is the refractive index at the absolute maximum of the input field spectrum located at frequency ω ₀. 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 ²(ω) is given by the 2s-pole Sellmeier equation [62–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 ₁ = 9.733, a ₂ = 0.936, ω ₁ = 1032.49 THz, and ω ₂ = 817.28 THz. These values provide an accurate fit of the refractive-index data for silicon in the wavelength range from 1.2 to 2 μm [65].

3.2. 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 ₂ = Δn _NL/I ≈ 6 × 10⁻⁵ cm²/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 two-photon 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–27, 33].

The general form of the susceptibility ${\tilde{\chi }}_{\kappa \lambda \mu \nu }^{\text{e}}(\omega ;{\omega }_1,{\omega }_2,{\omega }_3)$ required to calculate the electronic part of the nonlinear polarization, caused by an optical beam with an arbitrary spectrum, is not known. Fortunately, the response time of bound electrons is comparable to the duration of an optical cycle and can be considered to be negligible for a relatively narrow field spectrum [49, 55]. In this case, dispersion of the electronic susceptibility can be ignored, and the Kleinman's symmetry can be used to factorize the third-order susceptibility as [49]

where ${\tilde{\chi }}_{\alpha \alpha \alpha \alpha }^{\text{e}}(\omega )={\varepsilon }_2-{\eta }_{\text{TPA}}/(i\omega )$, ${\varepsilon }_2={\varepsilon }_0{n}_0^{2}c{n}_2$, and ${\eta }_{\text{TPA}}=\frac{1}{2}{\varepsilon }_0{c}^2{n}_0^2{\beta }_{\text{TPA}}$. The real and imaginary parts of ${\tilde{\chi }}_{\alpha \alpha \alpha \alpha }^{\text{e}}(\omega )$ represent the Kerr effect and TPA, respectively. All anisotropic effects are included through the anisotropy tensor ℰ_κλμν.

In the crystallographic coordinate system, the anisotropy tensor ℰ_klmn has the following well-known 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 ℰ_klmn according to Eqs. (7) and (8), we obtain the following nonzero elements of the anisotropy tensor in the FDTD coordinates:

where $A_{\vartheta }=\frac{1}{2}(\rho -1){sin}^2(2\vartheta )$ and $B_{\vartheta }=\frac{1}{4}(\rho -1)sin(4\vartheta )$.

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_{\kappa }^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) to incorporate the Kerr effect and TPA into the FDTD model.

3.3. 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–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 ℛ_κλμν describes the anisotropy of the Raman scattering and the function H̃(ω) represents the Raman gain profile,

Here, ξ_R = 2ε ₀ n ₀ c ² 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 (ℛ_kjjj = δ_kj + δ_kj − 2δ_kj = 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 ℛ_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ϑ) and $D_{\vartheta }=\frac{1}{2}sin(4\vartheta )$.

Using Eq. (12) and introducing a new function H (t) = F ⁻¹ [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

3.4. 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(\mathbf{\text{r}},t)=\frac{1}{2}c{\varepsilon }_0{n}_0{\left|\mathbf{\text{E}}\left(\mathbf{\text{r}},t\right)\right|}^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⁻⁴ K⁻¹, 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³ 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.

4. 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.

4.1. 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 (i + d ₁,J + d ₂,k + d ₃) because they are common to all functions. We continue to use this notational simplification in what follows.

To incorporate linear dispersion into the preceding equation, we calculate the value of $P_{\kappa }^{\text{LD}}$ at the time step n + 1, by independently updating its s components that correspond to different terms in the Sellmeier equation, i.e.,

Conversion into the time domain of the equation for the jth component, $({\omega }_j^2-{\omega }^2){\tilde{P}}_{\kappa j}^{\text{LD}}={\varepsilon }_0{a}_j{\omega }_j{\tilde{E}}_{\kappa }$, 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 ${\mathrm{\Xi }}_{\kappa \lambda }^R$ in Eq. (14) is governed by the differential equation

whose discretized version takes the form

The relations between the functions $S_{\kappa \lambda }^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+\frac{1}{2}$ can be found using its value at the time step $n-\frac{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:

where $D_{\kappa }^{\text{NL}}=D_{\kappa }-P_{\kappa }^{\text{LD}}$ and

Notice that each of the auxiliary functions $S_{\kappa }^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_{\kappa j}^{\text{LD}},{\mathrm{\Xi }}_{\kappa \lambda }^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.

4.2. 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_{\alpha }^R$ is not produced by the transverse components E_β and E_γ when the waveguide is fabricated along the crystallographic direction determined by the angle ${\vartheta }_j=\frac{\pi }{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 ₀ 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 ₀).

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 ${ℬ}_{\kappa }^{\pm }$ 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.

 

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

Download Full Size | PDF

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.

5. 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: ε ₂ = 1.72 × 10⁻¹⁹ m²/V², β _TPA= 0.9 cm/GW, n ₀ = 3.17, and α_L = 1 dB/cm.

5.1. 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 $\frac{1}{2}(1+\rho )\approx 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+\frac{1}{2}(\rho -1)$ sin²(2ϑ). Since the effects of TPA and FCA accumulate with propagation, the peak output intensity is lower in longer waveguides.

 

Fig. 5 Input–output characteristics for 100-, 200-, and 500-μm-long silicon waveguides pumped by TM- and TE-polarized, 1.4-ps Gaussian pulses. Other parameter values are given in the text.

Download Full Size | PDF

5.2. 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, ${\int }_{-\infty }^{+\infty }I(t)\text{d}t$ at the output and input ends of the waveguide.

 

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.

Download Full Size | PDF

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 ${\phi }_0={tan}^{-1}\left(1/\sqrt{2}\right)\approx 35°$. Interestingly, the value of φ ₀ does not depend on the material parameters or pulse characteristics. For an arbitrary fabrication direction ϑ, it can be found by noting that conservation of SOP requires $S_{\beta }^K=S_{\gamma }^K$ at any point along the waveguide. Using Eqs. (10b) and (10c), this condition yields

Thus, φ ₀ 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 ₀ and waveguide length L. However, due to nonlinear losses, there are optimal values of I ₀ and L that maximize the pulse transmittance [36].

5.3. 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).

 

Fig. 7 Poincaré-sphere representation of SOP variations along a 1.4-ps Gaussian pulse at the output of a 0.5-mm-long silicon waveguide. Numbers near the traces staring from the equator show φ values for the input SOPs (see Fig. 3); φ ₀ ≈ 35°. Traces starting from the poles correspond to circularly polarized input pulses. Here, I ₀ = 400 GW/cm²; other parameters are the same as in Fig. 6(a).

Download Full Size | PDF

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 φ ₀ in Eq. (22). When φ < φ ₀, 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 φ > φ ₀, 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 φ ₀ do not exhibit degenerate FWM and, therefore, do not interchange any energy.

5.4. 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 different 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 μm) by 350-fs Gaussian pulses whose spectrum peaks at 210 THz. The intensity of the CW beam is assumed to be 1 GW/cm² 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². 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)].

 

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_{p 0}. 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² and τ_c = 0.8 ns; other parameters are given in the text.

Download Full Size | PDF

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_{p 0} [36]. The existence of optimal values for L and I_{p 0} 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.

5.5. 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², 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].

 

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², respectively; carrier frequencies are 200 and 215.6 THz; g_R = 76 cm/GW; other parameters are given in the text.

Download Full Size | PDF

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.

6. Conclusions

We have presented a comprehensive three-dimensional FDTD model for studying nonlinear optical phenomena in silicon waveguides, which allows—for the 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.