Fully-Tensorial Elastic-Wave Mode-Solver in FEniCS for Stimulated Brillouin Scattering Modeling

TL;DR

This paper introduces a comprehensive elastic-wave mode solver in FEniCS that accounts for the full tensorial stiffness, enabling detailed simulations of stimulated Brillouin scattering in anisotropic waveguides.

Contribution

It presents a novel tensorial elastic-wave mode solver in FEniCS, including various approximations and coupling with electromagnetic models for Brillouin scattering analysis.

Findings

Validated the full tensorial model against approximations

Demonstrated the impact of elastic anisotropy on Brillouin gain

Provided an open-source tool for integrated elastic and optical simulations

Abstract

A framework for simulating the elastic-wave modes in waveguides, taking into account the full tensorial nature of the stiffness tensor, is presented and implemented in the open-source finite element solver, FEniCS. Various approximations of the elastic wave equation used in the stimulated Brillouin scattering literature are implemented and their validity and applicability are discussed. The elastic mode-solver is also coupled with an electromagnetic counterpart to study the influence of elastic anisotropies on Brillouin gain.

Tables

Table 1. Table I : Si Elastic Modes

Mode	Fund. x shear	Fund. x shear	Torsional
SBS	Forward	Forward	Forward
Propagation axis	[100]	[110]	[100]
Eigenvalue, this work	8.4064 GHz	9.2300 GHz	10.633 GHz
Eigenvalue, 3D COMSOL	8.4064 GHz	9.2300 GHz	10.633 GHz
Plot	Fig. 4(a)	-	Fig. 4(b)

Equations

$$\nabla \cdot \mathbf{\sigma }=\rho \frac{{\partial }^2\mathbf{u}}{{\partial }^{2}t},$$

$$\mathbf{\sigma }=\mathbf{C}:\mathbf{ϵ}\text{and}\mathbf{ϵ}=\frac{1}{2}[\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T].$$

$$\begin{array}{rl}& {\nabla }_T\cdot \mathbf{C}:\frac{1}{2}[{\nabla }_T\mathbf{u}+({\nabla }_T\mathbf{u}{)}^T]+i{q}_b\hat{\mathbf{z}}\cdot \mathbf{C}:\frac{1}{2}[{\nabla }_T\mathbf{u}+({\nabla }_T\mathbf{u}{)}^T]+\\ & {\nabla }_T\cdot \mathbf{C}:\frac{1}{2}[i{q}_b\mathbf{u}+(i{q}_b\mathbf{u}{)}^T]+i{q}_b\hat{\mathbf{z}}\cdot \mathbf{C}:\frac{1}{2}[i{q}_b\mathbf{u}+(i{q}_b\mathbf{u}{)}^T]=\\ & -\rho {\omega }^2\mathbf{u}.\end{array}$$

$$\begin{array}{rl}& -{\int }_{\Omega }{\mathbf{ϵ}}_T({\mathbf{v}}^{\ast }):\mathbf{C}:{\mathbf{ϵ}}_T(\mathbf{u})dx+{\int }_{\Omega }{\mathbf{ϵ}}_Z({\mathbf{v}}^{\ast }):\mathbf{C}:{\mathbf{ϵ}}_T(\mathbf{u})dx\\ +& {\int }_{\Omega }{\mathbf{ϵ}}_T({\mathbf{v}}^{\ast }):\mathbf{C}:{\mathbf{ϵ}}_Z(\mathbf{u})dx+{\int }_{\Omega }{\mathbf{ϵ}}_Z({\mathbf{v}}^{\ast }):\mathbf{C}:{\mathbf{ϵ}}_Z(\mathbf{u})dx\\ & +{\int }_{\partial \Omega }{\mathbf{v}}^{\ast }\cdot [\hat{\mathbf{n}}\cdot \mathbf{C}:{\mathbf{ϵ}}_T(\mathbf{u})]dx=-{\int }_{\Omega }\rho {\omega }^2{\mathbf{v}}^{\ast }\cdot \mathbf{u}dx\end{array},$$

$$\left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{yz}\\ {\sigma }_{xz}\\ {\sigma }_{xy}\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{21}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{31}& c_{32}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{41}& c_{42}& c_{43}& c_{44}& c_{45}& c_{46}\\ c_{51}& c_{52}& c_{53}& c_{54}& c_{55}& c_{56}\\ c_{61}& c_{62}& c_{63}& c_{64}& c_{65}& c_{66}\end{array}\right]\left[\begin{array}{c}{ϵ}_{xx}\\ {ϵ}_{yy}\\ {ϵ}_{zz}\\ 2{ϵ}_{yz}\\ 2{ϵ}_{xz}\\ 2{ϵ}_{xy}\end{array}\right].$$

