1. Introduction

In the last decades, silicon photonics has proved to be a suitable platform for high bandwidth and low cost on-chip communications and computing. Building blocks enabling light modulation are of paramount importance as they affect both the overall speed of communications and power consumption. In this context, silicon-based electro-optic modulators are actively investigated as they present the strong advantage of being monolithically integrated at low cost on silicon chip using CMOS-compatible processes.

However, due to the centrosymmetric nature of silicon, the refractive index change resulting from an applied electric field is extremely weak. Several methods are investigated to enhance electro-optical effects in silicon. For example, strained silicon waveguide allow breaking its inversion symmetry, hence enabling effective index change of Δn_eff ~10⁻⁵ for an applied voltage of 30 V [1]. Another technique consists in dynamically modifying the complex refractive index of silicon through electrical injection or removal of charge carriers. Using the carrier depletion effect in silicon, a modulation rate as high as 25 GHz was demonstrated [2]. However, high modulation speed requires high operating voltage and the modulation rate is limited by carrier recombination in silicon.

Alternatively, the monolithic integration of heterogeneous electro-optic materials on silicon appears promising. A prototypical example of such material is LiNbO₃ whose linear electro-optical effect (i.e. the Pockels effect) is already exploited to modulate light in fiber-optic communication systems. Other perovskite ferroelectric thin films, such as Pb(Zr_x,Ti_1-x)O₃ (PZT), BaTiO₃ (BTO) and KNbO₃ attract considerable interest due to their larger Pockels coefficients. In particular, bulk crystal BTO has the largest Pockels coefficient (r₄₂ ~800 pm V⁻¹) [3]. Promising results were obtained on BTO-based modulators: a bandwidth of 15 GHz was measured using a hybrid MgO/BTO/Si₃N₄ device [4]. Furthermore, recent advances in molecular beam epitaxy enable direct growth of monocrystalline BTO on silicon substrate through the use of SrTiO₃ (STO) templates [5]. Thanks to this technique, BTO thin films grown on silicon showing effective Pockels coefficients as high as r_eff = 148 pmV⁻¹ [6] and r_eff = 213 ± 49 pmV⁻¹ [7] were experimentally demonstrated.

To optimize the design of ferroelectric oxide-based electro-optic modulator, accurate modeling of such systems is required. Perovskite ferroelectric oxides are anisotropic crystals and therefore present both an ordinary and an extraordinary refractive indices. In such materials, the transformation of the index ellipsoid induced by the Pockels effect depends on the electric field direction respective to the crystal orientation. Furthermore, the deformation of the index ellipsoid will, in turn, affect the propagating modes and confinement in the waveguide. Therefore, an accurate model of such system needs to couple calculation of optical modes with the electrical response of such materials.

In the present work we develop a multi-physics model which couples an anisotropic full-vectorial finite-difference mode solver with a Laplace radio-frequency (RF) solver. It enables accurate calculations of refractive index modulation, propagation constant modification and electro-optical response induced by Pockels effect in anisotropic materials presenting nondiagonal permittivity tensor variation. To highlight the possibilities of this hybridized model, we use a BTO-based slot-waveguide as a testbed device to design and optimize an efficient electro-optic modulator (shown in Fig. 1). We further exploit the model to design index-matched travelling-wave electrodes that maximize the theoretical modulation bandwidth.

 

Fig. 1 Cross-section of the active zone of a BTO-based slot-waveguide electro-optical Mach-Zehnder modulator on a standard SOI substrate.

Download Full Size | PDF

2. Electro-optical effect in BaTiO₃

The Pockels effect is the linear variation of a material’s dielectric permittivity as a function of an applied electric field. This electro-optic effect only exists in noncentrosymmetric materials, i.e. materials whose structure do not possess inversion symmetry. An example of such material is BTO, which is a perovskite ferroelectric material with a 4mm symmetry. Its tetragonal structure has lattice parameters of a = 3.994 Å and c = 4.0335Å. In its bulk form, values of the first order electro-optic tensor elements (Pockels tensor elements) are r₁₃ = r₂₃ = 10 pm·V⁻¹, r₃₃ = 40 pm·V⁻¹, r₄₂ = r₅₁ = 820 pm·V⁻¹, and other r_ij = 0 [8]. When BTO is grown as a thin film, two different orientations can coexist: out-of-plane (c-axis) and in-plane (a-axis) oriented crystals. As the c/a ratio depends on numerous factors (growth techniques, conditions, nature of the substrate and thickness of the layer), for simplicity we restrict ourselves here to the c-axis case. Note that high quality fully c-axis-oriented monocrystalline BTO layers up to 20 nm can readily be grown on Si (001) substrate buffered by 5 nm-thick SrTiO₃ layer using molecular beam epitaxy (MBE) [9]. For the c-axis oriented thin films and taken into account the cross-section and axes shown in Fig. 1, the variation of the impermeability tensor in BTO is [10]:

We set$E_x^e=0$ here, as x is the propagation direction in the straight waveguide shown in Fig. 1. After some algebraic manipulation, the relative permittivity tensor in the BTO layer which is subjected to the applied electric field $E^e$is:

Interestingly, if we make a first order approximation on the permittivity tensor coefficients, we find the same formula as the one often used in literature [11]:

3. Optical and radio-frequency modeling

3.1 Anisotropic full-vectorial mode solver

Assuming a time-dependence as e^-jωt and a spatial x dependence of e^jßx, the Maxwell’s curl equations can be reduced to a system of two coupled eigenvalue equations for the transverse magnetic and electric fields components:

Note that in the following, the exact values of tensor components ε_xx, ε_yy, ε_zz, ε_yz, ε_zy, as displayed in Eq. (2), are used as inputs in Eqs. (4) and (5).

