Analysis of material susceptibility in silicon on insulator waveguides with combined simulation of four-wave mixing and linear mode coupling

We derive propagation equations modeling third-order susceptibility-induced nonlinear interaction and linear mode coupling in waveguides. We model material susceptibility with Raman and electronic response which include approximations suited for optical communications. We validate our model by comparing numerical integration of the propagation equations to continuous wave measurements of a silicon on insulator waveguide.


Introduction
All-optical signal processing is a promising technique for various applications. Two prominent examples are optical phase conjugation (OPC) for mitigating nonlinear impairments in optical fiber links (Vukovic et al. 2015) and wavelength conversion (WLC) which shows high potential for several scenarios. In a narrowband sense, WLC can enable all-optical routing in datacenters or metropolitan networks. In a broadband sense, WLC allows to jointly multiplex several fully-loaded C-bands into one fiber by shifting them to distinct bands. Different approaches based on highly non-linear fibers (Rademacher et al. 2020), LiNbO 3 waveguides (Minzioni et al. 2010), and silicon on insulator (SOI) waveguides (Gajda et al. 2018;Ronniger et al. 2020) have been investigated. We present a simulation for the latter, including nonlinear effects and linear mode coupling, with an accurate model of the susceptibility in silicon.
Third-order material nonlinearity-the origin of four-wave mixing (FWM)-is the basis for all-optical signal processing. Due to the high nonlinear potential of silicon, the application of an SOI waveguide demands a more detailed consideration of nonlinearity than it is provided by the optical nonlinear Schrödinger equation. Therefore, we derive a set of coupled differential equations that model linear and nonlinear interaction based on the third-order susceptibility (resembling pulse propagation in Agrawal (2013)). As a result of the high Raman gain coefficient in silicon-which is three to four orders of magnitude higher compared to silica (Jalali et al. 2006)-it is essential to include the effect of molecular vibrations (Raman response) besides the commonly considered electron vibrations (Kerr effect).
The input-output conversion efficiency (CE) defined as idler output power over signal input power is a critical measure of the nonlinear process's performance. For achieving a strong idler build-up, phase matching (PM) has to be performed and multi-mode operation is usually preferred since an additional degree of freedom is obtained. We refer to Kernetzky et al. (2020) for a detailed analysis of geometry optimization for PM in SOI waveguides. In Höfler et al. (2021) we described a model for the susceptibility calculation in silicon. This is extended in this work by applying the susceptibility model to the propagation equations of an SOI waveguide, as well as by comparing the results to a continuous wave measurement of the waveguide.

Modeling light propagation
For investigating the spatial development of the amplitudes of the interacting frequencies along the waveguide, the corresponding differential equations will be stated, starting from the wave equation The electric displacement field is defined as where the extended material permittivity matrix models material dispersion by r (x, y) = n 2 (x, y), linear perturbations by the matrix δ r (x, y, z), light attenuation by r (x, y), and where I is the identity matrix. Rearranging all perturbations into one vector yields with P = δ r E − j r E + P nl . Inserting D in Eq. 1 gives We model the nonlinear polarization vector as which contains the susceptibility tensor ↔ χ [3] ∈ C 3×3×3×3 , and where . . . represents the tensor product. Stating Eq. 6 in sum-notation gives with i, j, k, l ∈ {x, y, z}. The unmodulated total propagating electrical field can be written as a superposition of modes (m) at discrete positive frequencies , and propagation constants β (m) f i . In the following, we make the common assumption that P causes a z-dependence ofÊ, but does not have any effect on . Since the outer gradient in term A in Eq. 5 alters the transversal field profiles, ∇ · P = 0, ∀ P has to hold to fulfill the assumption that P does not affect . With that, evaluating Eq. 5 at the positive frequency f 0 (analytic signal) leads to with P := P +f 0 and E:= E| +f 0 . After quite some calculus, algebra and by assuming ∂ z r = 0, the x component of Eq. 9 becomes with β 0 = 2πf 0 /c 0 . Since the field amplitudeÊ (m) f 0 (z) only changes due to perturbations P , all addends with ∂ zÊ vanish in the homogeneous equation ( P = 0). This indicates that due to mode orthogonality all m summands needs to be equal to zero, i.e., B = 0 ∀ m. Based on the assumption that P does not alter , this also has to hold for P = 0. Consequently, by extending Eq. 10 by the y and z component one gets which also makes use of the slowly varying wave approximation ∂ 2 The nonlinear part of P is governed by the third-order material susceptibility ↔ χ [3] . This tensor generates a nonlinear perturbation at frequency f 0 by combining three interacting light waves. Inserting Eq. 8 into Eq. 6 and only considering frequency combinations (f ζ , f η , f ρ ) with two positive and one negative contribution (thus, considering processes such as OPC and Bragg scattering (BS)) gives where ↔ X [3] is the Fourier transform of ↔ χ [3] , i.e., F ↔ χ [3] = ↔ X [3] .
Inserting Eq. 12 in Eq. 11, multiplying with (a) * f 0 from the left and integrating over the cross section (making use of mode orthogonality), we obtain the propagation equation for mode (a) at frequency f 0 with The mode coupling (MC) coefficient C (m) couples modes at the same frequency whereas the nonlinearity coefficient N (a) (m 1 ,m 2 ,m 3 ) couples between modes at all possible frequency combinations. This leads to a set of coupled differential equations of the type of Eq. 13 which has to be solved numerically. While the source fields can propagate in different modes and at arbitrary frequencies, the efficiency of the linear and nonlinear processes is determined by energy conservation and by the phase mismatches β lin and β , respectively.

