Lithium Niobate Nonlinear Thermal Waveguide (FEEM)

Second-harmonic generation (SHG) in a Lithium Niobite - LiNbO3 (LNO) nanophotonic waveguide is studied using thermal tuning to achieve efficient phase matching. An analytic anisotropic material model allows us to sweep over temperature and wavelength, calculating the effective index of the fundamental and 2 ^nd harmonic mode using FEEM. In these complementary spectral domains, the dispersion curves cross satisfying a Type 1 phase-matching condition. The LNO indices have distinct temperature dependencies allowing the phase matching to be thermally tuned.

download example

Overview

Understand the simulation workflow and key results

Many integrated photonics applications require multiple coherent tunable sources. Introducing extra optical inputs from external lasers adds to the complexity of chip packaging; adding integrated active sources using III-V, SiN or LNO greatly increases fabrication complexity and cost. Sources generated through nonlinear harmonic interaction which do not require additional fabrication steps to produce gain, feedback, or external electric circuitry could support this essential functionality. In the following example, we demonstrate an LNO waveguide frequency conversion, based on Luo [1].

Step 1

Run initial simulation at the center wavelength and temperature to find the reference modes.

Step 2

Sweep over both spectral bands and across the operating temperature range

Step 3

Fit the data points to a function and determine the intercept of the dispersion curves at each temperature. The resulting temperature-dependent tuning slope of the phase matching is plotted as well as the conversion efficiency.

Run and results

Instructions for running the model and discussion of key results

Step 1: Material Properties and Modes

1. [[step]]Open [[linbo03_wg.ldevs]].
2. Run the script file [[waveguide_modesFEEM.lsf ]]
3. (Optional) run [[LNO_material.lsf ]] to save material data to text file

To perform the sweep in the next step we should look at the reference modes. To do this the script [[waveguide_modesFEEM.lsf]] is set up to solve for the fundamental TE0 mode and the second-harmonic TM2 mode at the center wavelengths and the mid-point operating temperature. The mode E-fields are then visualized for the user to inspect.

Step 2: Temperature, Wavelength Sweep

1. [[step]] (optional) Change your resources to solve sweep in parallel
2. Run the python file [[FEEM_TL_sweep.py]]

Wavelength and temperature range/resolution are set in at the top of the script. We have chosen 3 wavelength points and 3 temperatures [0,100C] for each spectral band. To ensure that the analysis is valid the wavelengths are defined such that lam_nir*2 = lam_tel. Once sweeps have been completed for both the fundamental and second harmonic the python file will call the script neff_fit.lsf. This file will fit the dispersion points to a 2 ^nd order polynomial, allowing us to readily find the intersection. Once this is calculated the curves at each temperature are plot and the results are saved to a json file.

Step 3: Phase matching and conversion efficiency

1. [[step]] Run [[post_plot.py]].

Simply run post plot from the script file editor and the bundled version of Python will be launched from a new terminal. The following script loads the json file written in the previous step and reshapes the variable back to their original form. Once we have the polynomial coefficients fit in step 2 we can find the roots by solving for the difference using the quadratic equation. Next the intercepts are plot and the linear tuning slope is calculated. Finally, the wavelength-dependent SHG efficiency is calculated and plotted in Python.

Important model settings

Description of important objects and settings used in this model

Step 1

Simulation region dimensions :

The solver region must be large enough to completely contain the mode. This includes the evanescent tail of the fields to a reasonable level. A good rule is to ensure the field has decayed to at least 1e-4 the peak intensity at the boundaries. This is easy to visualize by plotting the fields on a log scale. Higher accuracy tolerance will require lower levels and hence a larger simulation region.

FEEM Mesh: As a finite-element based solver, FEEM has two important global parameters to influence the mesh and therefore the accuracy of the solution: “edges per wavelength” and “polynomial order”. Since we expect the modes to be smooth and to not contain any sharp features or singularities, it makes sense to use a coarser mesh with a higher polynomial order. We start with “edges per wavelength” of just 1 and a “polynomial order” of 3. To check if our results are converged, we can repeat the calculation with “polynomial order” 4.

Boundary conditions :

If the fields are negligible at the boundary, as we ensured with simulation region dimension then the choice of using PML or PEC boundary should be significant. Metal boundaries are faster and minimize the simulation memory required. This allows us to sweep through temperature and wavelength more quickly. PML boundaries in FEEM are a new feature as of 2021R1.4, and will typically give more accurate results, particularly with lossy modes.