The coupling system of partial differential equations Eqs. (4) and (5) can be expressed as a sparse eigen system:

To numerically solve for the propagation constant, conventional mode solvers described in the literature only use the transverse magnetic field formulation to eliminate spurious modes and simplify the boundary conditions at interfaces [12]. A graded index approximation is sometimes necessary to mitigate the discontinuities at interfaces. Here, we use a finite difference method based on the 2nd order Taylor expansion with a nine-points-discretization scheme. This technique enables an accurate description of the discontinuities of the electric field at the waveguide’s corners without the need for a graded index approximation. This is particularly important in this study to precisely calculate the E_z component in a slot waveguide configuration. Thanks to that method, the same accuracy on the determination of the propagation constant is obtained for both electric and magnetic fields formulations. Validation of the straight mode solver was made with bench mark Si waveguide [13]. Applying the graded index approximation to our model, the obtained Δn_eff values are similar to that obtained by commercial tools (accuracy within 10⁻⁸). But when confronting modeling and experimental results, in the case of high-index-contrast silicon-based waveguides, a good agreement is only obtained when no graded index approximation is used in our model [14,15].

Transparent boundary condition for the edges of the calculation window is chosen to calculate the propagation loss [16].

The precise calculation of four key components of the electromagnetic field (E_y, E_z, H_y, and H_z) is then used to calculate the confinement factor in the active BTO layer (Γ) as:

3.2 Electro-optical and microwave models

We now want to calculate the change in the optical mode and propagation constant when the permittivity of the material is modified by the Pockels effect. According to the variation theorem [17], determination of the change Δβ in the propagation constant caused by a perturbation Δε is given by:

Because of the geometry of the electrode, the $E_z^e$ component of the electric field can be neglected (this assumption will be verified in section 4) and $\Delta \epsilon = -{\epsilon }_o{\epsilon }_e{r}_{42}{E}_y^e$ from Eq. (3).$E_x$is not considered here as we study the modulator device in 2 dimension (y-z). Thus, we obtain:

The interaction between the optical mode and the applied electric field can be characterized by the electro-optic overlap integral [18]:

As BTO thin film is the active EO layer in our waveguide geometry, the integral ${\displaystyle {\iint }_{BTO}{E}_y^e\left(E_z.E_y^{*}+E_y.E_z^{*}\right)dydz}$ is only calculated in the BTO layer. Note that because of the non-diagonal character of the permittivity tensor of the BTO layer, the electro-optic integral differs from that commonly used in literature ${\displaystyle {\iint }_{BTO}{E}_y^e{\left|E\right|}^{2}dydz}$ which is only valid for a diagonal matrix.

The effective Pockels coefficient r_eff can thus be extracted by [19]:

In the 2D case, the above equation is transformed into:

Numerical solutions of Eq. (13) are also obtained by finite difference method. Description of the method and values of the coefficients of the sparse linear system are given in Appendix 2. Solution of the sparse system of linear equations could be obtained using several iterative methods such as the Jacobi, Gauss–Seidel, as well as the successive over-relaxation method (SOR) [20,21]. As the value of the potential is known in the inner part of the metallic contacts, the corresponding nodes can be moved on the left part of the linear system. By this way, the problem is reduced to the inversion of the sparse matrix with a computation time of a few seconds compared to a few hours for the SOR method [22]. Transparent boundaries conditions are also used in order to reduce the size of the calculation window. Determination of the key RF parameters such as the microwave index n_m and the mapping of the electric field are derived from the calculation of the potential Φ [18,23].

4. Design of the BTO modulator

In this section, we use the theoretical framework described previously to model and simulate the optical mode and the RF electro-optical response of a BTO layer in a slot waveguide configuration. From these simulations, we design and optimize a hybrid slot-waveguide Mach-Zehnder Interferometer modulator that fully exploits the very high Pockels effect of BTO.

Geometrical parameters of the modulator are reported on Fig. 1. To simulate its electro-optical properties, we use realistic parameters based on state-of-the-art fabrication and technological processes. We start from a standard SOI wafer with a 2 µm buried SiO₂ and a 220-nm-thick top silicon layer. The active BTO layer is added on top of the SOI wafer. As explained before, we set its thickness to 20 nm in order to have a purely c-axis oriented monocrystalline BTO [9]. The refractive index of BTO (n_BTO ~2.3) being lower than that of silicon (n_Si~3.5), an additional amorphous silicon layer (a-Si) is deposited on top of the BTO layer to increase the optical power confinement in the active part of the modulator thanks to a slot waveguide effect. As the etching recipes of the BTO material are not well established in the literature, lateral confinement of the waveguide mode is realized by simply etching the a-Si layer. In this ridge-slot waveguide geometry, the optical beam propagates along the x-axis and the electric field is applied along the y-axis.

To design the modulator, we need to optimize several parameters such as the optical confinement in the active region, the propagation losses, the electro-optic overlap and the index-matched travelling wave electrodes. In the following, we split the design work in two subsections: optical design and microwave design.

4.1 Optical design

We first need to maximize the optical power in the active BTO layer. To do so, the degrees of freedom we have lie in the width (W) and the thickness of the a-Si ridge waveguide layer (T_a-Si). The refractive indices of the different materials used for the optical modelling are listed in Table 1 (a). The wavelength is fixed to 1.55 µm (standard telecommunication wavelength).

Table 1. Optical Constants of Thin Films of Materials

View Table

Figure 2 illustrates the variation of Fig. 2(a) the confinement factor Γ, Fig. 2(b) effective index $n_{eff}^o$ and Fig. 2(c) propagation losses versus waveguide geometry parameters W and T_a-Si for TM-like mode propagation.

 