Susceptibility in silicon
In the following, the third-order susceptibility for silicon as origin of nonlinear effects is analyzed. In silicon, there are two main parts that account for nonlinear processes, namely the electronic contribution (e) due to bound electrons and the Raman contribution (R) stemming from atomic lattice vibrations. Hence, it is reasonable to split the third-order nonlinear susceptibility tensor into its main parts, i.e., ↔ X [3] = ↔ X e + ↔ X R , and investigate each part separately . Throughout this section, we assume a waveguide fabricated on a (001) surface and parallel to the [110] direction.

Electronic susceptibility in frequency domain
Due to spatial symmetry, only 21 out of 81 entries of ↔ X e are nonzero, of which four are independent of each other (Boyd 2003). If only wavelengths λ > λ min = 1.10 μm are considered, the Kleinmann condition is satisfied and three of the four independent elements can be approximated to be equal Bristow et al. 2007;Osgood et al. 2009). Usually the electronic susceptibility is considered as nearly constant, since variations of the interacting frequencies lead to only small fluctuations of ↔ X e . For a wavelength range including the commonly used optical bands from O to L, the spacing between interacting frequencies is small enough to treat nonlinearity caused by the electronic contribution as being independent of the interacting frequencies, i.e., With this assumption and for λ ∈ [1.2 μm, 2.4 μm] (also including bands O to L), the last two independent entries of ↔ X e can be related to each other as well (Zhang et al. 2007). Furthermore, the real and imaginary part of the last independent entry is linked to the Kerr coefficient n 2 and the two-photon absorption coefficient β tpa , respectively. Altogether, this leads to the 21 nonzero elements with (Osgood et al. 2009;Hon et al. 2011), where the characteristics of n 2 and β tpa can be found in Bristow et al. (2007).