Material

The functional materials for LNO have been derived from references [2], [3], which are defined in the model set-up script. These functions have the following form.

$$n_{*}(\lambda, T)=a_{*}(\lambda)\left(T-T_{0}\right)+b_{*}(\lambda)\left(T^{2}-T_{0}^{2}\right)+S(\lambda)$$

Here S is based on experimental data of the wavelength dispersion from [2] fit to a 3 oscillator Sellemier equation. Thermo-optic coefficients a, and b are interpolated values from [3] that account for wavelength dispersion. See Appendix for more details.

As we can see the material indices of LiNbO3 at the center wavelength of 1550nm differs with its 2nd harmonic by ~0.04-0.05 or 2%. This results in a coherence length of \( L_{c,o} \) ~ 16um and \( L_{c,e} \) ~19um for ordinary and extraordinarily polarized beams. Over this length scale the pump beam and it’s second harmonic will propagate out of phase by a factor of \(\pi\). Interaction regions longer than the coherence length are counterproductive to SHG, but by engineering the modes such that the effective indices are almost equal it is possible to achieve phase matching and produce efficient devices of arbitrary length.

Across the operating temperature range, 0-100C, we can see that the temperature variation of the extraordinary index is about 0.2% greater than the ordinary index, see in the figure below. This will prove to be enough to produce tuning across most of the spectral bands.

Modes

The effective index will be highest for the fundamental mode and decrease with ascending order, due to less confinement in the higher index core region. By looking at higher order modes of the second harmonic SH, we should be able to identify a mode with a similar effective index to the fundamental mode at telecom. The geometric contributions to the effective index will have less importance at the second harmonic since the waveguide is twice the size in units of lambda. Furthermore, the modes co-propagate, but the anisotropy of the crystal means that differently polarized modes will experience different material indices. We can take advantage of different mode polarizations by aligning the SH mode along the extraordinary axis. LNO is a negative uniaxial crystal, so the 1.55um mode will be polarized along the ordinary axis(blue) and the 0.775 should primarily be polarized along the extraordinary axis(yellow) in lithium niobate index figure above.

As in [1] the modes of interest are the TM2 \(_{NIR}\) and TE0 \( _{Tel}\) since these have comparable effective indices for this particular waveguide geometry and crystal orientation. After running [[waveguide_modesFEEM.lsf]] to solve for modes at the center wavelength and temperature of the telecom band we identify the TE0 as mode1 and at then perform the same analysis at the second harmonic to identify the TM2 as mode7 below.

In the visualizer, one can look at the polarization, to confirm this.

Step 2:

Reference modes :

Now that we know the mode number for 1550nm and 775nm to be 1, and 7 we set this in the python file. Note that since python uses 0-based indexing these become 0, and 6 respectively.

Polynomial Fit:

The dispersion curves are fit to a Taylor expansion using [[poly_fit.lsf]]. A linear algebraic approach means that we could choose any order polynomial to solve for; however, we find that 2nd order Taylor expansion is sufficient here.

##This function defines a polynomial fitting of order n_int

function Taylor_exp(x, y, n_int){
#Function that fits x, and y data to a polynomial or order n_int
X = matrix(length(x),n_int+1);
#Set x data in matrix

for(i=0:n_int){ X(1:length(x),i+1) = x^i; }
#Linear algrebra to invert y = X x => x = (X^t X)^{-1} X^t y

a_fit = mult(mult(inv(mult(transpose(X),X)),transpose(X)),y);
return a_fit;
}

The results are then plotted at each temperature.

Step 3

Python:

The python script runs from the Lumerical script editor using a version of python that ships with Lumerical. It is not necessary to have an API license to run python scripts from Lumerical in this way, but a license is required when driving and driving Lumerical solvers.

Tuning Slope:

To find the intercept of the dispersion curves we analytically calculate the root using the quadratic formula at each temperature. Once the intercepts are known the phase-matched wavelength is fit to a line and plotted as a function of temperature.

Conversion Efficiency:

Second-harmonic generation efficiency is given in functional form as.

$$\Gamma \equiv \frac{P_{2}}{P_{1}^{2}}=\eta L^{2}\left[\frac{\sin (\Delta L / 2)}{\Delta L / 2}\right]^{2}$$