Fig. 2 Optimization of the slot waveguide geometry: (a) Optical confinement factor Γ versus W and T_a-Si; (b) Optical effective index versus W and T_a-Si; (c) Optical loss versus W and T_a-Si; (d) Optimized TM-like mode power profile

Download Full Size | PDF

We can notice that the propagation losses are close to zero when the effective index becomes greater than 2.8845. Interestingly, this value corresponds to the effective index of the fundamental TE mode for the surrounding slab part. The propagation losses for this kind of waveguide are only due to leakage in the surrounding BaTiO₃ and Si layers. The boundary between low and high losses areas is represented by dashed lines in Fig. 2(b) and Fig. 2(c). As shown in Fig. 2(a), the upper part of the boundary that corresponds to the minimum losses presents a high confinement factor. A maximal confinement is obtained for W = 0.76 µm and T_a-Si = 0.25 µm, indicated by a yellow star, where Γ = 12.16%. The optical loss is extremely weak (7.41 × 10⁻⁶ dB/cm), and $n_{eff}^o$ = 2.9121.

Mapping of the optical power in Fig. 2(d) shows that the parts of the mode in the surrounding slab layers are well confined under the a-Si ridge layer. A mono-modal behavior is obtained by rejection of higher order modes due to leakage effect [30,31]. The same calculations were made for the TE-like mode, which gives $n_{eff}^o$ = 3.0936 and Γ = 5.27%.

Optimization of the spacing between the optical waveguide and the metallic contacts (D) is a tradeoff between maximizing electro-optic effect while maintaining low propagation losses. To mitigate the propagation losses caused by the proximity of contacts, and to light the polarization fatigue due to metal contacts [32,33], Indium Tin Oxide (ITO) is used rather than metal for the electrodes. Figure 3 shows the evolution of propagation losses versus D.

 

Fig. 3 Evolution of propagation losses of slot mode versus electrode spacing D

Download Full Size | PDF

The propagation loss of the slot mode gradually decreases as the gap D is increased. For a spacing D of about 0.3 µm, propagation losses become lower than 1.5 dB/cm. As the optical mode is strongly confined under the a-Si layer, variation of the electrodes width W_e and thickness T_e have negligible influence on the optical losses. Thanks to that, the technological constraints on these two parameters are relaxed but their geometries are optimized in the following, for RF considerations.

4.2 Microwave model optimization

Next, we consider the geometry of transmission lines. The thickness (T_e) and width (W_e) of the CoPlanar Wave electrodes (CPW) have a strong influence on the modulation bandwidth (BW) and the half-wave voltage V_π. The input refractive index of each material values at 10 GHz in the RF range [25–29] are listed in Table 1(b).

Evolution of the electro-optic overlap integral Γ_eo (as calculated in Eq. (10)) and microwave effective index $n_{eff}^m$ versus the thickness T_e and the electrode width W_e are plotted on Fig. 4(a) and Fig. 4(b). An optimal design for the electrodes is obtained for thickness and width of about 1.0 and 6.3 µm, respectively. This geometry enables phase matching between RF and optical signals and an optimal Γ_eo of 11.34%. When the matching condition is totally satisfied, the modulation speed is then only limited by the conductor loss. In order to have the maximum modulation bandwidth, velocity match condition should be satisfied. In other words, microwave effective index $n_{eff}^m$ should be near optical effective index $n_{eff}^o$ = 2.9121.

 

Fig. 4 Evolution of (a) Γ_eo and (b) microwave effective index $n_{eff}^m$ versus the thickness T_e and width W_e of the electrode for a gap D of about 0.3 µm. The dashed line indicates geometries where value of $n_{eff}^m$matches that of $n_{eff}^o$.

Download Full Size | PDF

Cartographies of electrical field components $E_y^e$ and$E_z^e$, are shown in Fig. 5. Bias voltage of 20V is applied between the right and left electrodes. Zoom of the dashed box part shows the variation of the electric field in the ridge optical waveguide where the optical mode is confined. We thus verify here that this coplanar geometry for electrodes results in a y-oriented electric field under the a-Si strip.

 

Fig. 5 Electric voltage plot with 20 V applied across the electrodes and the gap spacing between the two electrodes G is 1.36 µm with D = 0.3µm.

Download Full Size | PDF

4.3 Electro-optical response with optimized modulator design

Coupling Eq. (2) with the calculated electric field, we can simulate the spatial variation of the permittivity tensor modifications induced by the Pockels effect throughout the active BTO. Variation of $n_{eff}^o$ is then calculated without approximation and directly give us the interaction length L from

To understand the physics of the system, we want to investigate how the applied electric field affects the individual elements of the permittivity tensor of BTO. The spatial distribution of modifications to the five tensor components Δε_xx, Δε_yy, Δε_zz, Δε_yz, Δε_zy, calculated from Eq. (2) with bias voltage of 20 V are shown in Fig. 6. Variations of the non-diagonal elements Δε_yz and Δε_zy (Fig. 6(a)) are strikingly larger than the diagonal ones. This is well explained by both the symmetry of the crystalline BTO layer and the direction of the electric field respective to this crystal orientation. Indeed, both the largest Pockels coefficient, r₄₂ (820 pm/V) and the strongest component of the electric field ($E_y^e$) are located in non-diagonal tensor elements. Retrospectively, we therefore verify that the coplanar electrodes geometry is well adapted to the c-axis oriented BTO.

 

Fig. 6 Variation of the five non-zero permittivity tensor elements Δε_xx, Δε_yy, Δε_zz, Δε_yz, Δε_zy within BTO layer as a result of a 20 V bias voltage and D = 0.3µm.

Download Full Size | PDF

