Secondharmonic 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 phasematching condition. The LNO indices have distinct temperature dependencies allowing the phase matching to be thermally tuned.
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 IIIV, 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 temperaturedependent 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
 [[step]]Open [[linbo03_wg.ldevs]].
 Run the script file [[waveguide_modesFEEM.lsf ]]
 (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 secondharmonic TM2 mode at the center wavelengths and the midpoint operating temperature. The mode Efields are then visualized for the user to inspect.
Step 2: Temperature, Wavelength Sweep
 [[step]] (optional) Change your resources to solve sweep in parallel
 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
 [[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 wavelengthdependent 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 1e4 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 finiteelement 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 setup script. These functions have the following form.
$$n_{*}(\lambda, T)=a_{*}(\lambda)\left(TT_{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. Thermooptic 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.040.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, 0100C, 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 copropagate, 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 0based 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 phasematched wavelength is fit to a line and plotted as a function of temperature.
Conversion Efficiency:
Secondharmonic 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 Zcut 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 phasematching device by designing the auxiliary temperature modulator in HEAT. A temperaturedependent 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
 R. Luo, Y. He, H. Liang, M. Li, and Q. Lin, "Highly tunable efficient secondharmonic generation in a lithium niobate nanophotonic waveguide," Optica 5, 1006 (2018). https://doi.org/10.1364/OPTICA.5.001006
 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 , 33193322 (1997) https://doi.org/10.1364/JOSAB.14.003319
 L. Moretti, M. Iodice, F. G. Della Corte, and I. Rendina, “Temperature dependence of the thermooptic 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 secondorder nonlinearity. 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 Efield 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 secondorder 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 singlecrystal 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.
Congruent LNO 
B1 
C1 
B2 
C2 
B3 
C3 
\( n_o \) 
2.6734 
0.01764 
1.229 
0.05914 
12.614 
474.6 
\( n_e \) 
2.9804 
0.02047 
0.5981 
0.0666 
8.9543 
416.08 
Table 1: Sellmeier coefficients
Temperaturedependent 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.
\( a_∗\) [10e5] 
\(b_∗\) [\(K^{1}\) 10e−8] 


\( m_a\) 
\( \alpha \) 
\( m_b\) 
\( \beta \) 
\( n_o \) 
0.89 
2.267 
4.377 
9.666 
\( n_e \) 
0 .0 
2.6 
2.918 
24.244 
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(TT_{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_2n} $$
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 walkoff.
Another phasematching 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 builtin polarization domains. By applying a strong DC Efield along the caxis 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 quasiphasematching. 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 firstorder 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 secondorder 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 nonsinusoidal, 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 higherorder contributions cancel out, and nonlinear phenomenon glass arises due to the inclusion of organic or metallic nanoparticles that mediate a nonlinear interaction differently.