Here P1 and P2 is the power in the fundamental and second harmonic respectively. L is the length of the waveguide and ∆ is the phase mismatch \( \Delta =\frac{4\pi}{\lambda}(n_2−n_1)\). The relative conversion efficiency has a sinc function dependence.

The power normalization is given by

$$\eta=\frac{8 \pi^{2}}{\epsilon_{0} c n_{1}^{2} n_{2} \lambda^{2}} \frac{\zeta^{2} d_{\mathrm{eff}}^{2}}{A_{\mathrm{eff}}}$$

Where the effective nonlinear effective area is \( A_{\mathrm{eff}} = \left(A_{1}^{2} A_{2}\right)^{\frac{1}{3}} \) where:

$$A_{i}=\frac{\left(\int_{\text{all}}\left|\vec{E}_{i}\right|^{2} \mathrm{d} x \mathrm{d} y\right)^{3}}{\left|\int_{\chi ^{(2)}} \left|\vec{E}_{i}\right|^{2} \vec{E}_{i} \mathrm{d} x \mathrm{d} y\right|^{2}}, \quad(i=1,2)$$

The spatial mode overlap is determined by integrating the fields over the \(\chi^{(2)} \) region meaning the LNO material

$$\zeta=\frac{\int_{\chi^{(2)}}\left(E_{1 y}^{*}\right)^{2} E_{2 z} \mathrm{d} y \mathrm{d} z}{\left|\int_{\chi^{(2)} } \left|\vec{E}_{1}\right|^{2} \vec{E}_{1} \mathrm{d} x \mathrm{d} y\right|^{\frac{2}{3}} \left| \int_{\chi^{(2)} } \left|\vec{E}_{2}\right|^{2} \vec{E}_{2} \mathrm{d} y \mathrm{d} y\right|^{\frac{1}{3}}}$$

We take the spatial mode overlap, and nonlinear effective indices directly from [1]. Conversion efficiency is calculated in Python and plotted below.

Updating the model with your parameters

Instructions for updating the model based on your device parameters

In this example we used Z-cut Lithium Niobate on Insulator (LNOI) such that the c vector was out of the wafer plane. The parameter ‘cut’ in model analysis can be changed to define a wafer with a different orientation.

Change the waveguide geometry for your process dependent dimensions. Effective indices are dependent on the geometry so the phase matching curve will vary with the actual fabrication parameters. The LNOI structure group has been parameterized to specify width, and sidewall angle,

Ignoring temperature, you could look at frequency mixing process in FDTD using a Chi 2 material. Create a sampled data material using the text file written by model setup script.

Taking the model further

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

Develop a phase-matching device by designing the auxiliary temperature modulator in HEAT. A temperature-dependent index perturbation simulation object can be used to model actual device performance.

Extracting the group velocity, from the dispersion curves - Dispersion Analysis - and export them to INTERCONNECT to see how the waveguide will work other components.

Perform the overlap integrals outlined in conversion efficiency for your own device.

Additional resources

Additional documentation, examples and training material

Related publications

1. R. Luo, Y. He, H. Liang, M. Li, and Q. Lin, "Highly tunable efficient second-harmonic generation in a lithium niobate nanophotonic waveguide," Optica 5, 1006 (2018). https://doi.org/10.1364/OPTICA.5.001006
2. D. Zelmon, D. Small, and D. Jundt, "Infrared corrected Sellmeier coefficients for congruently grown lithium niobate and 5 mol. % magnesium oxide–doped lithium niobate," J. Opt. Soc. Am. B 14 , 3319-3322 (1997) https://doi.org/10.1364/JOSAB.14.003319
3. L. Moretti, M. Iodice, F. G. Della Corte, and I. Rendina, “Temperature dependence of the thermo-optic coefficient of lithium niobate, from 300 to 515 K in the visible and infrared regions,” J. Appl. Phys. 98(3), (2005) https://doi.org/10.1063/1.1988987

See also

• Waveguide FDE
• Waveguide FEEM
• Thermally Tuned Waveguide FDE
• Thermally tuned waveguide (FEEM)
• Lithium Niobate Nonlinear Thermal Waveguide (FDE)

Related Ansys Innovation Courses

Appendix

Additional background information and theory

Material models