We now want to evidence the Pockels effect in our design by sweeping the bias voltage from 0 V to 20 V. TM-like mode and TE-like mode are respectively studied and compared. As displayed in Fig. 7(a), $\Delta n_{eff}^o$shows a quadratic increase with the bias voltage, due to the linear variation of the impermeability 1/ε with RF signal. In the meanwhile, the TM-like slot mode shows a more efficient electro-optical modulation than TE-like mode, with a larger $\Delta n_{eff}^o$for a same bias voltage. The evolution of a key parameter for modulators, V_πL, is shown as well in Fig. 7(b). Higher bias voltage decreases the interaction length L.

 

Fig. 7 (a) The optical effective index variation $\Delta n_{eff}^o$increases gradually with bias voltage from 0 to 20V, and TM-like mode shows a higher variation compared to TE-like mode; (b) Key optical modulator parameter Vπ·L decreases gradually with bias voltage from 0 to 20V, and TM-like mode shows a smaller Vπ·L value.

Download Full Size | PDF

With a bias voltage of 10 V for TM-mode, Δ$n_{eff}^o$ = 1.00 10⁻³ and V_πL = 0.78 V·cm. By applying a driving voltage $\pm$1V with 10 V DC bias, we obtain Δ$n_{eff}^o$ = 2.75 10⁻⁴.

As expected with the optical optimization, Δ$n_{eff}^o$is twice higher for the TM-like mode compared to the TE-like mode.

We then calculate the effective Pockels coefficient r_eff by introducing V_πL value shown in Fig. 7(b) and others optimized parameters in Eq. (11). The maximal value extracted is as high as r_eff = 276.59 pmV⁻¹ for a bias voltage of 20 V, which is comparable to recent experimental results [6,7]. Note that, this value was obtained by calculating the permittivity tensor modification without approximation.

5. Conclusion

We have developed a multi-physics model that hybridizes a full-vectorial mode solver with a Laplace RF mode solver. This numerical technique enables precise calculations of optical confinement, index ellipsoid modulation and electro-optic interaction in anisotropic materials. This is particularly important if one wants to leverage the promising electro-optic properties of anisotropic materials such as perovskite ferroelectric oxides. We highlighted the possibilities of our method by designing a structure that can exploit the large Pockels effect of ferroelectric BTO in an electro-optic modulator. We first optimized the optical mode confinement in the active BTO layer using a slot-waveguide geometry which supports a TM-like mode with very low theoretical propagation losses. Next, we optimized the electrodes geometry to simultaneously maximize the electro-optic response and the modulation bandwidth. To that aim, we rigorously calculated the modifications to the permittivity tensor of a c-axis oriented BTO thin film. This calculation reveals that the largest Pockels coefficient (r₄₂) is mainly located in off-diagonal tensor elements and excited by the y-component of an arbitrary electric field. We therefore designed coplanar electrodes that maximize the in-plane electric field in the BTO film while maintaining low propagation losses. In addition, we isolated the geometries needed for these CPW electrodes to have a microwave effective index matching the optical effective index of the slot mode. The proposed hybrid silicon-ferroelectric slot-waveguide modulator has a theoretical Δn_eff of ~2.8.10⁻³ and a r_eff of 276 pm.V⁻¹ at 20 V.

The described numerical method could be used to fully understand and exploit the promising electro-optic properties of a variety of anisotropic materials in useful photonic devices.

Appendix 1

The finite difference differential operators a_ij and b_ij could be expressed in term of the transverse electric and magnetic field at the center point P under consideration and at the eight neighboring grid points (see Fig. 8):

(15)$$a_{ij}={\displaystyle \sum _{X=P,N,S,W,E,NW,NE,SW,SE}{a}_{{}_{ij}}^X{H}_{{}_j}^X},$$

(16)$$b_{ij}={\displaystyle \sum _{X=P,N,S,W,E,NW,NE,SW,SE}{b}_{{}_{ij}}^X{E}_{{}_j}^X}$$

where i and j denote the y and z axes. The expressions of $a_{{}_{ij}}^X$and $b_{ij}^X$are as follows:

 

Fig. 8 Nomenclature convention of the grid point and meshing for the finite difference scheme

Download Full Size | PDF

$$a_{yy}^P=k_0^2{\epsilon }_{zz}^P-\left(\frac{2}{ew}\right)-\left(\frac{4s{\epsilon }_{zz}^P}{n^2\left(n+s\right)\left({\epsilon }_{xx}^N+{\epsilon }_{xx}^P\right)}+\frac{4n{\epsilon }_{zz}^P}{s^2\left(n+s\right)\left({\epsilon }_{xx}^S+{\epsilon }_{xx}^P\right)}-\frac{2{\epsilon }_{zz}^P{\left(n-s\right)}^2}{n^2{s}^2{\epsilon }_{xx}^P}\right)+\frac{\left(e-w\right).\left(n-s\right)}{e.w.n.s}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^P}\right).$$

$$a_{yy}^W=\left(\frac{2}{w\left(e+w\right)}\right)+\frac{e.\left(n-s\right)}{w.n.s.\left(e+w\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^W}\right).a_{yy}^E=\left(\frac{2}{e\left(e+w\right)}\right)-\frac{w.\left(n-s\right)}{e.n.s.\left(e+w\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^E}\right).$$

$$a_{yy}^S=\left(\frac{4n{\epsilon }_{zz}^P}{s^2\left(n+s\right)\left({\epsilon }_{xx}^S+{\epsilon }_{xx}^P\right)}-\frac{2\left(n-s\right){\epsilon }_{zz}^P}{s^2\left(n+s\right){\epsilon }_{xx}^P}\right)+\frac{n.\left(e-w\right)}{e.w.s.\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^P}\right).$$

