A generalized formalism is developed to model second-order nonlinear processes in finite-difference time-domain (FDTD) simulations. The method is capable of modeling frequency-conversion from all 18 elements of the second-order nonlinear tensor, where dispersion of the tensor elements is included at both the pump and generated frequencies. The model is validated by considering frequency-conversion in a LiNbO3 crystal, which has highly dispersive second-order nonlinear susceptibilities near the phonon resonances. The developed nonlinear formalism is able to model any arbitrary excitation polarization state and can be applied to investigate second-order nonlinear processes in type I or type II phase-matching. This generalized second-order nonlinear formalism represents an advancement for the FDTD computational technique and can provide more realistic modeling of second-order nonlinear interactions in nanoscale devices and waveguiding structures.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
The finite-difference time-domain (FDTD) technique has evolved to become one of the most powerful numerical methods in modeling electromagnetic phenomena and optical devices. In part, this is due to the simplicity of the numerical algorithm in calculating the fields’ evolution in real time. With recent advances in integrated nonlinear photonics, the FDTD method has proven to be ideal in modeling on-chip devices based on second or third order optical nonlinearities. However, such FDTD approaches have been restricted to a basic representation of the nonlinear material properties, such that intricate and complex nonlinear interactions have not been fully explored by these models. In particular, this is very important when dealing with a second-order nonlinear optical material having highly-dispersive nonlinearities and multiple non-vanishing nonlinear tensor components. Naturally, the next major advancement for FDTD modeling of second-order nonlinear effects [e.g. second harmonic generation, sum frequency generation, difference frequency generation (DFG), and optical rectification] is to incorporate realistic models for the nonlinear material.
Second-order nonlinear optical materials, especially those that are rich in optical resonances, exhibit high nonlinear dispersion across a wide spectral range. Second-order nonlinear susceptibility, , dispersion is prominent near the phonon resonance of a material , where crystals such as LiNbO3, LiTaO3, GaAs, GaP, InP, and CdTe exhibit enhancement of their tensor coefficients [1–5]. Another situation where dispersion occurs is near the band edge of a material. Here, can vary by as much as ~100 pm/V for point group materials (e.g. GaAs, GaP, ZnTe, ZnSe, ZnS, InAs, and InSb) [6–8]. Likewise, dispersion is crucial for organic molecules  and 2D materials , as their coefficients differ by >100% across the near-infrared and visible spectral ranges, respectively. Since a high-degree of complexity arises when including all the dispersive elements into FDTD formalisms, current FDTD methods implement a maximum of three dispersive tensor elements. Ignoring the remaining 15 dispersive elements lead to critical restrictions, as it precludes realistic modeling of second-order nonlinear effects for numerous material classes. In 3m point group materials, frequency-conversion can occur via eight nonlinear elements (i.e. , , , , , , , and ) . Even more involved are organic crystals, where 10 of the 18 nonlinear elements are non-zero for 2-methyl-4-nitroaniline .
Here, we present a generalized FDTD method for modeling dispersive second-order nonlinear effects. The model is based on Miller’s rule, which allows dispersion to be included for all 18 elements in the second order-nonlinear tensor. To validate the developed formalism, frequency-conversion is investigated in a LiNbO3 crystal, which has highly dispersive elements near its phonon modes. The developed method is able to model any arbitrary excitation polarization state and can be applied to investigate nonlinear processes in type I or type II phase-matched configurations. This generalized second-order nonlinear formalism represents a major advancement for the FDTD computation technique and is envisioned to have a strong impact in integrated optics.
2. Derivation of the second-order nonlinear current density
Light interaction with a non-centrosymmetric material can be described by the current density, , where and are the linear (i.e. first-order) and the second-order current densities, respectively. Although and its implementation into an FDTD formalism has been well-established , a generalized FDTD formalism needs to be developed for nonlinear light interactions involving . When a non-centrosymmetric material is excited using optical electric fields having the frequencies of ω and Ω-ω, a nonlinear optical current density, , is induced at the frequency Ω. The second-order nonlinear susceptibility tensor, utilizing contracted notation for the tensor elements, can be written as,1],Eq. (4) for , Eq. (3) becomes,Eq. (5) is the convolution operation, such that it can be recasted in a more compact form,Eq. (8), the complete set of second-order nonlinear polarizations along the x, y, and z directions are ,Equations (12-14) allow the full, dispersive second-order nonlinear susceptibility tensor to be implemented in the FDTD algorithm. Nonetheless, an accurate representation of the linear susceptibility is required. Since most physical processes permit the linear susceptibility (i.e. ) to be represented by a frequency-independent term and a summation of Lorentzian oscillators, the linear susceptibilities along the x, y, and z directions are,
3. Discretization of the second-order nonlinear current densityEq. (17), it is first transformed to the time-domain,Eq. (19) is then discretized for the time iteration of n to obtain,Eq. (18), the terms are first rearranged to,Eq. (21) can be transformed to its time-domain form,Eq. (16) to the time-domain and discretizing it for the time iteration of n + 1,Eq. (20) and (23), which is achieved by using the results from the current time iteration at the next time iteration. Next, using Eq. (15), in Eq. (7) can be recasted as,Equation (29) is solved by first converting it to the time-domain,Eq. (24). Equation (30) is solved by first rearranged it into the form,Eq. (33), we obtain the time-domain equation,Eq. (34) is discretized for the time iteration of n,Eq. (35) depends on the solution in Eq. (27). Now that and are obtained in their discretized form, a discretized solution can be attained for by converting Eq. (28) to the time-domain and discretizing it for the time iteration of n + 1,Eqs. (32) and (35), respectively. Finally, the discretization of [Eqs. (9-11)] for the time iteration of n + 1 is straightforward and results in,Eqs. (12-14)] for the time iteration of n + 1/2 provides,
4. Frequency-conversion in a LiNbO3 crystal
The derived FDTD method is evaluated by simulating the representative effects of optical rectification and DFG. These nonlinear effects are investigated using a LiNbO3 crystal layer having a length l, which is illustrated in Fig. 1. LiNbO3 serves as a prime material to evaluate the generalized second-order nonlinear method, since it has strong nonlinear dispersion in the terahertz (THz) frequency regime, contains a nonlinear tensor with numerous non-vanishing elements, and has sufficient nonlinear experimental data available. The LiNbO3 layer is excited using an electric field pulse having a central-wavelength of 1550 nm, a duration of 80 fs, and a polarization angle of θ with respect to the c-axis of the crystal (see Fig. 1). Here, it is assumed that an index-matched layer and free space are positioned at the input and output faces of the LiNbO3 layer, respectively. The c-axis of the crystal is oriented along the z-axis, the  crystal direction is aligned with the x-axis, and the y-axis is defined as the direction of propagation. Since a bulk crystal is being simulated, periodic boundary conditions are implemented for the boundaries normal to the x and z-directions, whereas perfectly matched layers are used for boundaries normal to the propagation direction (i.e. y-direction). A mesh size of 40 nm is used, which is sufficiently small for both the excitation and generated frequency components.
It is critical to determine the frequency dependence of the LiNbO3 refractive indices and extinction coefficients, along with its second-order nonlinear susceptibilities. The extraordinary, ne, and ordinary, no, refractive indices are shown in Fig. 2(a), 2(c), and 2(e) for the frequency ranges of interest. Furthermore, the extraordinary, кe, and ordinary, кo, extinction coefficients are shown in Fig. 2(b) and 2(d), respectively. The experimental data is obtained from references [15,16]. and the curve fits are achieved by using a superposition of Lorentzian oscillators. Clearly, the experimental data matches the fitted curves very well. ne and кe show a phonon resonance at 7.6 THz, corresponding to the lowest-frequency A1 mode of LiNbO3, whereas three E mode phonon resonances are observed in no and кo, which are located at 4.6, 7.9, and 9.7 THz . Since nonlinear dispersion is directly related to linear dispersion via Miller’s rule, the elements of the LiNbO3 crystal cannot be taken as constants. The LiNbO3 crystal belongs to the 3m point group symmetry class, such that its second-order nonlinear susceptibly tensor is written as ,Eq. (3), (9), (11), (49), and the fact that propagation is along the y-direction, it can be determined that the second-order nonlinear polarizations occurring along the x and z directions are,Figure 3(a) illustrates the highly-dispersive element, which is calculated using Miller’s rule [Eq. (4)], a proportionality constant of δ = 1.3 pm/V, and the experimental data in [17,18]. Similarly, the element [see Fig. 3(b)] is obtained by using a proportionality constant of δ = 0.21 pm/V and the experimental data in , whereas [see Fig. 3(c)] is described by δ = 0.71 pm/V and the experimental data from . Importantly, in comparison to the low-frequency nonlinear susceptibility values, all three of the nonlinear elements exhibit an enhancement of >7 times at their lowest-frequency phonon resonance.
To illustrate the implication of a dispersive second-order nonlinear susceptibility on the nonlinear frequency-conversion process, simulations are performed using l = 150 µm and θ = 0°, such that the 1550 nm excitation pulse is polarized along the c-axis. Here, the only contributing nonlinear element is . The simulations are conducted using the fully-dispersive element [see Fig. 3(a)], as well as the frequency-independent = 348 pm/V. Figure 4(a) depicts the spectral power generated along the z-direction, ξz, where the phase-mismatching effects appear as dips (e.g. f = 0.8 THz, 1.5 THz, etc.) in the power spectrum, in agreement with those predicted by the theoretical formula,Figure 4(b) shows the relative spectral power, which is defined as ξz obtained using the dispersive divided by ξz obtained using the frequency independent = 348 pm/V. While the discrepancy in the relative spectral power is moderate at frequencies <6 THz, it is more than 60 times higher at frequencies near 7.6 THz. Therefore, a dispersive nonlinear element is critical to accurately model frequency-conversion near the phonon resonance of a crystal.
The nonlinear algorithm can be used to model dispersive type II phase-matching for l = 50-150 µm and θ = 45°. Herein, the THz radiation generation occurs along the x-direction through the process of DFG. In DFG, phase-matching can occur between the pump frequency, fp, the signal frequency, fs, and either the forward, , or backward, , propagating idler frequencies. The forward propagation DFG coherence length, , is,Equation (53)], phase-matching occurs at = 3 THz, as shown in Fig. 5(a). This agrees with the FDTD model, where the x-component of the spectral power, ξx, exhibits THz radiation generation at = 3 THz, as seen in Fig. 5(b). Similarly, from Eq. (54) is displayed in Fig. 5(c), which shows that phase-matching occurs at the idler frequency of = 1.8 THz. By performing FDTD simulations and recording ξx near the input facet of the LiNbO3 layer, we confirm that = 1.8 THz [Fig. 5(d)]. Clearly, the versatility of this nonlinear FDTD approach arises from its ability in modeling excitation electric fields at arbitrary polarization states, along with the dispersive second-order nonlinear susceptibilities.
Next, by using l = 150 µm and θ = 0°-90°, we examine the two simultaneous nonlinear frequency-conversion processes of DFG and optical rectification in the THz regime. These two processes depend on the polarization state of the excitation electric field, and the generated THz radiation is polarized along different crystal directions. Here, the broadband THz radiation generated via optical rectification is polarized along the z-axis, whereas the narrowband THz radiation generated via DFG is polarized along the x-axis. Figure 6(a) depicts ξz, where the broadband THz radiation generation is highest when the excitation pulse is oriented along the c-axis (i.e. θ = 0°), since is maximum and the only contribution is from the term [see Eq. (51)]. THz radiation generation decreases as θ increases, since the contribution from is reduced; nonetheless, the weaker coefficient now contributes to . At θ = 90°, THz radiation generation is lowest, since is only influenced by . Figure 6(b) depicts ξx, where the THz radiation generation vanishes at both θ = 0° and 90°, since = 0 [see Eq. (50)]. The power spectra are similar in magnitude at the intermediate angels of θ = 22.5° and 67.5°, whereas THz radiation generation is highest when θ = 45°. The presented generalized FDTD method offers more accurate modeling of dispersive nonlinear frequency-conversion processes that strongly depend on the excitation polarization state and the generated electric field polarization.
A generalized nonlinear FDTD formalism is developed that considers all 18 elements of the second-order nonlinear tensor and incorporates dispersion of the nonlinear elements. The versatility of the method was demonstrated by modeling optical rectification and DFG in a LiNbO3 crystal, where the tensor elements are highly dispersive at the phonon resonances and the enhancement of the nonlinear interaction is more pronounced. When implementing the generalized nonlinear formalism, THz radiation generation near the phonon resonance is found to be 60 times higher than that obtained using frequency-independent nonlinear methods. Furthermore, the nonlinear formalism permits frequency-conversion for an arbitrary electric field polarization state. Although the implantation was focused on nonlinear frequency-conversion in the LiNbO3 crystal, the nonlinear algorithm can be applied to model more complex structures, such as nano-dimensional waveguides that support a longitudinal electric field, plasmonic structures that enhance frequency-conversion , or the generation of radiation in bulk materials that are excited by a focused laser beam.
Natural Sciences and Engineering Research Council of Canada.
This work was performed under NSERC’s Engage grant program in collaboration with Optiwave Systems Inc. The authors are grateful for the support provided by Optiwave Systems Inc., especially Scott Newman, Ahmad Atieh, and Kevin Chu, and for the access to Optiwave's OptiFDTD simulation software platform.
1. D. H. Auston and M. C. Nuss, “Electrooptical generation and detection of femtosecond electrical transients,” IEEE J. Quantum Electron. 24(2), 184–197 (1988). [CrossRef]
2. Y. J. Ding, “Efficient generation of high-frequency terahertz waves from highly lossy second-order nonlinear medium at polariton resonance under transverse-pumping geometry,” Opt. Lett. 35(2), 262–264 (2010). [CrossRef] [PubMed]
3. M. Cherchi, A. Taormina, A. C. Busacca, R. L. Oliveri, S. Bivona, A. C. Cino, S. Stivala, S. R. Sanseverino, and C. Leone, “Exploiting the optical quadratic nonlinearity of zinc-blende semiconductors for guided-wave terahertz generation: a material comparison,” IEEE J. Quantum Electron. 46(3), 368–376 (2010). [CrossRef]
5. D. E. Thompson and P. D. Coleman, “Step-tunable far infrared radiation by phase matched mixing in planar-dielectric waveguides,” IEEE Trans. Microw. Theory Tech. 22(12), 995–1000 (1974). [CrossRef]
6. J. Seres, “Dispersion of second-order nonlinear optical coefficient,” Appl. Phys. B 73(7), 705–709 (2001). [CrossRef]
7. H. P. Wagner, M. Kühnelt, W. Langbein, and J. M. Hvam, “Dispersion of the second-order nonlinear susceptibility in ZnTe, ZnSe, and ZnS,” Phys. Rev. B Condens. Matter Mater. Phys. 58(16), 10494–10501 (1998). [CrossRef]
8. M. I. Bell, “Frequency Dependence of Miller’s Rule for Nonlinear Susceptibilities,” Phys. Rev. B 6(2), 516–521 (1972). [CrossRef]
9. R. Morita, T. Kondo, Y. Kaneda, A. Sugihashi, N. Ogasawara, S. Umegaki, and R. Ito, “Dispersion of second-order nonlinear optical coefficient d11 of 2-Methyl-4-Nitroaniline (MNA),” Jpn. J. Appl. Phys. 27(6), L1131–L1133 (1988). [CrossRef]
10. M. Mokim, A. Card, B. Sah, and F. Ganikhanov, “Dispersion of the resonant second order nonlinearity in 2D semiconductors probed by femtosecond continuum pulses,” AIP Adv. 7(10), 105121 (2017). [CrossRef]
11. G. D. Boyd, R. C. Miller, K. Nassau, W. L. Bond, and A. Savage, “LiNbO3: An efficient phase matchable nonlinear optical material,” Appl. Phys. Lett. 5(11), 234–236 (1964). [CrossRef]
12. R. Vallée, P. Damman, M. Dosiére, E. Toussaere, and J. Zyss, “Nonlinear optical properties and crystalline orientation of 2-Methyl-4-nitroaniline layers grown on nanostructured poly(tetrafluoroethylene) substrates,” J. Am. Chem. Soc. 122(28), 6701–6709 (2000). [CrossRef]
13. A. Taflove, Computational electrodynamics: the finite-difference time-domain method (Artech House Publishers, 1995).
14. M. C. Gupta and J. Ballato, The Handbook of Photonics (CRC, 2007).
15. S. Kojima, K. Kanehara, T. Hoshina, and T. Tsurumi, “Optical phonons and polariton dispersions of congruent LiNbO3 studied by far-infrared spectroscopic ellipsometry and Raman scattering,” Jpn. J. Appl. Phys. 55(10S), 10TC02 (2016). [CrossRef]
16. A. S. Barker Jr. and R. Loudon, “Dielectric Properties and Optical Phonons in LiNbO3,” Phys. Rev. 158(2), 433–445 (1967). [CrossRef]
17. G. D. Boyd and M. A. Pollack, “Microwave nonlinearities in anisotropic dielectrics and their relation to optical and electro-optical nonlinearities,” Phys. Rev. B 7(12), 5345–5359 (1973). [CrossRef]
18. G. D. Boyd, T. J. Bridges, M. A. Pollack, and E. H. Turner, “Microwave nonlinear susceptibilities due to electronic and ionic anharmonicities in acentric crystals,” Phys. Rev. Lett. 26(7), 387–390 (1971). [CrossRef]
19. N. Courjal, M.-P. Bernal, A. Caspar, G. Ulliac, F. Bassignot, L. Gauthier-Manuel, and M. Suarez, “Lithium niobate optical waveguides and microwaveguides,” Emerging Waveguide Technology (IntechOpen, 2018).