$$v_l^2{\nabla }_T^2\tilde{\rho }-q_b^2{v}_l^2\tilde{\rho }=-{\omega }^2\tilde{\rho },$$

$$\mu {\nabla }_T^2{u}_z-q_b^2(2\mu +\lambda )u_z=-\rho {\omega }^2{u}_z,$$

$$\begin{array}{r}F(E)={\int }_{\Omega }\frac{1}{{\mu }_r}({\nabla }_T\times E_T)\cdot ({\nabla }_T\times E_T)-k_0^2{ϵ}_r{E}_T\cdot E_T+\\ {\gamma }^2[\frac{1}{{\mu }_r}({\nabla }_T{E}_{z,\gamma }+E_T)\cdot ({\nabla }_T{E}_{z,\gamma }+E_T)-k_0^2{ϵ}_r{E}_z\cdot E_z]dx,\end{array}$$

$${\sigma }_{ij}^{ES}=-\frac{1}{4}{ϵ}_0{n}^4{p}_{ijkl}(E_{pk}{E}_{sl}^{\ast }+E_{pl}{E}_{sk}^{\ast }),$$

$${\mathbf{f}}^{\mathbf{ES}}=-{\nabla }_T\cdot {\mathbf{\sigma }}^{\mathbf{ES}}-i{q}_b\hat{\mathbf{z}}\cdot {\mathbf{\sigma }}^{\mathbf{ES}}.$$

$$f_i^{BES}=({\sigma }_{1ij}^{ES}-{\sigma }_{2ij}^{ES})n_j,$$

$$T_{ij}=\frac{1}{2}{ϵ}_0{ϵ}_r((E_{pi}{E}_{sj}^{\ast }+E_{pj}{E}_{si}^{\ast })-{\delta }_{ij}(E_{pk}{E}_{sk}^{\ast })),$$

$$f_i^{RP}=(T_{2ij}-T_{2ij})n_j.$$

$$G=\frac{{\omega }_{opt}{Q}_{mech}}{4P_{mech}{P}_s{P}_p}|{\int }_{\Omega }\mathbf{f}\cdot {\mathbf{u}}^{\ast }dx{|}^2,$$

Full text

I Introduction

Stimulated Brillouin scattering (SBS) has been extensively studied in optical fibers, albeit it is generally considered a nuisance for the long-haul optical communication application [1]. Recently, there has been renewed interest in SBS in order to harness it in integrated photonic devices [2, 3, 4, 5]. The narrow bandwidth of the SBS response leads to a host of applications in microwave photonics [6], such as the construction of tunable bandpass [7] and notch [8] filters, phase shifters [9] and microwave synthesizers [10]. Also, long phonon lifetime has been used to store optical pulses in the acoustic domain [11]. SBS has been, furthermore, utilized in wavelength-selective amplifiers [12] and narrow-linewidth lasers [13]. Cascaded SBS process has been proposed as a multiwavelength source for optical communications [14].

Unlike silica, the common material in optical fibers, most of the materials used in integrated photonics are not isotropic. A prime example is silicon with a cubic lattice structure [15], on which large Brillouin gain has been demonstrated [3]. To fully explore the whole parameter space of integrated acoustooptic devices, knowledge of all acoustic modes is necessary, especially given that the acoustic and optical modes can be tailored independently, while retaining coupling between the two classes of modes as in the case of a Brillouin laser on silicon membranes [16].

The finite-element method (FEM) is a versatile modeling choice for implementing an elastic-wave mode-solver that takes into account the tensorial stiffness of materials. The method has been originally developed in the 1960s and 70s to model the problems of structural mechanics [17], including the elastic-wave equation [18], sometimes also called the seismic wave equation. Unlike the finite-difference method, FEM can easily handle complex geometries. Among other modeling problems, FEM has also been used to find the acoustic modes of structures in the ultrasound and the GHz ranges. Examples are the scalar pressure model in chalcogenide waveguides [19], the scalar elastic model for SBS in fibers [20], isotropic beams [21], and the full-viscoelastic response of waveguides [22] at ultrasound frequencies. The mode-solvers fall under the general scheme of the semi-analytical finite element method (SAFE), which reduces the three-dimensional (3D) wave-propagation problems to 2D variants, by assuming periodicity in the direction of propagation [18]. In this paper, we develop and implement the SAFE method for fully-tensorial elastic-wave equation, described in Section III, using an open-source element-solver, called FEniCS [23].