$$a_{yy}^N=\left(\frac{4s{\epsilon }_{zz}^P}{n^2\left(n+s\right)\left({\epsilon }_{xx}^N+{\epsilon }_{xx}^P\right)}+\frac{2\left(n-s\right){\epsilon }_{zz}^P}{n^2\left(n+s\right){\epsilon }_{xx}^P}\right)-\frac{s.\left(e-w\right)}{e.w.n.\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^P}\right).$$

$$a_{yy}^{NE}=\frac{w.s}{e.n.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^E}\right).a_{yy}^{SE}=-\frac{w.n}{e.s.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^E}\right).$$

$$a_{yy}^{NW}=-\frac{e.s}{w.n.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^W}\right).a_{yy}^{SW}=\frac{e.n}{w.s.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^W}\right).$$

$$a_{yz}^P=-k_0^2{\epsilon }_{zy}^P+\left(\frac{4w{\epsilon }_{zy}^P}{e^2\left(e+w\right)\left({\epsilon }_{xx}^E+{\epsilon }_{xx}^P\right)}+\frac{4e{\epsilon }_{zy}^P}{w^2\left(e+w\right)\left({\epsilon }_{xx}^W+{\epsilon }_{xx}^P\right)}-\frac{2{\epsilon }_{zy}^P{\left(e-w\right)}^2}{e^2{w}^2{\epsilon }_{xx}^P}\right)+\frac{\left(e-w\right).\left(n-s\right)}{e.w.n.s}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^P}\right).$$

$$a_{yz}^W=-\left(\frac{4e{\epsilon }_{zy}^P}{w^2\left(e+w\right)\left({\epsilon }_{xx}^W+{\epsilon }_{xx}^P\right)}-\frac{2\left(e-w\right){\epsilon }_{zy}^P}{w^2\left(e+w\right){\epsilon }_{xx}^P}\right)+\frac{e.\left(n-s\right)}{w.n.s.\left(e+w\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^P}\right).$$

$$a_{yz}^E=-\left(\frac{4w{\epsilon }_{zy}^P}{e^2\left(e+w\right)\left({\epsilon }_{xx}^E+{\epsilon }_{xx}^P\right)}+\frac{2\left(e-w\right){\epsilon }_{zy}^P}{e^2\left(e+w\right){\epsilon }_{xx}^P}\right)-\frac{w.\left(n-s\right)}{e.n.s.\left(e+w\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^P}\right).$$

$$a_{yz}^S=\frac{n.\left(e-w\right)}{e.w.s.\left(n+s\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^S}\right).a_{yz}^N=-\frac{s.\left(e-w\right)}{e.w.n.\left(n+s\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^N}\right).$$

$$a_{yz}^{NE}=\frac{w.s}{e.n.\left(e+w\right)\left(n+s\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^N}\right).a_{yz}^{SE}=-\frac{w.n}{e.s.\left(e+w\right)\left(n+s\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^S}\right).$$

$$a_{yz}^{NW}=-\frac{e.s}{w.n.\left(e+w\right)\left(n+s\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^N}\right).a_{yz}^{SW}=\frac{e.n}{w.s.\left(e+w\right)\left(n+s\right)}\left(1-\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^S}\right).$$

$$b_{yy}^P=k_0^2{\epsilon }_{yy}^P-\left(\frac{4w{\epsilon }_{yy}^P}{e^2\left(e+w\right)\left({\epsilon }_{xx}^E+{\epsilon }_{xx}^P\right)}+\frac{4e{\epsilon }_{yy}^P}{w^2\left(e+w\right)\left({\epsilon }_{xx}^W+{\epsilon }_{xx}^P\right)}-\frac{2{\epsilon }_{yy}^P{\left(e-w\right)}^2}{e^2{w}^2{\epsilon }_{xx}^P}\right)-\left(\frac{2}{ns}\right)+\frac{\left(e-w\right).\left(n-s\right)}{e.w.n.s}\left(\frac{{\epsilon }_{zy}^P}{{\epsilon }_{xx}^P}\right).$$

$$b_{yy}^W=\left(\frac{4e{\epsilon }_{yy}^W}{w^2\left(e+w\right)\left({\epsilon }_{xx}^W+{\epsilon }_{xx}^P\right)}-\frac{2\left(e-w\right){\epsilon }_{yy}^W}{w^2\left(e+w\right){\epsilon }_{xx}^P}\right)+\frac{e.\left(n-s\right)}{w.n.s.\left(e+w\right)}\left(\frac{{\epsilon }_{zy}^W}{{\epsilon }_{xx}^W}\right).$$

$$b_{yy}^E=\left(\frac{4w{\epsilon }_{yy}^E}{e^2\left(e+w\right)\left({\epsilon }_{xx}^E+{\epsilon }_{xx}^P\right)}+\frac{2\left(e-w\right){\epsilon }_{yy}^E}{e^2\left(e+w\right){\epsilon }_{xx}^P}\right)-\frac{w.\left(n-s\right)}{e.n.s.\left(e+w\right)}\left(\frac{{\epsilon }_{zy}^E}{{\epsilon }_{xx}^E}\right).$$

$$b_{yy}^S=\left(\frac{2}{s\left(n+s\right)}\right)+\frac{n.\left(e-w\right)}{e.w.s.\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^S}{{\epsilon }_{xx}^P}\right).b_{yy}^N=\left(\frac{2}{n\left(n+s\right)}\right)-\frac{s.\left(e-w\right)}{e.w.n.\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^N}{{\epsilon }_{xx}^P}\right).$$

$$b_{yy}^{NE}=\frac{1}{\left(e+w\right)\left(n+s\right)}\left(\frac{w.s.{\epsilon }_{zy}^{NE}}{e.n.{\epsilon }_{xx}^E}\right).b_{yy}^{SE}=-\frac{w.n}{e.s.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^{SE}}{{\epsilon }_{xx}^E}\right).$$

