Thin Film Lithium Niobate Electro-Optic Phase Modulator

In this example, we simulate a thin film Lithium Niobate (LN) phase modulator based on the work of Mercante et al . [1]. The RF induced capacitive electric fields (E-fields) are calculated in CHARGE taking advantage of the anisotropic DC dielectric permittivity feature introduced in 2023 R1.2. The electric fields are then used to calculate the electro-optic index perturbation in LN via the Pockels effect at telecom wavelengths. Optical modes of the perturbed LN waveguide are then calculated in FEEM, and the voltage dependent phase modulation metrics of the fundamental TE mode are calculated including loss, \(L_{\pi}\), and \(V_{\pi}L\).

download example

Overview

Optical transceivers convert electrical signals into the optical domain. All computation begins and ends electronically, but by translating from the electrical to optical domain we open up more channels, with greater bandwidth, that can be passed long distances with significantly less attenuation. These devices play essential role for long haul segments of the internet and are increasingly moving upstream in the network as traffic and latency requirements grow. It is helpful to think of the optical transceiver and its complementary device the photodetector as the on and off ramps to the superhighways of the internet.

These devices typically use phase modulators in a Mach-Zehnder Interferometer configuration with the carrier wave split into two arms and recombined at the output. Changing the phase of light in the arms in an unbalanced way, via electrical signals, results in constructive or destructive interference at the output. MZIs often serve as very sensitive optical instruments, but in this capacity the phase of the light is intentionally modulated and they are referred to as Mach-Zehnder Modulators (MZM) in this capacity. Several material platforms and physical effects have been used to achieve this functionality. In this example we focus on Pockel's effect in LN.

Most physical mechanisms for phase modulation are quite weak resulting in very large devices. Alternatively a number of exotic materials that are lossy and/or challenging to integrate with other optical and electronic components. LN is transparent over a large bandwidth and has large anisotropies, allowing for low loss and high modulation efficiencies. Conventional bulk-LN has been used widely; however, traditional fabrication methods did not allow for high index contrast waveguides. Recent developments in fabrication techniques have made the thin-film LNOI platform an excellent candidate for ultracompact and high-performance integrated photonic components.

In this article we demonstrate how to simulate the electro-optic modulation in LNOI using our Finite Element IDE. The simulations performed as part of this work consist of two main stages: i. Electrical and ii. Optical. The schematic of the simulated modulator is shown below.

Step 1: Electrical Simulation

In Step 1, we use the CHARGE solver to simulate the electric field distribution in the LN rib waveguide resulting from applied voltage bias. The bias is applied through gold electrodes in ground-signal-ground configuration.  Voltages ranging from 0 V to 5 V in intervals of 0.5 V are applied to the signal electrode with a fixed voltage of 0 V applied to the ground electrodes. The electric field results are used to calculate the index perturbation in the LN material through the Pockels effect.

Step 2: Optical Simulation

From the index perturbation calculations performed in Step 1 a nk material simulation object is created and apply to the LN waveguide structure. Then the FEEM solver is used to calculate the modes of the waveguide at a wavelength of 1.55 um. These actions are performed within a for loop, where each iteration corresponds to a single voltage point. We track the fundamental TE mode and plot the variation of the effective index with applied voltage. We also calculate the associated losses (in dB/cm) and the voltage-length product \(V_{\pi} L\) for different voltages.

Run and Results

Instructions for running the model and discussion of key results

Step 1: Electrical Simulation

1. Open [[LN_phase_modulator.ldev]] using finite Element IDE.
2. Open [[LN_phase_modulator_CHARGE.lsf]] and run the script.
3. The script performs the following actions:
   • Runs the CHARGE simulation and collects the electrostatics results from the monitor.
   • Calculates the index perturbation using Pockels effect.
   • Opens the visualizer to view the calculated results such as electric fields, n, dn stored in an unstructured dataset “electro”.

Electrostatic Results

The dataset electrostatics provide a number of values, including the important results of the CHARGE simulation, the E-field (\(E\)) between the capacitor plates. The open visualizer displays these and allows us to inspect them. Clicking through the V_signal parameters, figure below, we can see how these results change as a function of Voltage.

CHARGE uses unstructured datasets which may be unfamiliar to some users. The area and ID are properties of the constituent triangles of the mesh and may be useful for further calculations. For more information on how to work with unstructured datasets see - Dataset builder utility overview . Here we simply used \(E\) on the mesh generated by CHARGE to calculate our index results on all vertices.