In Section IV on isotropic materials, we perform a thorough review of various approximations of the elastic-wave equation in the SBS literature and discuss their validity in the context of integrated waveguides. As an example, a material system consisting of a chalcogenide glass, As2S3 embedded in thermal oxide is studied [24]. So far, the consensus for such isotropic materials, and following the optical fiber literature [25], has been to assume that the displacement in waveguides follows the scalar acoustic-wave model [19]. However, we find that this model overestimates the eigenmodes of the system by as much as 0.5 GHz in submicron-sized waveguides. In addition, the acoustic-wave equation assumes close to plane-wave propagation. A better approximation is provided in reference [20], which we later refer to as the scalar elastic equation. In both cases, nonetheless, these models do not appropriately capture the mode profile near the material interfaces, when compared with the most general model used in this work.

In Section V on anisotropic materials, we couple the elastic mode-solver with an electromagnetic counterpart to explore the effect of elastic anisotropies on Brillouin gain, in addition to the variation of the mechanical resonance frequency. Silicon is chosen as an example, which is ubiquitous in integrated photonics and provides high Brillouin gain in suspended structures [3, 4]. Here, researchers have resorted to full-3D simulations with the Floquet boundary conditions [26, 27] usually implemented in the commercial FEM solver, COMSOLTM. Following the micro-electro-mechanical systems (MEMS) literature [15], some authors have chosen to use the simplified isotropic model for silicon [28, 29]. While qualitatively acceptable, this assumption can lead to large discrepancies quantitatively. We show that in extreme cases, where the waveguide is aligned along the [100] or the [110] crystalline axis, the difference in eigenfrequencies can be as high as 0.8 GHz, which highlights the need for a fully-tensorial formulation. This is directly translated to the Brillouin gain, which scales with the inverse of the mechanical resonance-frequency squared. For back-scattering, we show that, for an arbitrary crystal orientation in silicon, the elastic modes do not need to have the same symmetry as the optical modes, which greatly affects the Brillouin coupling coefficient. While, in principle, the anisotropic behavior of silicon can be simulated in COMSOLTM, the use of the SAFE method in this work leads to substantial - over two orders of magnitude - computational-time improvement. Furthermore, since the SAFE method reduces the equations to a 2D problem, it requires less computational memory than the full-3D simulation, without neglecting any physical modeling features. Finally, it should be noted that recently a mode-solver for isotropic materials utilizing the finite-difference method has been presented [30], but our work provides faster convergence through the use of higher-order finite elements.

II Governing Equations

The constitutive equation of motion for an elastic medium is [17]

[TABLE]

where $\boldsymbol{u}$ is the displacement vector, $\rho$ is the density and the $\boldsymbol{\sigma}$ is the stress. Stress is linearly related to strain, $\boldsymbol{\epsilon}$, via the stiffness tensor, $\boldsymbol{C}$, as follows:

[TABLE]

To obtain the equations describing the modal distribution, the Fourier transforms of t and z are employed. We also introduce the transverse gradient operator $\nabla_{T}$ such that $\nabla=\nabla_{T}+iq_{b}\boldsymbol{\hat{z}},$ where $q_{b}$ is the acoustic-wave propagation constant. Afterwards, the elastic wave equation reads [18, 31]

[TABLE]

To obtain the weak form of Equation 3, we multiply it by a test function $\boldsymbol{v^{*}}$, integrate and apply the Green’s theorem in analogy to the 3D case [32], to obtain

[TABLE]

where the symmetry of the stiffness tensor $\boldsymbol{C}_{ijkl}=\boldsymbol{C}_{jikl}$, the notations $\boldsymbol{\epsilon}_{T}(\boldsymbol{u})=\frac{1}{2}[\nabla_{T}\boldsymbol{u}+(\nabla_{T}\boldsymbol{u})^{T}]$, as well as $\boldsymbol{\epsilon}_{Z}(\boldsymbol{u})=\frac{1}{2}[iq_{b}\boldsymbol{u}+(iq_{b}\boldsymbol{u})^{T}]$, have been utilized. Also, $\Omega$ denotes the whole space, $\partial\Omega$ the boundary and $\boldsymbol{\hat{n}}$ the normal to that boundary. Finding the elastic modes is equivalent to finding the eigenvectors of Equation 4, where $\omega^{2}$ is the eigenvalue.