$$b_{yy}^{NW}=-\frac{e.s}{w.n.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^{NW}}{{\epsilon }_{xx}^W}\right).b_{yy}^{SW}=\frac{e.n}{w.s.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zy}^{SW}}{{\epsilon }_{xx}^W}\right).$$

$$b_{yz}^P=k_0^2{\epsilon }_{yz}^P-\left(\frac{4w{\epsilon }_{yz}^P}{e^2\left(e+w\right)\left({\epsilon }_{xx}^E+{\epsilon }_{xx}^P\right)}+\frac{4e{\epsilon }_{yz}^P}{w^2\left(e+w\right)\left({\epsilon }_{xx}^W+{\epsilon }_{xx}^P\right)}-\frac{2{\epsilon }_{yz}^P{\left(e-w\right)}^2}{e^2{w}^2{\epsilon }_{xx}^P}\right)+\frac{\left(e-w\right).\left(n-s\right)}{e.w.n.s}\left(\frac{{\epsilon }_{zz}^P}{{\epsilon }_{xx}^P}-1\right).$$

$$b_{yz}^W=\left(\frac{4e{\epsilon }_{yz}^W}{w^2\left(e+w\right)\left({\epsilon }_{xx}^W+{\epsilon }_{xx}^P\right)}-\frac{2\left(e-w\right){\epsilon }_{yz}^W}{w^2\left(e+w\right){\epsilon }_{xx}^P}\right)+\frac{e.\left(n-s\right)}{w.n.s.\left(e+w\right)}\left(\frac{{\epsilon }_{zz}^W}{{\epsilon }_{xx}^W}-1\right).$$

$$b_{yz}^E=\left(\frac{4w{\epsilon }_{yz}^E}{e^2\left(e+w\right)\left({\epsilon }_{xx}^E+{\epsilon }_{xx}^P\right)}+\frac{2\left(e-w\right){\epsilon }_{yz}^E}{e^2\left(e+w\right){\epsilon }_{xx}^P}\right)-\frac{w.\left(n-s\right)}{e.n.s.\left(e+w\right)}\left(\frac{{\epsilon }_{zz}^E}{{\epsilon }_{xx}^E}-1\right).$$

$$b_{yz}^S=\frac{n.\left(e-w\right)}{e.w.s.\left(n+s\right)}\left(\frac{{\epsilon }_{zz}^S}{{\epsilon }_{xx}^P}-1\right).b_{yz}^N=-\frac{s.\left(e-w\right)}{e.w.n.\left(n+s\right)}\left(\frac{{\epsilon }_{zz}^N}{{\epsilon }_{xx}^P}-1\right).$$

$$b_{yz}^{NE}=\frac{w.s}{e.n.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zz}^{NE}}{{\epsilon }_{xx}^E}-1\right).b_{yz}^{SE}=-\frac{w.n}{e.s\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zz}^{SE}}{{\epsilon }_{xx}^E}-1\right).$$

$$b_{yz}^{NW}=-\frac{e.s}{w.n.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zz}^{NW}}{{\epsilon }_{xx}^W}-1\right).b_{yz}^{SW}=\frac{e.n}{w.s.\left(e+w\right)\left(n+s\right)}\left(\frac{{\epsilon }_{zz}^{SW}}{{\epsilon }_{xx}^W}-1\right).$$

We can notice that the other coefficients ${\text{a}}_{\text{zy}}^{\text{x}}$, ${\text{a}}_{\text{zz}}^{\text{x}}$, ${\text{b}}_{\text{zy}}^{\text{x}}$ and ${\text{b}}_{\text{zz}}^{\text{x}}$ are obtained respectively with the formulas of coefficients ${\text{a}}_{\text{yz}}^{\text{x}}$ and ${\text{a}}_{\text{yy}}^{\text{x}}$, ${\text{b}}_{\text{yz}}^{\text{x}}$ and ${\text{b}}_{\text{yy}}^{\text{x}}$ with a permutation of subscripts $\text{y}\to \text{z}$, $\text{z}\to \text{y}$, $\text{E}\to \text{N}$, $\text{W}\to \text{S}$, $\text{N}\to \text{E}$, $\text{S}\to \text{W}$, $\text{e}\to \text{n}$, $\text{w}\to \text{s}$, $\text{n}\to \text{e}$, $\text{s}\to \text{w}$.

Appendix 2

Eq. (12) could be expressed as:

(17)$${\displaystyle \sum _{X=P,N,S,W,E,NW,NE,SW,SE}{a}_x{\Phi }_X}=0$$

$$\begin{array}{l}{a}_P = -\left(\frac{w\left({\epsilon }_{yy}^E+{\epsilon }_{yy}^P\right)}{e^2\left(e+w\right)}+\frac{e\left({\epsilon }_{yy}^W+{\epsilon }_{yy}^P\right)}{w^2\left(e+w\right)}-\frac{{\left(e-w\right)}^2{\epsilon }_{yy}^P}{e^2{w}^2}\right)-\left(\frac{s\left({\epsilon }_{zz}^N+{\epsilon }_{zz}^P\right)}{n^2\left(n+s\right)}+\frac{n\left({\epsilon }_{zz}^S+{\epsilon }_{zz}^P\right)}{s^2\left(n+s\right)}-\frac{{\left(n-s\right)}^2{\epsilon }_{zz}^P}{n^2{s}^2}\right)\\ +\frac{\left(e^2-w^2\right)\left(n^2-s^2\right)\left({\epsilon }_{yz}^P+{\epsilon }_{zy}^P\right)}{\left(e+w\right)\left(n+s\right)ewns}.\end{array}$$