Raman susceptibility in frequency domain
The Raman contribution ↔ X R emerges from the interaction of light with lattice vibrations (phonons) of the material. If the difference of two incident light waves coincides with the frequency of the lattice vibration (resonance), the atom is excited to a higher vibrational eigenstate. The susceptibility elements induced by Raman scattering can be stated as Dim-itropoulos et al. (2003); Jalali et al. (2006); Lin et al. (2007) with and i, j, k, l ∈ {x, y, z}, where = 105 GHz is the FWHM-bandwidth, f v = 15.6 THz the vibrational eigenstate frequency, Z 0 = μ 0 0 , g R the Raman gain coefficient, and R n , n ∈ {1, 2, 3} the three Raman matrices with Each Raman matrix corresponds to the respective displacement of the phonons along the crystallographic directions of the medium and reflects its crystal symmetry. The terms 3 n=1 (R i j ) n · (R kl ) n and 3 n=1 (R il ) n · (R jk ) n determine the 18 nonzero elements of ↔ X R . All denoted frequency values apply only at room temperature, i.e., T ≈ 300 K. The formula and the required parameters for the approximation of g R can be found in Jalali et al. (2006); ; Ralston and Chang (1970); Renucci et al. (1975).
As indicated by Eq. (16), the entries of the Raman susceptibility consist of two contributions, X 1 ijkl (f 2 − f 1 ) and X 1 ijkl (f 2 − f 3 ). This can be derived from the fact that there exist two potential ways to promote atoms from the ground state to a higher vibrational eigenstate. Figure 1 illustrates the Raman process f 0 = f 1 − f 2 + f 3 in the special case of resonance. A photon of energy hf 1,3 excites the atom from the ground state to a virtual energy state E 1 . Then, stimulated by a photon of frequency f 2 , part of the energy is used to promote the atom to the vibrational eigenstate E v = hf v , the other part is emitted as a photon with frequency f 2 . Thus, photons with energy hf 1,3 can be understood as the driving force for providing atoms to the vibrational eigenstate E v . If a photon of energy hf 3,1 is absorbed by an atom located at the vibrational eigenstate E v , the atom is excited to the virtual state E 2 and falls back to the ground state immediately while emitting a photon of frequency f 0 . Considering f 1 as the frequency that implicitly provides atoms to the vibrational eigenstate, X i jkl 1 is differently pronounced depending on how precisely the resonance frequency f v is hit by the frequency difference f 1 − f 2 . In the case of resonance, i.e., f 1 − f 2 = f v , the maximal value of X i jkl 1 is obtained. Analogously, X i jkl 2 considers the possibility of resonance between f 2 and f 3 , by regarding f 3 as the frequency that implicitly supplies atoms to the vibrational eigenstate.

Numerical impelementation
In our simulation, we launch three waves (two pumps and one signal) with discrete frequencies into an SOI waveguide. By linear MC, waves at all frequencies couple into all  Table 1 The considered set of coupled differential equations which has to be solved simultaneously available modes at the same frequencies. Due to the nonlinear interaction, waves at all possible frequency combinations are generated in all available modes. Both interactions occur differently pronounced depending on phase mismatches β lin and β and coefficients C and N . Consequently, the number of possible combinations for interaction quickly rises. We limit the simulation by only considering light propagation in the guided modes of the waveguide (TE 0 , TE 1 , TE 2 ) and by only taking frequencies into account that will be generated with non-negligible efficiency. Thus, the frequencies of the three input waves, i.e., f P 1 , f P 2 , f S , and the two generated frequencies f BS = −f P 1 + f S + f P 2 and f OPC = +f P 1 − f S + f P 2 (arising from the BS and OPC FWM processes) are considered. Although this limitation drastically reduces complexity, still 15 coupled differential equations have to be solved simultaneously as summarized in Table 1.

The transversal field distributions
(m) f i and their corresponding propagation constants β (m) f i for the considered waveguide were calculated with an own full-vectorial mode solver based on Fallahkhair et al. (2008). After determining the required ↔ X [3] (f ζ , f η , f ρ ) (for all possible frequency combinations accounting to f i ), the nonlinear coupling coefficients N (a) (m 1 ,m 2 ,m 3 ) can be calculated. The influence of linear coupling and thus the magnitude of the linear coupling coefficients C (a) (m) is determined by the coupling matrix δ r which models deviations from the unperturbed refractive index profile, e.g., caused by waveguide roughness or mechanical stress. Since δ r depends on various factors it is difficult to find analytic expressions for its entries. Therefore, a heuristic approach is often preferred that uses δ r to match the simulation and measurement results of linear coupling.