III Implementation

III-A Description

FEniCS is an open-source FEM solver that automates significant portion of the finite-element assembly. The package contains the unified form language (UFL) [33], a domain-specific language for declaration of variational forms with syntax that follows mathematical notation. The core of our code consists of Equation 4, implemented in UFL. The notable difference is that UFL does not support complex numbers, thus Equation 4 is split into the real and imaginary parts that form a linear system of coupled equations. The Dolfin package [23] provides a high-level interface to various linear algebra packages that are need for efficient solution of the eigenvalue problem. The FIAT package [34] enables quick testing of different finite elements. For the tensorial elastic mode-solver, we have used the Lagrange finite elements and found them to be stable.

When the eigenvectors are expressed in terms of the basis functions of the finite elements, $\boldsymbol{u}=u_{i}\phi_{i}$ in Equation 4 becomes a generalized eigenvalue problem of the from $S(\phi_{i},\phi_{j})u_{i}=-M_{ij}\omega^{2}u_{i}$, where the eigenvalues are weighted by the mass matrix, M, in our case being the material density. After splitting the complex coefficients into the real and imaginary parts, the weak form (Equation 4) forces the stiffness matrix, S, to be symmetric, which ensures that $\omega$ is real. The stiffness matrix is sparse, since $S(\phi_{i},\phi_{j})$ is zero for nonadjacent finite elements. While solving the eigenvalue problem, we are interested only in a few modes close to the fundamental mode, therefore there is no need to find all the eigenvectors. Given these requirements, the appropriate linear algebra package is recognized to be SLEPc [35] with the Krylov-Shur algorithm and a spectral shift-and-invert preconditioner. SLEPc enables the parallel computation of eigenvalues, which is necessary for higher resolution. On a personal computer, with 16 GB of memory, we were able to compute eigenmodes for meshes up to 180 second-order Lagrange elements per facet.

Both the stress and the strain tensors are symmetric, therefore to minimize memory usage we use the Voigt notation [17], where the two tensors are mapped onto six-dimensional vectors, i.e., the stress-strain relation becomes

[TABLE]

Better numerical performance is obtained by using small numerical values, so dimensions of the parameters used in the solver are $[\rho]=$ g/cm${}^{3}=10^{-12}$ g/$\mu$m3, $[\boldsymbol{C}]=$ GPa $=10^{-12}$ g/(ns${}^{2}\mu$m) , $[x,\boldsymbol{u}]=\mu$m, which result in eigenvalues $[\omega]=$ GHz after cancellation of prefactors. It is finally reminded that the mode-solver finds two modes for each eigenvalue, corresponding to waves traveling in opposite directions [18].

III-B Testing

In the case of homogeneous isotropic materials, Equation 1 can be rewritten as $-\rho\omega^{2}\boldsymbol{u}=(2\mu+\lambda)\nabla(\nabla\cdot\boldsymbol{u})-\mu\nabla\times\nabla\times\boldsymbol{u}$. Conveniently, for shear waves $\nabla\cdot\boldsymbol{u}=0$ and the problem reduces to the Maxwell’s eigenvalue problem. If facets at $\pm a$ are clamped, then one of the solutions is $u_{y}=\cos(k_{x}x)\exp(iq_{b}z)$, with $u_{x},u_{z}=0$ and for the given boundary conditions $k_{x}=2\pi/a$. We use this fundamental mode to test the accuracy of the solver against a known solution. For the numerical simulations, we choose a 2$\mu$m $\times$ 5$\mu$m waveguide of thermal oxide and $q_{b}=11.67\mu$m*-1*.

The error is estimated as $|\omega^{2}_{analytical}-\omega^{2}_{numerical}|$. This error is plotted against the resolution of the mesh for Lagrange elements of different polynomial orders [36] in Fig. 1. As expected, for third-order Lagrange elements, the convergence rate exceeds the second-order finite-difference scheme [30].

IV Isotropic materials:chalcogenides