Given the precise phase matching requirements in harmonic generation, accounting for the thermal, and wavelength dependence becomes important. To build this model we need to understand the properties of Lithium Niobates crystal structure which is composed of trigonal unit cells that have 3m point symmetry -meaning threefold rotation and a plane of mirror symmetry. This underlying crystal symmetry is the origin of linear anisotropy and second-order non-linearity. LNO is a negative uniaxial crystal meaning that light polarized along the crystal axis c, which lacks inversion symmetry, and experiences an extraordinary index of refraction less than the ordinary index. Light with zero E-field component projected along the c axis experienced the normal index. The offset between the ordinary and extraordinary refractive index is dn ~ 0.1, which is 10 times that of crystalline SiO2, and 100 times larger than water ice. Linear anisotropy is associated with birefringence where different polarizations produce a double refraction due to different indices they experience at the interface.

Lithium Niobate is a synthetic inorganic crystal since does not appear naturally, but it can be grown in a similar fashion to more common semiconductors. Its adoption in telecommunications is due to the broad transmission window and the high second-order nonlinearity. The transparency window and index are consistent with glass and so it appropriate to apply a Sellmeier model; however, this needs to be done for the ordinary and extraordinary axes separately. Like silica glass the actual dispersion will depend on the temperature and dopants present; furthermore, the optical properties of single-crystal LNO will depend on the stoichiometry due to the presence of crystal defects and vacancy sites.

To describe the material dispersion of LNO we have employed the three oscillator Sellemier model based on experimental coefficients fit by Zelmon et al [2] for congruent lithium niobite. In its most general form, the Sellmeier model describes n resonant oscillators with strength \( B_i \) and wavelength \( \sqrt(C_i)\}.

$$n^{2}(\lambda)=1+\sum_{i} \frac{B_{i} \lambda^{2}}{\lambda^{2}-C_{i}}$$

The values used in our models are shown below.

Table 1: Sellmeier coefficients

Temperature-dependent material parameters are taken from Morretti et al [3]. In this paper, the shift of the ordinary, and extraordinary indices with respect to temperature is measured at 632 and 1523nm. We can think of these as lines a long 632nm and 1555nm in 3D space (\( dn / dT, \lambda, T \)), with a linear dependence on temperature.

At optical frequencies these were found to be.

$$\begin{aligned}
&\frac{d n_{o}}{d T} (\lambda = 632nm)=-1.7+6.9 \times 10^{-3} T\left(10^{-5} \mathrm{K}^{-1}\right)\\
&\frac{d n_{e}}{d T} (\lambda = 632nm)=-2.6+22.4 \times 10^{-3} T\left(10^{-5} \mathrm{K}^{-1}\right)
\end{aligned}$$

And at telecom they were found.

$$\begin{aligned}
&\frac{d n_{o}}{d T} (\lambda = 1523nm)=-0.9+3.0 \times 10^{-3} T\left(10^{-5} \mathrm{K}^{-1}\right)\\
&\frac{d n_{e}}{d T}(\lambda = 1523nm)=-2.6+19.8 \times 10^{-3} T\left(10^{-5} \mathrm{K}^{-1}\right)
\end{aligned}$$

These are skew lines in 3D, meaning we cannot fit a 2D plane solution for dn/dT; however, we can approximate the wavelength dependence by fitting each coefficient separately. Since the equations from Moretti given above are of the form

$$\frac{d n_{*}}{d T}= a_∗+b_∗T $$

We postulate the functions satisfy.

$$ \frac{d n_{*}}{d T}(\lambda,T)=a_∗(\lambda)+b_∗(\lambda)T $$

Given this information, we fit a and b to linear functions of lambda using the known coefficients at 632nm, and 1523nm.

Table 2: Temperature coefficients

Thus we have two polynomials, ordinary and extraordinary, of the form

$$p(\lambda, T)=c_{00}+c_{01} T+c_{10} \lambda+c_{11} \lambda \mathrm{T}$$

Integrating dn/dT

$$n_{*}(\lambda, T)=a_{*}(\lambda)(T)+b_{*}(\lambda)\left(T^{2}\right)+C(\operatorname{\lambda})$$

Ensuring the dispersion is equivalent to that defined in [2] at \(T_0\) = 293.3K = 23C we determine the constants of integration with S the Sellmeier equation.

$$n_{*}(\lambda, T)=a_{*}(\lambda)\left(T-T_{0}\right)+b_{*}(\lambda)\left(T^{2}-T_{0}^{2}\right)+S(\lambda)$$

