Second harmonic generation (SHG) in a Lithium Niobite  LiNbO3 (LNO) nanophotonic waveguide is studied using temperature modulation 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 2nd harmonic mode using FDE. 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 and external electric circuitry could support this essential functionality. In the following example, we demonstrate a 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 [[LNO_waveguide.lms]].
 Run the script file [[waveguide_modes.lsf]]
 (Optional) edit model to plot material and/or save to text file
In order to perform the sweep in the next step, we need to find the reference modes. To do this the script [[waveguide_modes.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.
Step 2: Temperature, Wavelength Sweep
 [[step]] (optional) Change your resources to solve sweep in parallel
 Run the sweep Telecom_temp_sweep
 Change the mode tracking to TM2 in model analysis
 Run the sweep NIR_temp_temp_sweep .
 Run [[poly_fit.lsf]].
Wavelength and temperature range/resolution are set in the sweep objects Telecom_temp_sweep, and NIR_temp_temp_sweep. We have chosen 4 wavelength points and 3 temperatures [0,100C] for each spectral band. To ensure that the analysis is valid the wavelengths should be complementary meaning \( \lambda_{NIR} = \frac{1}{2}*\lambda_{Tel}\); furthermore, the same temperature points should be used. Run the telecom sweep extracting the mode profiles and effective indices. Next, update the ModeName parameter in the model analysis script to TM2 and run the NIR_temp_sweep sweep. Once both sweeps have completed run the [[poly_fit.lsf]]. This file will fit a 2nd order polynomial to the dispersion curves at each temperature, plot, and save the results 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, including it's evanescent tails. Ensure the field has decayed to at least 1e4 near the boundaries. This is easy to visualize by plotting the fields on a log scale.
Boundary conditions :
Use metal boundaries. If the fields are negligible at the boundary (as recommended in the previous comment), then the choice of boundary is not important. Metal boundaries are fastest and minimize the simulation memory. PML boundaries a required for sensitive leakage and bending loss calculations, but this done with much greater computational cost.
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.
Fig 1: LiNbO\(_3\) Material indices
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 Fig. 2 below. This will prove to be enough to produce tuning across most of the spectral bands.
Fig 2: LiNbO\(_3\) Index temperature dependence
Modes
The effective index is highest for the fundamental mode and decreases with higherorder, due to less confinement in the core region. By looking for a higher mode at the second harmonic, we should be able to identify a mode with a similar effective index to the fundamental mode at telecom. The geometric contributions will have less importance at the second harmonic (SH) 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. LNO is a negative uniaxial crystal, so the 1.55um mode will be polarized along the ordinary axis(blue) and the SH mode 0.775um should primarily be polarized along the extraordinary axis(yellow) in Fig 1.
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_modes.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 Fig 4. These reference modes are saved to the dcard for use in the next step.
Fig 3: Fundamental TE0 mode at 1.55\(\mu\)m
Fig 4: Secondorder TM2 mode at 0.75\(\mu\)m
In the visualizer, one can look at the polarization, to confirm for themselves the above arguments.
Step 2:
Reference modes :
The reference modes we calculated in step 1 are saved to the dcard so that we can choose the mode with the best overlap while sweeping temperature and wavelength. You can adjust the model analysis script parameter ModeName to change the tracked mode.
Fig 5: Model analysis window
Sweep:
The total sim number \( N = n_T \times n_{\lambda} \times 2\) = 3*4*2 = 24 these are defined as two nested sweeps that will need to be run separately.
Fig 6: Temperature sweep
Fig 7: Wavelength sweep
Manage your resources to improve throughput. Since all the simulations in this workflow can be performed concurrently you can update your resource capacity. By running each simulation on a subset of CPU cores you will see are marginal slow down in each simulation, but the total time to run all simulations will be reduced. The following configuration will run 6 simulations simultaneously on the local computer using 1 process and 2 threads each.
Fig 8: Resource configuration window
For best results, you should ensure that the # processes *# threads = #logical_cores on your machine. This is set up to work on a machine with 12 logical cores. If in doubt run a configuration test on the current resource configuration using Run tests button.
In previous versions of Lumerical the capacity option was not available, but you could still perform multithreading by duplicating the defined resource. To achieve the same setup in Fig 7 you would duplicate 5 times to define 6 local resources.
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.
Fig 9: Dispersion curve fitting
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.
Fig 10: Temperaturedependent phasematching slope
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.
Fig 11: Conversion efficiency for different temperatures
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. In our coordinate system this corresponds to 'Y' cut. The parameter ‘cut’ in model analysis can be changed to define a wafer with a different orientation. For crystal axes not aligned with the principle, axes use a grid attribute to change the orientation of the underlying crystal lattice.
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
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.
Fig 12: Phase matching window
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 \).
Fig 13: Anharmonic potential well
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.