As the first example, the material parameters of a chalcogenide glass composition, As2S3, embedded in thermal oxide are implemented [24]. Both materials are isotropic, hence the non-zero components of their stiffness tensor elements are $c_{11}=c_{22}=c_{33}=2\mu+\lambda$, $c_{44}=c_{55}=c_{66}=\mu$. Also, the non-diagonal elements $c_{12},c_{13},c_{23}=\lambda$ together with their mirrored components. For thermal oxide, $\mu=29.9$ GPa, $\lambda=$15.4 GPa and $\rho=2.2$ g/cm3 [37], while As2S3 is expectedly softer with $\mu=6.2$ GPa, $\lambda=9.78$ GPa and $\rho=3.2$ g/cm3 [24, 38].

In the below comparison, the focus is on the fundamental longitudinal mode, $\boldsymbol{u}\approx u_{z}$. Following the optical fiber literature [25], the consensus so far has been to assume that the displacement in waveguides follows the scalar acoustic model [19]

[TABLE]

where the longitudinal velocity is related by $\rho v_{l}^{2}=2\mu+\lambda$ to the stiffness constants, and $\tilde{\rho}$ refers to the change of material density. Since the displacement is only in the z direction, $\tilde{\rho}$ is proportional to $\rho u_{z}$, which implies that Equation 6 can be rewritten as $(2\mu+\lambda)\nabla^{2}_{T}u_{z}-q_{b}^{2}(2\mu+\lambda)u_{z}=-\rho\omega^{2}u_{z}$. The eigenvalues of this and preceding equation are indeed identical. The issue with the acoustic model, however, is that it was originally developed for fluids and therefore neglects shear waves, or components of thereof. A better approach is to start with Equation 5 and enforce $\boldsymbol{u}=u_{z}$, upon which we arrive at

[TABLE]

which is known as the scalar elastic equation and has been previously used to model SBS in fibers [20].

Figure 2(a) presents the fundamental longitudinal eigenvalues of chalcogenide waveguides, based on the three models (namely, Equations 3, 6 and 7) on the same mesh with second-order Lagrange elements. In all cases, the acoustic propagation constant is fixed to $q_{b}=18.23$ $\mu$m*-1*, hence the eigenfrequency increases with mode confinement. The figure clearly shows that the scalar acoustic model is too simplistic for integrated devices with a discrepancy of 550 MHz for waveguides that are 0.5 $\mu$m wide. However, the scalar elastic equation remains accurate to within 100 MHz of the tensorial model for waveguides wider than 1 $\mu$m. The modal cross-sections normalized to unity are shown in Figure 2(b), which shows a large discrepancy between the scalar models and the tensorial model near the material interfaces.

Based on the detailed model presented later in section V B, we also calculate the Brillouin gain for chalcogenide waveguides. The physical parameters used in the simulation are as follows: photoelastic constants of $p_{11}$ = 0.25 $p_{12}$ = 0.24 for As2S3, and $p_{11}$ = 0.12 $p_{12}$ = 0.27 for SiO2 [39], waveguide dimensions of 4 $\mu$m by 0.85 $\mu$m and experimental $Q_{mech}$ = 226 [24]. The phase-matching for back-scattering of the fundamental transverse-electric (TE) mode at 1544 nm dictates $q_{b}$ = 18.37 $\mu$m*-1*. Using the mode profile calculated using our solver, we obtain a gain coefficient of 322 W*-1m-1*, which agrees well with the scalar formulation, 321 W*-1m-1*, and the experimentally measured value of 311 W*-1m-1* [24]. We note that for such large waveguides, the coupling is dominated by eletrostriction and the radiation pressure is negligible, that is less than 1 W*-1m-1*.

Finally, the chalcogenide glass material system is used to highlight the computational benefits of using the SAFE method against a fully-3D simulation in COMSOLTM with Floquet boundary conditions. Figure 3 summarizes how this work provides two orders of magnitude decrease in computational time.

V Anisotropic materials:silicon

V-A Resonance frequency

Due to its high refractive index, it is anticipated that silicon should provide high Brillouin gain. Until recently, observing the gain has been elusive, because the acoustic velocity in the silica layers, typically used for cladding the waveguides, is lower than the acoustic velocity of silicon. In other words, the acoustic is leaky in standard silicon-on-insulator (SOI) waveguides. The solution to this problem is to use a suspended membrane [3] or a waveguide supported by a nanoscale strut [4]. It should be stressed that, due to the enhancement coming from radiation pressure, the strongest Brillouin gain in silicon arises from forward-propagating modes [28].