Simulation (I)
First, an SOI waveguide with rib width 1800 nm, slab height 100 nm, SOI height 220 nm, propagation length 2 cm and one dip of width 400 nm and depth 70 nm is considered. Figure  2 shows the waveguide's geometry, which is the same as in Kernetzky et al. (2020). Table 2 summarizes the input/output wavelengths and their corresponding modes, optimized for best efficiency of the OPC process. We select realistic values for input powers, and mode-dependent attenuation coefficients, which are also included in Table 2. The resulting power evolution along the waveguide-without linear MC (C = 0)-is shown in Fig. 3a. It can be seen, that the propagation of the input waves (marked as , and ) is dominated by the linear waveguide attenuation; the power decrease caused by the nonlinear power transfer to the OPC and BS idler is hardly visible in the log-domain. The OPC process ( ) with best PM is created with highest efficiency, while the BS process ( ) with slightly worse PM is created with less efficiency. The additional decrease in BS idler power compared to the OPC idler is caused by the larger phase mismatch β . All other potential FWM products at other modes and frequencies are not shown, since they are highly phase-mismatched, and thus have negligible power.
For the analysis of MC, we heuristically set   to introduce MC. This breaks the orthogonality of modes in the computation of C in Eq. 14, by not taking into account the small but non-negligible y and z components of the transversal field profiles. We emphasize again, that more meaningful entries in this matrix need to be found by matching simulation with experimental data. Figure 3b shows the power evolution with additional linear MC. In principle, all waves couple into all considered modes. We exemplarily show one pump and one idler ( , ). The oscillation periods depend on the corresponding β lin values as L Oscillation = 2π | β lin | and the coupled powers on both, C and β lin . For instance, the magnitude of the normalized coupling coefficient of the OPC idler from TE 1 to TE 0 ( ) is C The corresponding phase mismatches | β lin | are |β and |β 2 × 10 5 m −1 ( ), respectively. While the coupling coefficient C from TE 1 to TE 2 is greater than from TE 1 to TE 0 (large values induce more coupled power), the phase mismatch β lin from TE 1 to TE 2 is also greater than from TE 1 to TE 0 (large values induce less coupled power). This leads to similar power levels of the OPC idler in TE 0 and TE 2 ( , ). Figure 3c shows a closeup of I OPC in mode TE 0 and reveals two oscillations. The slower one with period of roughly the plotted range is caused by direct coupling from TE 1 . The fast oscillation arises from second-order linear coupling from TE 2 and is therefore only weakly pronounced.

Simulation (II) with experimental verification
For verifying the simulation results we use measured data from an SOI waveguide with rib width 1672 nm, slab height 100 nm, SOI height 220 nm, and propagation length 11.3 mm ; Kernetzky et al. (2021)). Table 3 summarizes the simulation parameters optimized for best efficiency of the BS process. Since mode multiplexer and grating coupler at the input and output of the waveguide are lossy, the simulation parameters have to be adjusted accordingly which is reflected in a reduced input power, i.e., P P 1 = 13.92 dBm, P P 2 = 18.75 dBm, and P S = 6.28 dBm. A suitable means to evaluate the generation process of the BS idler is the input-output conversion efficiency (CE): CE = P I BS /P S (or CE = P I BS − P S in log-domain), where the signal and idler power is defined at the input (before mode multiplexer and grating coupler) and at the output (after mode demultiplexer and grating coupler) of the chip, respectively. Figure 4a shows that the simulated idler power is P I BS = −22.5 dBm at the end of the waveguide. This results in a CE of CE sim = −22.5 dBm − 5 dB − 11.28 dBm ≈ −39 dB, where -5 dB includes the losses caused by demultiplexer and grating coupler at the output. From Fig. 4b, the measured CE for λ I BS = 1303.44 nm is CE meas ≈ −42 dB. This is a very good match between measurement and simulation, since the remaining CE difference of 3 dB is very small considering measurement imperfections, material parameter, temperature variations, etc.

Conclusions
We modeled electronic and molecular parts ↔ X e and ↔ X R of the silicon susceptibility ↔ X [3] . We linked the components of ↔ X e to material parameters and presented a closed-form solution for ↔ X R . Although the approximations we applied limit the usable wavelength region, they are valid for commonly used optical transmission bands.
We derived and numerically integrated a set of coupled differential equations, which model linear and nonlinear light evolution along waveguides. Each frequency component in each mode is modeled by an equation that includes linear attenuation, linear MC, and FWM nonlinearity caused by the material susceptibility ↔ X [3] . We verified the validity of the proposed model by measured data of an SOI waveguide.