EO Index Perturbation

Using the \(E\) values we calculate n_EO and dn the total spatial vector indices, and the difference from the base indices respectively. The figure below graphically displays the dn results of a voltage sweep from 0 to 5 V with 0.5 V spacing, which provides insight into the devices operation, but is not used in the next step.

To view and compare dn , move the slider to 1.5V and then lock the color bar. Then we can move the slider in either direction to get a sense of the index change induced by RF fields. It is important to note that the fields and the index perturbation calculations are only done in the LN waveguide structure.

Step 2: Optical Simulation

1. Open [[LN_phase_modulator_FEEM.lsf]] and run.
2. The script performs the following actions:
   • Adds the anisotropic index perturbation to nk simulation object.
   • Calculates the modes at each voltage and finds the fundamental TE mode.
   • Obtains the effective index of the tracked mode and calculates the loss and \(V_{\pi}L\), a common modulation metric.

Optical Modes

The fundamental TE mode of the LN rib waveguide at a voltage of 5 V are shown above. The mode profile of the fundamental TE mode shows evanescent tails extending underneath the metal electrodes. It should be noted that the  \(x\)-cut LN is susceptible to mode hybridization as the modal plane experiences two different refractive indices [2]. Furthermore, it is also likely that the TE and TM modes of the waveguide hybridize with plasmon modes at the dielectric-Au electrodes interfaces creating these long evanescent tails. The presence of numerous modes around \(n_o\) is another observation indicative of mode hybridization. Therefore, we use a TE polarization fraction of 95% to carefully identify the fundamental TE mode.

Modulation Metrics

The effective index of the fundamental TE mode is directly collected from the FEEM results during the voltage sweep. The linear dependence of the effective index on the applied voltage, characteristic of Pockels effect is shown below.

The change in the real part of the effective index with respect to the real part of index at 0 V gives us \(L_{\pi}\), which is the required modulator length to create a phase change of \(\pi\). This can be multiplied by the voltage to get \(V_{\pi}L\), a crucial metric for measuring modulator performance. The exact expression used for the calculation of \(L_{\pi}\) is:

$$L_{\pi} = \displaystyle\frac{\lambda}{2\times \Re(n_{\mathrm{eff}}(0 V) - n_{\mathrm{eff}}( \mathrm{V}))}$$

The plot for \(V_{\pi}L\) vs Voltage is shown below. We obtain a constant \(V_{\pi}L\) value of 3.18 V-cm. In comparison, the cited literature [1] reported a DC half-wave voltage of 3.8 V per cm. This difference can be due to incomplete material properties reported in the literature, process-specific material properties of LN, and experimental non-idealities.

The losses (in dB/cm) suffered by the fundamental TE mode of the device is calculated using the following equation.

$$\mathrm{Loss} = -0.20*\mathrm{log}_{10}(e^{-2\pi\kappa/\lambda})$$

Where, \(\kappa\) is the imaginary part of the effective index of the fundamental TE mode. For more information on the relationship between the imaginary part of the effective index and the material losses, please refer: Working with lossy modes and dB/m to \(\kappa\) conversion .

The losses in the modulator are almost constant over the voltages as shown above. The losses can be be mostly attributed to the Au electrodes. As a result, the losses do not depend upon the applied voltage unlike the other modulator metrics. The key factor influencing the losses is the separation between the electrodes which can be optimized for better performance to loss ratio.

Important Model Settings

Description of important objects and settings used in this model

Simulation Region Dimension and Boundary Conditions: CHARGE

A 2D \(y\)-normal simulation region is used for the CHARGE and FEEM simulations.  At the external (open) boundaries of the simulation domain, homogeneous Neumann boundary conditions are applied. This corresponds to the normal component of the \(E\) field being set to zero at the edge of the simulation; which is physically interpreted as completely insulating box. Symmetry considerations ensure that this device will have no \(E\) flux through the line \(x = 0\). The other open boundaries of the simulation region do not have this same symmetry, and to obtain a finite solver domain we need to truncate the fields. This will have the effect of confining the fields and will affect the overall profile but should not strongly affect the results in the the waveguide core. The choice of \(x\)-min and \(z\) span is therefore purely based on the simulation span required to attain convergence of the DC \(E\) field in the waveguide core, which will be discussed in a subsequent section. A z-span of 9 \(\mu\)m is used based on convergence testing performed in CHARGE.