Herein, the solutions of our solver are compared against the 3D finite-element solution of COMSOLTM. A rectangular silicon waveguide (450 nm $\times$ 230 nm) suspended in air is utilized, which means that $\int_{\partial\Omega}dx=0$ in Equation 4. Also, the phase matching condition requires that the elastic wavevector be equal to the difference between the optical wavevectors, but to a good approximation the phonons oscillate in the xy-plane, thus we can look for modes with $q_{b}=0$. Due to the symmetries of the stiffness tensor of silicon, a cubic material, the following identities hold true: $c_{11}=c_{22}=c_{33}=164$ GPa, $c_{44}=c_{55}=c_{66}=79$ GPa and the non-diagonal $c_{12}=c_{13}=c_{23}=c_{21}=c_{31}=c_{32}=64$ GPa, while the rest of elements of the stiffness tensor are zero [4]. The density of silicon is 2.328 g/cm3.

The anisotropic nature of silicon necessitates the use of a fully-tensorial description of elasticity, as implemented in this work. As mentioned before and following the MEMS literature [15], some authors have chosen to use the simplified isotropic model for silicon [28]. While it is qualitatively fair, the isotropic model can lead to large discrepancies in quantitative values. In extreme cases, when the waveguide is aligned along the [100] or the [110] crystal directions, the differences in eigenfrequencies is as much as 0.8 GHz, as shown in Table I. The modal profile is included for completeness in Fig. 4(a). Similar elastic anisotropy is expected from other cubic materials, such as germanium [40].

V-B Brillouin Forces and Gain

The presented elastic-wave mode-solver can be combined with an electromagnetic mode-solver, implemented in FEniCS and described in detail in [36], to calculate the forces acting on the opttical waveguide. The electromagnetic mode-solver is based on minimizing the functional [36]

[TABLE]

which comes from the wave equation with the $E=(E_{T}+\gamma E_{z,\gamma})\exp(-\gamma z)$ ansatz. Since for propagating modes, $\gamma$ is imaginary, the transverse field is real and the z-field is imaginary. Also, $\epsilon_{r}$ and $\mu_{r}$ are assumed to be real and scalar.

There are three types of forces that contribute to the Brillouin gain [29]. First, there is bulk electrostriction, which is the dominant force in optical fibers. Given that there are two fields present in waveguides, described by $\boldsymbol{E}=\frac{1}{2}(\boldsymbol{E_{p}}\exp(ik_{p}z-\omega_{p}t)+\boldsymbol{E_{s}}\exp(ik_{s}z-\omega_{s}t))+c.c.$, the phase-matched component of the stress, $\sigma_{ij}^{ES}\exp(iq_{b}z-\omega_{mech}t)$, is given by

[TABLE]

where $p$ is the electrostriction tensor and the force is

[TABLE]

Second, there is a boundary electrostriction force at the interface of dielectrics 1 and 2, given by

[TABLE]

where $n_{j}$ is a normal pointing from material 1 to material 2. The third force, arising from radiation pressure, also appears only on the boundary and can be computed from Maxwell’s stress tensor. Again, retaining only the phase-matched components, $T_{ij}\exp(iq_{b}z-\omega_{mech}t)$, with

[TABLE]

the corresponding force from radiation pressure is given by

[TABLE]

Distribution of the forces and the shape of the elastic modes is necessary for the computation of the Brillouin gain coefficient. Since Brillouin coupling coefficient can be related to the photon generation rate through particle flux conservation [28] and the phonon generation rate is proportional to the power generated by optical forces, $\int_{\Omega}\boldsymbol{f}\cdot\partial_{t}\boldsymbol{u^{*}}dx$, the gain coefficient is expressed as [29]

[TABLE]

where $\omega_{opt}$ and $\omega_{mech}$ refer to optical and mechanical frequencies of propagating waves, $Q_{mech}$ is the mechanical quality-factor of the mode, $P_{s,p}=\frac{c_{0}}{2n_{g}}\epsilon_{0}\int_{\Omega}\boldsymbol{E_{s,p}}\cdot\boldsymbol{\epsilon_{r}E^{*}_{s,p}}dx$ are the optical powers of the pump scattered field and $P_{mech}=\frac{1}{2}\omega_{mech}^{2}\int_{\Omega}\rho\boldsymbol{u}\cdot\boldsymbol{u^{*}}dx$ is the mechanical power. The forces in the numerator are added coherently, i.e., $\boldsymbol{f}=\boldsymbol{f^{ES}}+\boldsymbol{f^{RP}}$.