Phase Matching

Transparent materials follow a normal dispersion such that the material indices decrease with increasing wavelength. In order to achieve decent conversion efficiency, we need the pump and the second harmonic light to be in phase. In nonlinear optics the coherence length is used to refer to the distance over which the two different signals will move out of phase by a factor of \(\pi\).

$$ L_c = \frac{\pi}{\Delta k} = \frac{1}{2} \frac{\lambda}{n_2-n} $$

One needs to be careful as this has no relation with the usual meaning of coherence length, but in this article, we are strictly referring to the nonlinear coherence length. Beyond the coherence length, a relative (\{pi\) phase shift means that frequency conversion becomes ineffective, and making a device longer than the coherence length does not make sense. At 1550nm we have a coherence length of \ L_c \) ~ 16um, and 19um for ordinary and extraordinarily polarized beams. SHG efficiency is not high enough to produce a usable signal over this distance.

In bulk LNO it is possible to achieve phase matching by rotating the angle of the beam, known as angle tuning. Reducing the effective index is accomplished by reducing the projection of the E field along the ordinary axis and increasing the projection along the extraordinary axis. By choosing the direction of propagation, and polarization we can tune the index of light to be anywhere between \(n_e\) and \(n_o\). This is realized in parametric oscillators/amplifiers and in tilted pulse THz generation. In the image below it is shown that by aligning the polarization of the telecom slightly off from \(n_o\) and its second harmonic slightly off \(n_e\) there is a green region which would allow for phase matching. Achieving momentum conservation in this way is highly sensitive to alignment, and this simple argument does not consider nonlinear efficiency and beam walk-off.

Another phase-matching scheme of practical importance for Lithium Niobate is accomplished by periodically changing the polarity of the ferroelectric domains. Analogous to a permanent magnet, or ferromagnetic material, Lithium Niobate is ferroelectric meaning it has built-in polarization domains. By applying a strong DC E-field along the c-axis it is possible to flip the polarity of these domains. Fabricating electrodes on the same length scale as the \(L_c\) one can periodically flip the polarization and achieve quasi-phase-matching. These types of devices are commonly referred to as periodically poled lithium niobite PPLN and are employed for efficient frequency conversion devices in telecommunications; however, the interaction length required for such devices is typically on the order of 10cm.

Effective nonlinear coefficients and Overlap calculations

The optical properties of most materials can be expressed as an expansion of the polarization in terms of the electric field and their coefficients given as electric susceptibilities which are tensors of the same order as the exponential.

$$P(\omega)=\varepsilon_{0} \chi_{e}^{(1)}(\omega) E(\omega)+\varepsilon_{0} \chi_{e}^{(2)}(\omega) E^{2}(\omega)+\varepsilon_{0} \chi_{e}^{(3)}(\omega) E^{3}(\omega)+\cdots $$

In most cases the first-order term dominates, and we have the familiar linear relationship between the material polarization and the electric field. Physically the polarization is the macroscopic response of the constituent dipoles to an applied external electric field. With no underlying preferential orientation, the dipoles will oscillate in the plane of the electric field. The linear response is due to the electron orbitals bound by a harmonic spherically symmetric restoring force, i.e. Coulomb potential. Coupling between the optical fields and the electrons in the harmonic potential produces a sine wave at the same frequency as the driving field, but with the phase velocity altered by the strength of the linear susceptibility. Although any potential can be approximated as parabolic for small amplitude oscillations, when the field amplitudes are high enough the anharmonicity becomes important. In LNO the strength of the second-order nonlinear coefficients is due to the difference in electronegativity of the Lithium and Niobium ions. Once the harmonic nature is broken the oscillatory behavior becomes non-sinusoidal, but the output can be broken into a superposition of sinusoids via Fourier analysis. These sinusoids will have angular frequencies \( n*\omega\) where \( n = 0, 1 , 2, 3,\cdots \).

Furthermore, when the dipoles experience different restoring forces there can easily be mixing between field components. This is simplified when the crystal axes are aligned with the principle coordinates. Given that amorphous solids like glass are isotropic meaning they have molecules aligned randomly the higher-order contributions cancel out, and nonlinear phenomenon glass arises due to the inclusion of organic or metallic nanoparticles that mediate a nonlinear interaction differently.