COMSOL 6.3 - Stress–Optical Effects in a Photonic Waveguide

Stress–Optical Effects in a Photonic Waveguide

Planar photonic waveguides in silica (SiO2) have great potential for use in wavelength routing applications. The major problem with these kinds of waveguides is birefringence. Anisotropic refractive indices result in fundamental mode splitting and pulse broadening. The goal is to minimize birefringence effects by adapting materials and manufacturing processes. One source of birefringence is the use of a silicon (Si) wafer on which the waveguide structure is deposited. After annealing at a high temperature (approximately 1000°C), mismatch in thermal expansivity between the silica and silicon layers results in thermally induced stresses in the structure at the operating temperature (typically room temperature around 20°C).

This model requires the Wave Optics Module and the Structural Mechanics Module.

Model Definition

The Stress–Optical Effect

The general linear stress–optical relation can be written, using tensor notation, as

where Δnij = nij − n0 Iij, nij is the refractive index tensor, n0 is the refractive index for a stress-free material, Iij is the identity tensor, Bijkl is the stress–optical tensor, and Skl is the stress tensor. The number of independent parameters in the stress–optical tensor that characterizes this constitutive relation is reduced by symmetry. Because nij and Skl are both symmetric, Bijkl = Bjikl and Bijkl = Bijlk. In many cases it is possible to further reduce the number of independent parameters. The model at hand considers only two independent parameters, B1 and B2. The stress–optical relation then simplifies to

where nx = n11, ny = n22, nz = n33, Sx = S11, Sy = S22, and Sz = S33.

This translates to

Using the two parameters B1 and B2, the model assumes that the nondiagonal parts of nij are negligible.

The first part of this model utilizes the plane strain approximation available in the Solid Mechanics interface. The resulting birefringent refractive index is computed using expression variables and can be considered a postprocessing step of the plane strain model. The refractive index tensor is used as material data for the second part of the model, the mode analysis.

The Application Library model Stress–Optical Effects with Generalized Plane Strain (the model name is stress_optical_generalized) demonstrates a computation, for which the structure is free to bend in the z direction, using a formulation called Generalized Plane Strain.

Electromagnetic Waves

For a given frequency ν, or equivalently, free-space wavelength λ0 = c0/ν, the Wave Optics Module’s Electromagnetic Waves, Frequency Domain interface can be used for the mode analysis. In this model the free-space wavelength is 1.55 μm. The simulation is set up with the electric field components E = (Ex, Ey, Ez) as dependent variables. The wave is assumed to have the form

and the effective mode index, neff = β/k0, is obtained from the eigenvalues.

For propagating modes it must hold that

The optical core and planar waveguide layers are made of Silica (SiO2) which is deposited unto a silicon (Si) wafer. The material properties are shown in the following table:


Name	Expression	Description
nSi	3.5	Refractive index, silicon (Si)
nSiO2	1.445	Refractive index, silica (SiO2)
deltan	0.0075	Relative index difference:
nCore	nClad/sqrt(1-2*deltan)	Refractive index, core
alphaSi	2.5e-6[1/K]	Coefficient of thermal expansion Si
alphaSiO2	0.35e-6[1/K]	Coefficient of thermal expansion, SiO2
ESi	110[GPa]	Young’s modulus, Si
ESiO2	78[GPa]	Young’s modulus, SiO2
nuSi	0.19	Poisson’s ratio, Si
nuSiO2	0.17	Poisson’s ratio, SiO2
B1	0.65e-12[m^2/N]	First stress optical coefficient
B2	4.2e-12[m^2/N]	Second stress optical coefficient
T1	20[degC]	Operating temperature
T0	1000[degC]	Reference temperature

The computational domain can be reduced significantly for the optical mode analysis because the energy of the fundamental modes is concentrated in the core region and the energy density decays rapidly in the cladding and buffer regions. Thus, you do not need to model the air domain.

Results and Discussion