From the field distributions shown in the above figures, it is observed that the \(E_x\)-component of the electric field is much larger than the \(E_z\) component due to the invariance in the propagation direction. The device is designed like a parallel plate capacitor, with a constant \(E\)-field in the region of the optical waveguide core. The \(E_x\) component therefore plays the most critical role in the determination of index perturbation. The crystal orientation of the LN is chosen such that the \(r_{33}\) component of the \(\chi_2\) tensor aligns with the \(E\)-field in this direction.

Electrical Boundary conditions are applied to the “Signal” and “Ground” electrodes. This device operates with a central electrode modulated with a positive bias and two outer grounded electrodes. As mentioned before, a voltage sweep ranging from (0-5) V in intervals of 0.5 V is applied to the signal (central) electrode and a voltage of 0 V is applied to the ground (left) electrode.

Materials: CHARGE

The CHARGE Poisson solver, includes the diagonal anisotropy of DC Permittivity which was added to the CHARGE solver in the 2023 R1.2 release. In the material list, one can see that the \(x/y\)-cut and the \(z\)-cut LN are both available as insulator and semiconductor type materials. The \(x\)-cut LN used in this work has an extraordinary permittivity of 27.9 along the \(x\)-axis and an ordinary permittivity of 44.3 along the \(y\) and \(z\) axis.

When setting materials as insulators they are treated as boundary conditions in CHARGE simulations. This is helpful when we want to look at the charge transport; however, we are mostly interested in the electrostatics; therefore, we select the semiconductor type material. The Poisson equation in electrostatics is only solved in semiconductor regions. The large bandgap of LN and negligible carrier concentration ensure that there will not be charge transport and only the electrostatic calculations will be used. Edge effects around the electrode's sharp corners amplify the field intensity. In realistic process these will be rounded and we would not see such high peak fields.  This should not be a significant concern as the RF fields in the LN optical core of the waveguide will have the greatest contribution the TE mode effective index.

nk material Simulation Object

The index perturbation in LN along the \(x\), \(y\), and \(z\) directions due to the applied voltage is incorporated into the modal calculations using a nk material object. We create an unstructured dataset called “nkmaterial” using the mesh from the electrostatics results in the CHARGE. The perturbed index n_EO is stored as an attribute for different voltages and wavelength(s) (parameters). We import the dataset and apply its attribute “nk” to the LN waveguide overriding the underlying optical material properties.

Simulation Region Dimension and Boundary Conditions: FEEM

In FEEM, the simulation region should be large enough to have the electric fields decayed sufficiently at the boundaries. This limits the mode interaction with the PML or PEC. Convergence testing on the \(x\)-span and the \(z\)-span of the simulation region was performed by varying the z-span from 5 \(\mu\)m to 13 \(\mu\)m. Due to the presence of metal electrodes, a much larger range of values is chosen for \(x\)-span, 10 \(\mu\)m to 30 \(\mu\)m. The effective index and the losses of the fundamental TE mode are tracked over these values. Based on the results, an \(x\)-span of 28 \(\mu\)m and a \(z\)-span of 9 \(\mu\)m are used.

We also observe that the type of boundary (PML/PEC) has negligible effect on the properties of the fundamental TE mode. In this simulation, we use PML BC. However, if users wish to change any simulation parameters or look at higher order modes, they may consider performing convergence testing to confirm the accuracy of the results.

Materials: FEEM

Calculation of \(\Delta n\)

LN is a negative uniaxial crystal with an extraordinary refractive index of 2.14 and ordinary refractive index of 2.21 at a wavelength of 1550 nm. The linear electro-optic effect in LN is expressed as \(\Delta \left(\frac{1}{n^{2}}\right) = \sum_k r_{ijk}E_k\), where \(i\), \(j\), \(k\) have values 1, 2, 3 corresponding to the \(x\), \(y\), and \(z\) axis, and \(r_{ijk}\) is the Pockels tensor. From the electrostatic simulations, it is shown that \(E_y = 0\) and \(E_x >> E_z\). These observations only leave the \(E_x\) component of the electric field for consideration. The perturbed index of the LN material can then be expressed as a diagonal matrix as shown below:

$$\Delta n_{p}^2 = {\left( {\begin{bmatrix}n_e^2 & 0 & 0 \\ 0 &n_o^2 & 0\\ 0 & 0 &n_o^2 \end{bmatrix}}^{-1} + \begin{bmatrix}r_{33}E_x & 0 & 0 \\ 0 & r_{13}E_x & 0\\ 0 & 0 & r_{13}E_x \end{bmatrix} \right)}^{-1}$$

Here, \(n_p\) is the perturbed index. On further simplification, the perturbed index in each direction can be reformulated as (although we directly use the matrix equations for our calculations):

$$n_x = n_e - \frac{n_e^3}{2}r_{33}E_x$$

$$n_y = n_o - \frac{n_o^3}{2}r_{13}E_x$$

$$n_z = n_o - \frac{n_o^3}{2}r_{13}E_x$$

The calculated index perturbation is applied to the " LiNbO3 WG " object as a nk material. We also perform calculations considering the off-diagonal elements of the Pockels tensor to confirm that the diagonal approximation can be successfully applied without affecting the accuracy of the results. For more details on the diagonal approximation and the off-diagonal tensor elements, please refer to Appendix C.

Mesh

We use 4 "edges per wavelength" and a polynomial order of 2 for meshing the modulator in FEEM. The users may refer to: Graded-index fiber with FEEM to obtain more information on convergence testing of results based on the two parameters mentioned above.

Lock Mesh: Enabled

The users might note that for the FEEM simulation, the "Lock Mesh" functionality is enabled. This ensures that we use the same mesh for each iteration of the FEEM sweep. This is done to avoid slight variations in meshing from iteration to iteration. Although the loss of precision due to mesh variations is typically insignificant, the \(V_\pi L\) result is very sensitive to effective index change (\\Delta n\). Locking the mesh allows us to minimize numerical noise and produces a smooth result.

Mode Tracking

The FEEM simulation is run in a for loop with index perturbation for each voltage point in step 1 being applied. We only track the fundamental TE mode of the device by first filtering through the polarization fraction for mode with TE polarization fraction greater than 95% then taking the mode with the highest index. An alternative approach to tracking the modes would be to calculate the unperturbed modes of the device first (without the nk material) and then select the perturbed mode that has the highest overlap with the former. MODE solver has methods for doing this for various sweeps but becomes a bit more manual when changing the underlying material properties.

Mode tracking is challenging in this case due to the presence of the Au electrodes. This result in numerous plasmon modes guided in the partly etched thin film LN. Furthermore, the modes that are guided by the waveguide core show significant interaction with the electrodes resulting in increased loss. This can be seen in the Optical Modes results from the “Run and Results” section. For these reasons of the Au electrodes present challenges in convergence testing this simulation. One could simplify the analysis by choosing to disable them; however, the interaction is physical and important for valid results as they play a significant role in the loss. This highlights a major design trade-off where the proximity of the electrodes allows for more effective modulation in the RF while leading to increased loss for the optical mode.

Updating the Model With Your Parameters

Instructions for updating the model based on your device parameters

• Material: The user may change the material of any of the structures depending upon their requirements. As the CHARGE and FEEM simulations are performed in a single simulation (.ldev) file, the customer will have to ensure that every material has both their electrical and optical properties linked to them. This is done by first adding the material from the “Electrical and Thermal” database and then adding the same material from the “Optical Database” to the existing material. It is important to keep in mind that if the user wishes to change LN into a different material, they will also have to update the equations used to formulate Pockels effect and the Pockels coefficients in the script.
• Wavelength and crystal orientation: The wavelength of the FEEM simulation can be changed from the “Modal Analysis” tab of FEEM. The users must make sure that they update the index/permittivity of dispersive materials for more accurate results. Any change in the crystal orientation will require the users to apply this rotation to the formulation of Pockels effect. However, if the change in orientation results in off-diagonal elements in the index perturbation matrix then users must use FDE (instead of FEEM) for more accurate results. This is because FEEM does not support the “Matrix Transform” grid attribute which allows the users in MODE and FDTD to input the full transformation matrix.
• Electrode Spacing: The electrode spacing is a key structural parameter in modulator design that can be optimized to balance the modulator performance and the losses. The users can change the distance directly by editing the geometry of the individual electrodes.

Taking the Model Further

Information and tips for users that want to further customize the model

Optical Design