$$a_E = \left(\frac{w{\epsilon }_{yy}^E+e{\epsilon }_{yy}^P}{e^2\left(e+w\right)}\right)+\frac{w\left(n^2-s^2\right)\left({\epsilon }_{yz}^E+{\epsilon }_{zy}^P\right)}{ens}.a_W = \frac{e{\epsilon }_{yy}^W+w{\epsilon }_{yy}^P}{w^2\left(e+w\right)}-\frac{e\left(n^2-s^2\right)\left({\epsilon }_{yz}^W+{\epsilon }_{zy}^P\right)}{wns}.$$

$$a_N = \left(\frac{s{\epsilon }_{zz}^N+n{\epsilon }_{zz}^P}{n^2\left(n+s\right)}\right)+\frac{s\left(e^2-w^2\right)\left({\epsilon }_{yz}^P+{\epsilon }_{zy}^N\right)}{ewn}.a_s = \left(\frac{n{\epsilon }_{zz}^S+s{\epsilon }_{zz}^P}{s^2\left(n+s\right)}\right)-\frac{n\left(e^2-w^2\right)\left({\epsilon }_{yz}^P+{\epsilon }_{zy}^S\right)}{ews}.$$

$$a_{NE} = \frac{ws}{en}\frac{{\epsilon }_{yz}^E+{\epsilon }_{zy}^N}{\left(e+w\right)\left(n+s\right)}.a_{SE} = -\frac{wn}{es}\frac{{\epsilon }_{yz}^E+{\epsilon }_{zy}^S}{\left(e+w\right)\left(n+s\right)}.a_{NW} = -\frac{es}{wn}\frac{{\epsilon }_{yz}^W+{\epsilon }_{zy}^N}{\left(e+w\right)\left(n+s\right)}.a_{SW} = \frac{en}{ws}\frac{{\epsilon }_{yz}^W+{\epsilon }_{zy}^N}{\left(e+w\right)\left(n+s\right)}.$$

Determination of the B vector is illustrated on Fig. 9

 

Fig. 9 Mesh grid of a metallic contact (i₁ = 7, i₂ = 12 and j₁ = 5, j₂ = 8).

Download Full Size | PDF

The boundaries of the metallic contact are in blue on Fig. 9. Determination of the B vector in Eq. (17) is made with the points of the blue area on Fig. 9.

If the potential of the metallic contact is equal to V_c,

• - at point (j₁-1, i₁-1), $B\left(j,i\right)= -a_{NE}{V}_c$,
• - at point (j₁-1, i₁), $B\left(j,i\right)= -\left(a_{NE}+a_N\right)V_c$,
• - at points between (j₁-1, i₁ + 1) and (j₁-1,i₂-1), $B\left(j,i\right)= -\left(a_{NE}+a_N+a_{NW}\right)V_c$,
• - at point (j₁, i₁), $B\left(j,i\right)= -\left(a_{NE}+a_E+a_N+a_P\right)V_c$,
• - at points between (j₁, i₁ + 1) and (j₁,i₂-1), $B\left(j,i\right)= -\left(a_{NE}+a_N+a_{NW}+a_P+a_W+a_E\right)V_c$.

The part of the linear system where all the coefficients are constant corresponding to the red area of Fig. 9 is equal to $V_i^P = V_c$

Acknowledgment

This work was supported by the European STREP program “FP7-ICT-2013-11-619456-SITOGA”

References and links

1. B. Chmielak, M. Waldow, C. Matheisen, C. Ripperda, J. Bolten, T. Wahlbrink, M. Nagel, F. Merget, and H. Kurz, “Pockels effect based fully integrated, strained silicon electro-optic modulator,” Opt. Express 19(18), 17212–17219 (2011). [CrossRef]   [PubMed]  

2. D. J. Thomson, F. Y. Gardes, J.-M. Fedeli, S. Zlatanovic, Y. Hu, B. P. P. Kuo, E. Myslivets, N. Alic, S. Radic, G. Z. Mashanovich, and G. T. Reed, “50-Gb/s silicon optical modulator,” IEEE Photon. Technol. Lett. 24(4), 234–236 (2012). [CrossRef]  

3. M. Zgonik, P. Bernasconi, M. Duelli, R. Schlesser, P. Günter, M. Garrett, D. Rytz, Y. Zhu, and X. Wu, “Dielectric, elastic, piezoelectric, electro-optic, and elasto-optic tensors of BaTiO3 crystals,” Phys. Rev. B 50(9), 5941–5949 (1994).

4. B. W. Wessels, “Ferroelectric epitaxial thin films for integrated optics,” Annu. Rev. Mater. Res. 37(1), 659–679 (2007). [CrossRef]  

5. G. Niu, S. Yin, G. Saint-Girons, B. Gautier, P. Lecoeur, V. Pillard, G. Hollinger, and B. Vilquin, “Epitaxy of BaTiO3 thin film on Si (001) using a SrTiO3 buffer layer for non-volatile memory application,” Microelectron. Eng. 88(7), 1232–1235 (2011). [CrossRef]  

6. S. Abel, T. Stöferle, C. Marchiori, C. Rossel, M. D. Rossell, R. Erni, D. Caimi, M. Sousa, A. Chelnokov, B. J. Offrein, and J. Fompeyrine, “A strong electro-optically active lead-free ferroelectric integrated on silicon,” Nat. Commun. 4, 1671 (2013). [CrossRef]   [PubMed]  

7. C. Xiong, W. H. Pernice, J. H. Ngai, J. W. Reiner, D. Kumah, F. J. Walker, C. H. Ahn, and H. X. Tang, “Active silicon integrated nanophotonics: ferroelectric BaTiO3 devices,” Nano Lett. 14(3), 1419–1425 (2014). [CrossRef]   [PubMed]  

8. R. Syms and J. Cozens Syms, Optical Guided Waves and Devices (McGraw-Hill, 1992).