Figure 1 shows the von Mises stress distribution together with the deformed shape of the waveguide. Notice that the stress varies slowly in the horizontal direction. This means that the significant influence on the stress-induced changes in the refractive index comes from the stress variations in the vertical direction. This is expected because the extension of the domains in the x direction is chosen to minimize the effects of the edges.

Figure 1: von Mises equivalent stress.

Figure 2 shows the birefringence, nx − ny, in the optical computation domain. The birefringence is nearly constant along any horizontal line, so the influence of the side edges is indeed reduced to a minimum.

Figure 2: Birefringence.

Figure 3 shows the effective mode indices for the first four propagating modes. The value para = 0 corresponds to the case without stress–optical effect, while this effect is included for para = 1.

The difference is significant, which shows that the shift in the effective mode indices due to the stress–optical effect is indeed resolved.

The computations show a shift in effective mode index due to the stress-induced change in refractive index. The birefringence causes the otherwise two-fold degenerate fundamental mode to split.

Figure 3: Effective mode indices for the first four propagating modes without (para = 0) and with (para = 1) the stress–optical effect.

Figure 4 shows the out-of-plane component (z component) of the Poynting vector for four eigenmodes computed without the stress–optical coupling.

Figure 4: Eigenmodes computed without the stress–optical coupling.

Visual inspection of the higher eigenmodes indicates that they have a larger portion of energy leaking into the cladding and buffer, and are thus more affected than the fundamental modes by the distance to the air and silicon layers. Because of this leakage, the boundary condition affects the higher eigenmodes more than the fundamental mode.

The change in the eigenmodes shape due to the stress–optical effect is shown in Figure 5.

Figure 5: Eigenmodes computed with the stress–optical effect taken into account.

Although the fundamental modes have converged to 5 decimal places, the known modeling errors make the exactness of the numbers uncertain. One major modeling error is due to the fact that the model contains a plane strain assumption in a case where the real-world model does not necessarily conform to this deformation state. This modeling error is reduced in the refined model Stress–Optical Effects with Generalized Plane Strain (the model name is stress_optical_generalized).

References

1. H.P. Schriemer and M. Cada, “Modal birefringence and power density distribution in strained buried-core square waveguides,” IEEE Journal of Quantum Electronics, vol. 40, pp. 1131–1139, 2004.

2. J. Stone, “Stress-Optic Effects, Birefringence, and Reduction of Birefringence by Annealing in Fiber Fabry-Perot Interferometers,” J. Lightwave Technol., vol. 6, pp. 1245–1248, 1988.

Wave_Optics_Module/Waveguides/stress_optical

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file stress_optical_parameters.txt.

Add a parameter to switch the stress effects on and off.

In the table, enter the following settings:


Name	Expression	Value	Description
para	1	1	1: stress-optical coupling, 0: no coupling

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 320.

In the text field, type 83.

Locate the section. In the text field, type -160.

In the text field, type -100.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 320.

In the text field, type 14.

Locate the section. In the text field, type -160.

In the text field, type -17.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 320.

In the text field, type 16.

Locate the section. In the text field, type -160.

In the text field, type -3.

Rectangle 4 (r4)

In the toolbar, click

In the window for , locate the section.

In the text field, type 6.

Locate the section. In the text field, type -3.

In the text field, type -3.

Rectangle 5 (r5)

In the toolbar, click

In the window for , locate the section.

In the text field, type 20.

Locate the section. In the text field, type -10.

In the text field, type -10.

Click

The last rectangular region encloses the optical computational domain. It can be enlarged if needed for validating the results. The region should be chosen large enough so that the computed propagation constants do not change significantly if the region is enlarged.

Solid Mechanics (solid)

Linear Elastic Material 1

In the window, under > click .

Thermal Expansion 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the T list, choose . In the associated text field, type T1.

Click

for .

Global Definitions

Default Model Inputs

In the window, under click .

In the window for , locate the section.