The interaction between the optical mode and the electrodes creates some plasmonic characteristics in the mode, resulting in greater loss. It would be possible to move the electrodes further away to minimize loss, but this would reduce the modulation efficiency. Balancing these metrics is a critical part of the phase modulator design.

RF Analysis

We have assumed that RF operation is essentially static compared to the optical modes. This is reasonable as a first step in order to focus on the optical performance of the device. But the non-idealities in the RF operation should be considered as well, using a travelling wave model. Ansys HFSS can be used to perform the RF analysis and extract the S-parameters of the interaction region. It can also be used to evaluate the capacitance of the modulator, which affects its time-constant/switching time.

Circuit Level Validation

Simulations using INTERCONNECT incorporate OOK (On-off keying) at the modulation efficiency given that both arms will be modulated. It is also important to evaluate how the losses affect the bit error rate down the line.

Additional Resources

Additional documentation, examples and training material

Related Publications

1. Mercante, Andrew J., et al. "Thin film lithium niobate electro-optic modulator with terahertz operating bandwidth." Optics express 26.11 (2018): 14810-14816.
2. Wang, Jingyi, et al. "Polarization coupling of \(X\)-cut thin film lithium niobate based waveguides." IEEE Photonics Journal 12.3 (2020): 1-10.
3. Zhu, Di, et al. "Integrated photonics on thin-film lithium niobate." Advances in Optics and Photonics 13.2 (2021): 242-352.

See Also

• CHARGE solver introduction
• Spatial (n,k) data - Simulation object
• Finite Element EigenMode (FEEM) solver introduction
• Matrix Transformation - Simulation object
• Ferroelectric modulator
• Thermally tuned waveguide (FEEM)
• Thermally tuned waveguide (FDE)
• Working with lossy modes and dB/m to κ conversion
• Graded-index fiber with FEEM

Appendix

Additional background information and theory

Pockels Effect: Full Tensor Calculation

The Pockels tensor of LiNbO3 also consists of off-diagonal terms. In the presence of electric fields, the change in the inverse permittivity \(\frac{1}{n^2}\) (impermeability tensor) is given as:

$$ \left(\frac{1}{n^2}\right) = {\left( {\begin{bmatrix}n_e^2 & 0 & 0 \\ 0 &n_o^2 & 0\\ 0 & 0 &n_o^2 \end{bmatrix}}^{-1} + \begin{bmatrix}r_{33}E_x & r_{42}E_y & r_{42}E_z \\ r_{42}E_y & r_{13}E_x & -r_{22}E_z \\  r_{42}E_z & -r_{22}E_z & -r_{22}E_y + r_{13}E_x \end{bmatrix} \right)}$$

The perturbed refractive index is therefore given as:

$$ n_{p}^2 = {\left( {\begin{bmatrix}n_e^2 & 0 & 0 \\ 0 &n_o^2 & 0\\ 0 & 0 &n_o^2 \end{bmatrix}}^{-1} + \begin{bmatrix}r_{33}E_x & r_{42}E_y & r_{42}E_z \\ r_{42}E_y & r_{13}E_x & -r_{22}E_z \\  r_{42}E_z & -r_{22}E_z & -r_{22}E_y + r_{13}E_x \end{bmatrix} \right)}^{-1}$$

The comparison (on log scale) between the diagonal and the off-diagonal components of the permittivity tensor at three different locations is shown below. As the values of the diagonal components of the tensor are close on the log-scale, they cluster together at the top of the graph. The off-diagonal components of the tensor in the waveguide region under the signal and ground electrodes are at least 4 orders of magnitude smaller than the diagonal components. They are at least 7 orders of magnitude lower at the center of the waveguide.

To further validate our diagonal approximation, we consider the full permittivity tensor and diagonalize it to determine the index along the \(x\), \(y\), and \(z\) directions. We then calculate the relative percentage error with respect to the values under the diagonal approximation and plot it on the log scale as shown below.

The error is extremely small with the maximum error being 4.6e-11% along the \(x\)-direction and 1.55e-06% along the \(y\)- direction, and \(z\)-direction at the center of the waveguide. Therefore, the rotation imparted by the off-diagonal elements will be negligible and can be safely ignored. If the user wishes to consider the off-diagonal elements, they are encouraged to look at our Application Gallery Example: Ferroelectric Modulator. As mentioned before, FEEM does not offer the Matrix Transformation grid attribute. The users will be required to use FDE for their calculations.