The simulations in Fig. 6 are performed for a rectangular silicon waveguide (450 nm $\times$ 230 nm), suspended in air, and for the fundamental TE mode at the 1550 nm wavelength. We have assumed the photoelastic constants of p11 = -0.09, p12 = 0.017, p44 = -0.051 and $Q_{mech}$ = 249 [4]. It is important to distinguish two cases. For the forward Brillouin scattering, $\boldsymbol{E_{s}}=\boldsymbol{E_{p}}$ and $q_{b}\approx 0$, since the frequency shift is small in comparison to the optical carrier frequency. In this case (Figs. 6(a) and 6(b)), all forces are real. The boundary electrostriction (ES) is weaker (11% when comparing the maxima) and in opposite direction to the radiation pressure (RP). For the [110] direction of propagation the dominant x-components of bulk electrostriction (1684 W*-1m-1*) and radiation pressure on the boundary (151 W*-1m-1*) add up constructively resulting in a total gain of 3099 W*-1m-1*. This value is in good agreement with the experimental direct characterization of the gain yielding 3218 W*-1m-1* and indirect characterization through cross-phase modulation giving 3055 W*-1m-1* [4]. The total simulated gain for the [100] direction is 3896 W*-1m-1*. This is to be expected from the $1/\omega_{mech}^{2}$ dependence of gain in Equation 14 and the aforementioned difference of 0.8 GHz in resonance frequency between the two directions of propagation. Evidently, using the described fully-tensorial model is critical in modeling SBS in devices constructed from anisotropic materials.

For the backward Brillouin scattering, $\boldsymbol{E_{s}}=\boldsymbol{E_{p}^{*}}$ and $q_{b}=2k_{0}n_{eff}$. In the software, the employed sign convention is negative wavevector, $-\gamma$, for the pump and positive, $+\gamma$, for the scattered field, which also implies positive wavevector, $+q_{b}$, for the elastic wave. In Fig. 6(d), the boundary electrostriction (ES) is weaker (29% when comparing the maxima) and in opposite direction to the radiation pressure (RP). Since the dominant x-component of radiation pressure in Fig. 6(d) acts in antiphase to the bulk electrostriction in Fig. 6(e), the total Brillouin gain is smaller than the electrostriction component alone, as shown in Fig. 6(f). The value of total gain for the [100] direction is 380 W*-1m-1*.

The features of backward SBS is even more complex than forward SBS and is best explained with the aid of a dispersion diagram from Fig. 5, where the resonant frequencies are plotted versus the crystal orientation with respect to the direction of propagation. For the gain coupling calculation, we follow the fundamental symmetric mode marked with a blue ”S” in Fig. 5 at 24.513 GHz and $Q_{mech}$ = 249. We note that, in general, the elastic mode is not symmetric for cubic materials such as silicon [41], unlike the optical mode. Thus, the initially symmetric mode acquires a degree of asymmetry as the crystal orientation is tilted, to the point that at 3o, the symmetric- and the antisymmetric-mode branches merge. This is the origin of the first discontinuity in the gain diagram in Fig. 6(f). The next discontinuity at 42o occurs when the common branch splits again into the symmetric and antisymmetric modes. The optical mode is symmetric, thereby the optical force distribution is symmetric, so the coupling to the asymmetric modes is null. Therefore, the Brillouin gain drops as the crystal orientation is rotated away from the [100] or [110] directions, since the elastic modes acquire a degree of asymmetry.

VI Conclusions

In summary, we have presented a fully-tensorial elastic-wave mode-solver and discussed its implementation in the open-source finite element solver, FEniCS. The source code is available at https://github.com/MarcinJM/PySBS. The use of 2D SAFE method in this work leads to computational times smaller by two orders of magnitude than the commercial 3D finite-element solver, COMSOLTM. We have also performed a thorough review and comparison of a multitude of approximations used in simulating the elastic modes in the SBS literature and compared them against the present fully-tensorial model. In some cases, the discrepancy in the calculated eigenvalues was found to be as large as 0.8 GHz. We have shown that for silicon with arbitrary crystal orientations, the elastic modes do not need to have the same symmetry as the optical modes, which greatly affects the Brillouin coupling coefficient. The modeling tool is made publicly available to other researchers and is expected to be useful in understanding and tailoring SBS in integrated devices.