Find the subsection. In the text field, type T0.

Definitions

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
N	nCore		Refractive index for stress-free material

Locate the section. From the list, choose .

Select Domain 6 only.

Variables 2

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
N	nSiO2		Refractive index for stress-free material

Locate the section. From the list, choose .

Select Domains 4 and 5 only.

Variables 3

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Nx	N-para(B1solid.sx+B2*(solid.sy+solid.sz))		Refractive index, x-component
Ny	N-para(B1solid.sy+B2*(solid.sx+solid.sz))		Refractive index, y-component
Nz	N-para(B1solid.sz+B2*(solid.sx+solid.sy))		Refractive index, z-component

Locate the section. From the list, choose .

Select Domains 4–6 only.

Solid Mechanics (solid)

All regions have free boundaries, which also is the default boundary condition. However, these conditions will not be sufficient for creating a unique solution because the computational domain is allowed to move and rotate freely. The problem becomes well-posed by adding constraints to restrain such rigid body movements.

Rigid Motion Suppression 1

In the toolbar, click

and choose .

Select Domain 1 only.

Electromagnetic Waves, Frequency Domain (ewfd)

The computational domain is reduced significantly for the optical mode analysis.

In the window, under click .

Select Domains 4–6 only.

Wave Equation, Electric 1

In the window, under > click .

In the window for , locate the section.

From the n list, choose . From the list, choose .

Specify the n matrix as


Nx	0	0
0	Ny	0
0	0	Nz

From the k list, choose .

Materials

In the window, under right-click and choose .

In the window for , type Si in the text field.

Select Domain 1 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	ESi	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	nuSi	1	Young’s modulus and Poisson’s ratio
Density	rho	rhoSi	kg/m³	Basic
Coefficient of thermal expansion	alpha_iso ; alphaii = alpha_iso, alphaij = 0	alphaSi	1/K	Basic

SiO2

Right-click and choose .

In the window for , type SiO2 in the text field.

Select Domains 2–6 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	ESiO2	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	nuSiO2	1	Young’s modulus and Poisson’s ratio
Density	rho	rhoSiO2	kg/m³	Basic
Coefficient of thermal expansion	alpha_iso ; alphaii = alpha_iso, alphaij = 0	alphaSiO2	1/K	Basic

Mesh 1

Free Triangular 1

In the toolbar, click

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domains 4–6 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 0.2.

Size

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Click

Study 1

Step 2: Mode Analysis

In the toolbar, click

and choose > .

In the window for , locate the section.

From the list, choose .

In the text field, type 1.46.

Select the checkbox. In the associated text field, type 4.

In the text field, type c_const/lambda0_ewfd.

These settings make the eigenmode solver search for the 4 eigenmodes with effective mode indices closest to the value 1.46. This value is an estimate of the effective mode index for the fundamental mode.

Exclude Solid Mechanics from the Mode Analysis step.

Locate the section. In the column of the table, under , clear the checkbox for .

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
para (1: stress-optical coupling, 0: no coupling)	0 1

In the toolbar, click

Results

Stress (solid)

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Click the

button in the toolbar.

To visualize the details of the eigenmodes, you first set up a view that includes the optical computation domain only.

Click the

button in the toolbar.

In the dialog, in the tree, select the checkbox for the node > .

Click .

View 2D 2

In the window, under right-click and choose .

In the window for , locate the section.

Select the checkbox.

Axis

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type -11.

In the text field, type 11.

In the text field, type -11.

In the text field, type 11.

Electric Field (ewfd)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Surface 1

In the window, expand the node, then click .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > > .

Electric Field (ewfd)

In the window, click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

. This creates a visualization of the power flow, also called optical intensity or the Poynting vector, in the z direction (out-of-plane direction). Continue to plot the power flow for the other computed eigenmodes, with and without stress effects.

2D Plot Group 3

Next, plot the birefringence Nx-Ny in the waveguide.

In the toolbar, click