9. S. Abel, M. Sousa, C. Rossel, D. Caimi, M. D. Rossell, R. Erni, J. Fompeyrine, and C. Marchiori, “Controlling tetragonality and crystalline orientation in BaTiO3 nano-layers grown on Si,” Nanotechnology 24(28), 285701 (2013). [CrossRef]   [PubMed]  

10. J. M. Liu, Photonic Devices (Cambridge University, 2005.

11. J. D. Bull and N. A. F. Jaeger, “Parasitic Mode Conversion in Z-Propagating Lithium-Niobate Waveguides,” J. Lightwave Technol. 25(1), 387–393 (2007). [CrossRef]  

12. B. Fallahkhair, K. S. Li, and T. E. Murphy, “Vector Finite Difference Mode solver for Anisotropic Dielectric Waveguides,” J. Lightwave Technol. 26(11), 1423–1431 (2008). [CrossRef]  

13. P. Bienstman, S. Selleri, I. Rosa, H. P. Uranus, W. C. Hopman, R. Costa, A. Melloni, C. Andreani, J. P. Hugonin, P. Lalanne, D. Pinto, S. S. A. Obayya, M. Dems, and K. Panajotov, “Modelling leaky photonic wires: A mode solver comparison,” Opt. Quantum Electron. 38(9-11), 731–759 (2007). [CrossRef]  

14. M. Masi, R. Orobtchouk, G. F. Fan, and L. Pavesi, “Towards a Realistic Modelling of Ultra-Compact Racetrack Resonators,” J. Lightwave Technol. 28, 3233–3242 (2010).

15. R. Orobtchouk, P. Labeye, X. Hu, S. Malhouitre, Ph. Grosse, and J.-M. Fedeli, “Design, realization, and characterization of a silicon photonics coherent mixer for PDM-QPSK optical communications,” presented at the Photonics Europeen Optical Society Annual Meeting TOM 2-Silicon, Scotland, 25–28 September 2012.

16. G. R. Hadley, “Transparent boundary condition for the beam propagation method,” IEEE J. Quantum Electron. 28(1), 363–370 (1992). [CrossRef]  

17. H. Kogelnik, “Theory of optical waveguides,” in Guided Wave Optoelectronics, 2nd Ed., T. Tamir ed. (Springer, 1990,), pp. 7–88.

18. M. Minakata, “Recent progress of 40-GHz high-speed LiNbO3 optical modulator,” ITCom 2001: Intern. Symp. on the Convergence of IT and Communications. International Society for Optics and Photonics (2001).

19. P. Tang, D. J. Towner, A. L. Meier, and B. W. Wessels, “Low-voltage, polarization-insensitive, electro-optic modulator based on a polydomain barium titanate thin film,” Appl. Phys. Lett. 85(20), 4615–4617 (2004). [CrossRef]  

20. J. H. Mathews, Numerical Methods for Mathematics, Science and Engineering (Prentice-Hall, 1992).

21. T. E. Shoup, A Practical Guide to Computer Methods for Engineers (Prentice-Hall, 1979).

22. S. Rosloniec, Fundamental Numerical Methods for Electrical Engineering (Springer, 2008).

23. G. Fan, Y. Li, B. Han, Q. Wang, X. Liu, and Z. Zhen, “Study of an Electro-Optic Polymer Modulator,” J. Lightwave Technol. 30(15), 2482–2487 (2012). [CrossRef]  

24. E. D. Palik, Handbook of Optical Constants of Solids, (Academic press, 1998, Vol. III).

25. M. N. Afsar and K. J. Button, “Precise Millimeter-Wave Measurements of Complex Refractive Index, Complex Dielectric Permittivity and Loss Tangent of GaAs, Si, SiO₂, A1₂O₃, BeO, Macor, and Glass,” in Proceedings of IEEE Conference on Microwave Theory and Techniques 31, 217–223 (1983).

26. C. Rhodes, M. Cerruti, A. Efremenko, M. Losego, D. E. Aspnes, J. P. Maria, and S. Franzen, “Dependence of plasmon polaritons on the thickness of indium tin oxide thin films,” J. Appl. Phys. 103(9), 093108 (2008). [CrossRef]  

27. M. D. Janezic, D. F. Williams, V. Blaschke, A. Karamcheti, and C. S. Chang, “Permittivity characterization of low-k thin films from transmission-line measurements,” in Proceedings of IEEE Conference on Transactions on Microwave Theory and Techniques 51, 132–136 (2003). [CrossRef]  

28. D. G. Sun, Z. Liu, Y. Huang, S. T. Ho, D. J. Towner, and B. W. Wessels, “Performance simulation for ferroelectric thin-film based waveguide electro-optic modulators,” Opt. Commun. 255(4-6), 319–330 (2005). [CrossRef]  

29. T. Hamano, D. J. Towner, and B. W. Wessels, “Relative dielectric constant of epitaxial BaTiO₃ thin films in the GHz frequency range,” Appl. Phys. Lett. 83(25), 5274–5276 (2003). [CrossRef]  

30. S. P. Pogossian, L. Vescan, and A. Vonsovici, “The single-mode condition for semiconductor rib waveguides with large cross section,” J. Lightwave Technol. 16(10), 1851–1853 (1998). [CrossRef]  

31. O. Powell, “Single-Mode Condition for Silicon Rib Waveguides,” J. Lightwave Technol. 20(10), 1851–1855 (2002). [CrossRef]  

32. H. M. Duiker, P. D. Beale, J. F. Scott, C. A. Paz de Araujo, B. M. Melnick, J. D. Cuchiaro, and L. D. McMillan, “Fatigue and switching in ferroelectric memories - theory and experiment,” J. Appl. Phys. 68(11), 5783–5791 (1990). [CrossRef]  

33. J. F. Scott, Ferroelectric Memories (Springer, 2000).