A modeling of dispersive tensorial second-order nonlinear effects for the finite-difference time-domain method

B. N. Carnio; A. Y. Elezzabi

doi:10.1364/OE.27.023432

1. Introduction

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, $χ^{(2)}$ , dispersion is prominent near the phonon resonance of a material [1], where crystals such as LiNbO₃, LiTaO₃, GaAs, GaP, InP, and CdTe exhibit enhancement of their $χ^{(2)}$ tensor coefficients [1–5]. Another situation where $χ^{(2)}$ dispersion occurs is near the band edge of a material. Here, $χ^{(2)}$ can vary by as much as ~100 pm/V for $\bar{4} 3 m$ point group materials (e.g. GaAs, GaP, ZnTe, ZnSe, ZnS, InAs, and InSb) [6–8]. Likewise, $χ^{(2)}$ dispersion is crucial for organic molecules [9] and 2D materials [10], as their $χ^{(2)}$ 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 $χ^{(2)}$ elements into FDTD formalisms, current FDTD methods implement a maximum of three dispersive tensor elements. Ignoring the remaining 15 dispersive $χ^{(2)}$ 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. $χ_{15}^{(2)}$ , $χ_{16}^{(2)}$ , $χ_{21}^{(2)}$ , $χ_{22}^{(2)}$ , $χ_{24}^{(2)}$ , $χ_{31}^{(2)}$ , $χ_{32}^{(2)}$ , and $χ_{33}^{(2)}$ ) [11]. Even more involved are organic crystals, where 10 of the 18 nonlinear elements are non-zero for 2-methyl-4-nitroaniline [12].

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 LiNbO₃ crystal, which has highly dispersive $χ^{(2)}$ elements near its phonon modes. The developed method is able to model any arbitrary excitation polarization state and can be applied to investigate nonlinear $χ^{(2)}$ 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, $J = J^{(1)} + J^{(2)}$ , where $J^{(1)}$ and $J^{(2)}$ are the linear (i.e. first-order) and the second-order current densities, respectively. Although $J^{(1)}$ and its implementation into an FDTD formalism has been well-established [13], a generalized FDTD formalism needs to be developed for nonlinear light interactions involving $J^{(2)}$ . When a non-centrosymmetric material is excited using optical electric fields having the frequencies of ω and Ω-ω, a nonlinear optical current density, $J^{(2)}$ , is induced at the frequency Ω. The second-order nonlinear susceptibility tensor, utilizing contracted notation for the tensor elements, can be written as,

\bar{\bar{χ^{(2)}}} (Ω, ω, Ω - ω) = [\begin{matrix} \begin{matrix} χ_{11}^{(2)} (Ω, ω, Ω - ω) & χ_{12}^{(2)} (Ω, ω, Ω - ω) & χ_{13}^{(2)} (Ω, ω, Ω - ω) \\ χ_{21}^{(2)} (Ω, ω, Ω - ω) & χ_{22}^{(2)} (Ω, ω, Ω - ω) & χ_{23}^{(2)} (Ω, ω, Ω - ω) \\ χ_{31}^{(2)} (Ω, ω, Ω - ω) & χ_{32}^{(2)} (Ω, ω, Ω - ω) & χ_{33}^{(2)} (Ω, ω, Ω - ω) \end{matrix} & \begin{matrix} χ_{14}^{(2)} (Ω, ω, Ω - ω) & χ_{15}^{(2)} (Ω, ω, Ω - ω) & χ_{16}^{(2)} (Ω, ω, Ω - ω) \\ χ_{24}^{(2)} (Ω, ω, Ω - ω) & χ_{25}^{(2)} (Ω, ω, Ω - ω) & χ_{26}^{(2)} (Ω, ω, Ω - ω) \\ χ_{34}^{(2)} (Ω, ω, Ω - ω) & χ_{35}^{(2)} (Ω, ω, Ω - ω) & χ_{36}^{(2)} (Ω, ω, Ω - ω) \end{matrix} \end{matrix}] .

However, in the following derivations, it is more intuitive to implement the uncontracted tensor element representation,

\bar{\bar{χ^{(2)}}} (Ω, ω, Ω - ω) = [\begin{matrix} \begin{matrix} χ_{x x x}^{(2)} (Ω, ω, Ω - ω) & χ_{x y y}^{(2)} (Ω, ω, Ω - ω) & χ_{x z z}^{(2)} (Ω, ω, Ω - ω) \\ χ_{y x x}^{(2)} (Ω, ω, Ω - ω) & χ_{y y y}^{(2)} (Ω, ω, Ω - ω) & χ_{y z z}^{(2)} (Ω, ω, Ω - ω) \\ χ_{z x x}^{(2)} (Ω, ω, Ω - ω) & χ_{z y y}^{(2)} (Ω, ω, Ω - ω) & χ_{z z z}^{(2)} (Ω, ω, Ω - ω) \end{matrix} & \begin{matrix} χ_{x y z}^{(2)} (Ω, ω, Ω - ω) & χ_{x x z}^{(2)} (Ω, ω, Ω - ω) & χ_{x x y}^{(2)} (Ω, ω, Ω - ω) \\ χ_{y y z}^{(2)} (Ω, ω, Ω - ω) & χ_{y x z}^{(2)} (Ω, ω, Ω - ω) & χ_{y x y}^{(2)} (Ω, ω, Ω - ω) \\ χ_{z y z}^{(2)} (Ω, ω, Ω - ω) & χ_{z x z}^{(2)} (Ω, ω, Ω - ω) & χ_{z x y}^{(2)} (Ω, ω, Ω - ω) \end{matrix} \end{matrix}],

where

χ_{h j k}^{(2)}

is the second-order nonlinear susceptibility representing generation along the h coordinate due to excitations along the j and k coordinates. Notably, the subscripts h, j, and k can be either x, y or z. Each individual element in the nonlinear tensor produces a second-order nonlinear polarization according to the following equation,

P_{h : j, k}^{(2)} (Ω) = ε_{0} \int_{- \infty}^{\infty} χ_{h j k}^{(2)} (Ω, ω, Ω - ω) E_{j} (ω) E_{k} (Ω - ω) d ω,

where ε₀ is the permittivity of free space and

E_{j, k}

are the optical excitation electric fields having polarizations along the j and k coordinates, respectively. Moreover, the subscript h:j,k signifies that excitation frequencies oriented along the j and k coordinates produce a second-order nonlinear polarization along the h coordinate. To incorporate second-order nonlinear dispersive effects, the FDTD model implements

χ_{h j k}^{(2)}

coefficients according to Miller’s rule [1],

χ_{h j k}^{(2)} (Ω, ω, Ω - ω) = δ_{h j k} χ_{h} (Ω) χ_{j} (ω) χ_{k} (Ω - ω),

where

δ_{h j k}

is Miller’s proportionality constant and

χ_{h, j, k}

are the linear susceptibilities along the h, j, and k coordinates, respectively. It should be noted that Miller’s formalism inherently links the nonlinear dispersion of an optical material to the linear material dispersion, at both the excitation frequencies (i.e. ω and Ω-ω) and the generation frequency (i.e. Ω). Using the definition in Eq. (4) for

χ_{h j k}^{(2)}

, Eq. (3) becomes,

P_{h : j, k}^{(2)} (Ω) = ε_{0} δ_{h j k} χ_{h} (Ω) \int_{- \infty}^{\infty} S_{j} (ω) S_{k} (Ω - ω) d ω,

where,

S_{j, k} (ω) = χ_{j, k} (ω) E_{j, k} (ω) .

The integral in Eq. (5) is the convolution operation, such that it can be recasted in a more compact form,

P_{h : j, k}^{(2)} (Ω) = ε_{0} δ_{h j k} χ_{h} (Ω) {S_{j} * S_{k}} (Ω) .

By transforming

P_{h : j, k}^{(2)}

to the time-domain, we obtain,

P_{h : j, k}^{(2)} (t) = ε_{0} δ_{h j k} {χ_{h} * (S_{j} S_{k})} (t),

where t is time. Using Eq. (8), the complete set of second-order nonlinear polarizations along the x, y, and z directions are [14],

P_{x}^{(2)} (t) = P_{x : x, x}^{(2)} (t) + P_{x : y, y}^{(2)} (t) + P_{x : z, z}^{(2)} (t) + 2 P_{x : y, z}^{(2)} (t) + 2 P_{x : x, z}^{(2)} (t) + 2 P_{x : x, y}^{(2)} (t),

P_{y}^{(2)} (t) = P_{y : x, x}^{(2)} (t) + P_{y : y, y}^{(2)} (t) + P_{y : z, z}^{(2)} (t) + 2 P_{y : y, z}^{(2)} (t) + 2 P_{y : x, z}^{(2)} (t) + 2 P_{y : x, y}^{(2)} (t),

P_{z}^{(2)} (t) = P_{z : x, x}^{(2)} (t) + P_{z : y, y}^{(2)} (t) + P_{z : z, z}^{(2)} (t) + 2 P_{z : y, z}^{(2)} (t) + 2 P_{z : x, z}^{(2)} (t) + 2 P_{z : x, y}^{(2)} (t) .

The second-order nonlinear polarizations are related to the second-order nonlinear current densities through,

J_{x}^{(2)} (t) = \frac{d P_{x}^{(2)} (t)}{d t},

J_{y}^{(2)} (t) = \frac{d P_{y}^{(2)} (t)}{d t},

J_{z}^{(2)} (t) = \frac{d P_{z}^{(2)} (t)}{d t} .

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.

χ_{h, j, k}

) 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,

χ_{x, y, z} (ω) = χ_{s}^{x, y, z} + \sum_{m = 1}^{N^{x, y, z}} \frac{{(ω_{m}^{x, y, z})}^{2} χ_{m}^{x, y, z}}{{(ω_{m}^{x, y, z})}^{2} + i γ_{m}^{x, y, z} ω - ω^{2}},

where

χ_{s}^{x, y, z}

are the frequency-independent linear susceptibilities,

χ_{m}^{x, y, z}

are the Lorentz susceptibilities of the m^th Lorentzian oscillator,

ω_{m}^{x, y, z}

are the resonant frequencies of the m^th Lorentzian oscillator,

γ_{m}^{x, y, z}

are the damping factors of the m^th Lorentzian oscillator, and

N^{x, y, z}

are the number of Lorentzian oscillators used to describe the linear susceptibilities.

3. Discretization of the second-order nonlinear current density

To discretize the $J_{h}^{(2)}$ terms in Eqs. (12-14), we must first discretize the $S_{j, k}$ and $P_{h : j, k}^{(2)}$ terms in Eq. (6) and Eq. (7), respectively. Using Eq. (15), $S_{j, k}$ can be written as,

S_{j, k} (ω) = G_{s}^{j, k} (ω) + \sum_{m = 1}^{N^{j, k}} G_{m}^{j, k} (ω),

where,

G_{s}^{j, k} (ω) = χ_{s}^{j, k} E_{j, k} (ω),

and

G_{m}^{j, k} (ω) = \frac{{(ω_{m}^{j, k})}^{2} χ_{m}^{j, k}}{{(ω_{m}^{j, k})}^{2} + i γ_{m}^{j, k} ω - ω^{2}} E_{j, k} (ω) f o r m = 1, 2, ... N^{j, k} .

To solve the term in Eq. (17), it is first transformed to the time-domain,

G_{s}^{j, k} (t) = χ_{s}^{j, k} E_{j, k} (t) .

Using temporal averaging of the

G_{s}^{j, k}

term, Eq. (19) is then discretized for the time iteration of n to obtain,

G_{s}^{j, k} (n + 1) = 2 χ_{s}^{j, k} E_{j, k} (n) - G_{s}^{j, k} (n - 1),

where n is related to the time step, Δt, via the relationship t = nΔt. To solve Eq. (18), the terms are first rearranged to,

\begin{array}{l} (^{ω_{m}^{j, k}) 2} G_{m}^{j, k} (ω) + i γ_{m}^{j, k} ω G_{m}^{j, k} (ω) - ω^{2} G_{m}^{j, k} (ω) \\ = (^{ω_{m}^{j, k}) 2} χ_{m}^{j, k} E_{j, k} (ω) ​ ​ ​ ​ f o r m = 1, 2, ... N^{j, k} . \end{array}

Next, using the fact that

i ω G_{m}^{j, k} (ω) \Leftrightarrow \frac{d G_{m}^{j, k} (t)}{d t}

and

- ω^{2} G_{m}^{j, k} (ω) \Leftrightarrow \frac{d^{2} G_{m}^{j, k} (t)}{d t^{2}}

, Eq. (21) can be transformed to its time-domain form,

\begin{array}{l} (^{ω_{m}^{j, k}) 2} G_{m}^{j . k} (t) + γ_{m}^{j, k} \frac{d G_{m}^{j, k} (t)}{d t} + \frac{d^{2} G_{m}^{j, k} (t)}{d t^{2}} \\ = (^{ω_{m}^{j, k}) 2} χ_{m}^{j, k} E_{j, k} (t) ​ ​ ​ ​ ​ f o r m = 1, 2, ... N^{j, k} . \end{array}

Using central-differencing techniques, this equation is discretized for the time iteration n to obtain,

\begin{array}{l} G_{m}^{j, k} (n + 1) = \frac{2 Δ t^{2} {(ω_{m}^{j, k})}^{2} χ_{m}^{j, k}}{γ_{m}^{j, k} Δ t + 2} E_{j, k} (n) + \frac{4 - 2 Δ t^{2} {(ω_{m}^{j, k})}^{2}}{γ_{m}^{j, k} Δ t + 2} G_{m}^{j, k} (n) \\ + \frac{γ_{m}^{j, k} Δ t - 2}{γ_{m}^{j, k} Δ t + 2} G_{m}^{j, k} (n - 1) ​ ​ ​ ​ ​ f o r m = 1, 2, ... N^{j, k} . \end{array}

Now that both

G_{s}^{j, k}

and

G_{m}^{j, k}

have been retrieved in their discretized form, a solution for

S_{j, k}

is obtained by converting Eq. (16) to the time-domain and discretizing it for the time iteration of n + 1,

S_{j, k} (n + 1) = G_{s}^{j, k} (n + 1) + \sum_{m = 1}^{N^{j, k}} G_{m}^{j, k} (n + 1) .

Importantly,

S_{j, k} (n)

must also be obtained which leads to the following equations,

G_{s}^{j, k} (n) = χ_{s}^{j, k} E_{j, k} (n),

\begin{array}{l} G_{m}^{j, k} (n) = \frac{2 Δ t^{2} {(ω_{m}^{j, k})}^{2} χ_{m}^{j, k}}{γ_{m}^{j, k} Δ t + 2} E_{j, k} (n - 1) + \frac{4 - 2 Δ t^{2} {(ω_{m}^{j, k})}^{2}}{γ_{m}^{j, k} Δ t + 2} G_{m}^{j, k} (n - 1) \\ + \frac{γ_{m}^{j, k} Δ t - 2}{γ_{m}^{j, k} Δ t + 2} G_{m}^{j, k} (n - 2) ​ ​ ​ ​ ​ f o r m = 1, 2, ... N^{j, k}, \end{array}

S_{j, k} (n) = G_{s}^{j, k} (n) + \sum_{m = 1}^{N^{j, k}} G_{m}^{j, k} (n) .

Alternatively,

G_{s}^{j, k} (n)

and

G_{m}^{j, k} (n)

can instead be obtained from 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),

P_{h : j, k}^{(2)}

in Eq. (7) can be recasted as,

P_{h : j, k}^{(2)} (Ω) = K_{s}^{h j k} (Ω) + \sum_{m = 1}^{N^{h}} K_{m}^{h j k} (Ω),

where,

K_{s}^{h j k} (Ω) = ε_{0} δ_{h j k} χ_{s}^{h} {S_{j} * S_{k}} (Ω),

and

K_{m}^{h j k} (Ω) = ε_{0} δ_{h j k} \frac{{(ω_{m}^{h})}^{2} χ_{m}^{h}}{{(ω_{m}^{h})}^{2} + i γ_{m}^{h} Ω - Ω^{2}} ​ ​ ​ ​ ​ {S_{j} * S_{k}} (Ω) f o r m = 1, 2, ... N^{h} .

Equation (29) is solved by first converting it to the time-domain,

K_{s}^{h j k} (t) = ε_{0} δ_{h j k} χ_{s}^{h} S_{j} (t) S_{k} (t),

and then discretizing it for the time iteration of n + 1 to obtain,

K_{s}^{h j k} (n + 1) = ε_{0} δ_{h j k} χ_{s}^{h} S_{j} (n + 1) S_{k} (n + 1) .

Clearly, this depends on the solution to

S_{j, k}

, which is presented in Eq. (24). Equation (30) is solved by first rearranged it into the form,

\begin{array}{l} (^{ω_{m}^{h}) 2} K_{m}^{h j k} (Ω) + i γ_{m}^{h} Ω K_{m}^{h j k} (Ω) - Ω^{2} K_{m}^{h j k} (Ω) \\ = ε_{0} δ_{h j k} (^{ω_{m}^{h}) 2} χ_{m}^{h} {S_{j} * S_{k}} (Ω) ​ ​ ​ ​ ​ f o r m = 1, 2, ... N^{h} . \end{array}

Utilizing

i Ω K_{m}^{h j k} (Ω) \Leftrightarrow \frac{d K_{m}^{h j k} (t)}{d t}

and

- Ω^{2} K_{m}^{h j k} (Ω) \Leftrightarrow \frac{d^{2} K_{m}^{h j k} (t)}{d t^{2}}

with Eq. (33), we obtain the time-domain equation,

\begin{array}{l} (^{ω_{m}^{h}) 2} K_{m}^{h j k} (t) + γ_{m}^{h} \frac{d K_{m}^{h j k} (t)}{d t} + \frac{d^{2} K_{m}^{h j k} (t)}{d t^{2}} \\ = ε_{0} δ_{h j k} (^{ω_{m}^{h}) 2} χ_{m}^{h} S_{j} (t) S_{k} (t) ​ ​ ​ ​ ​ f o r m = 1, 2, ... N^{h} . \end{array}

Using central-differencing techniques, Eq. (34) is discretized for the time iteration of n,

\begin{array}{l} K_{m}^{h j k} (n + 1) = ε_{0} δ_{h j k} \frac{2 Δ t^{2} {(ω_{m}^{h})}^{2} χ_{m}^{h}}{γ_{m}^{h} Δ t + 2} S_{j} (n) S_{k} (n) + \frac{4 - 2 Δ t^{2} {(ω_{m}^{h})}^{2}}{γ_{m}^{h} Δ t + 2} K_{m}^{h j k} (n) \\ + \frac{γ_{m}^{h} Δ t - 2}{γ_{m}^{h} Δ t + 2} K_{m}^{h j k} (n - 1) ​ ​ ​ ​ ​ f o r m = 1, 2, ... N^{h}, \end{array}

where Eq. (35) depends on the

S_{j, k}

solution in Eq. (27). Now that

K_{s}^{h j k}

and

K_{m}^{h j k}

are obtained in their discretized form, a discretized solution can be attained for

P_{h : j, k}^{(2)}

by converting Eq. (28) to the time-domain and discretizing it for the time iteration of n + 1,

P_{h : j, k}^{(2)} (n + 1) = K_{s}^{h j k} (n + 1) + \sum_{m = 1}^{N^{h}} K_{m}^{h j k} (n + 1) .

The above equation is the final, discretized form for

P_{h : j, k}^{(2)}

, where

K_{s}^{h j k}

and

K_{m}^{h j k}

are given in Eqs. (32) and (35), respectively. Finally, the discretization of

P_{h}^{(2)}

[Eqs. (9-11)] for the time iteration of n + 1 is straightforward and results in,

\begin{array}{l} P_{x}^{(2)} (n + 1) = P_{x : x, x}^{(2)} (n + 1) + P_{x : y, y}^{(2)} (n + 1) + P_{x : z, z}^{(2)} (n + 1) \\ + 2 P_{x : y, z}^{(2)} (n + 1) + 2 P_{x : x, z}^{(2)} (n + 1) + 2 P_{x : x, y}^{(2)} (n + 1), \end{array}

\begin{array}{l} P_{y}^{(2)} (n + 1) = P_{y : x, x}^{(2)} (n + 1) + P_{y : y, y}^{(2)} (n + 1) + P_{y : z, z}^{(2)} (n + 1) \\ + 2 P_{y : y, z}^{(2)} (n + 1) + 2 P_{y : x, z}^{(2)} (n + 1) + 2 P_{y : x, y}^{(2)} (n + 1), \end{array}

\begin{array}{l} P_{z}^{(2)} (n + 1) = P_{z : x, x}^{(2)} (n + 1) + P_{z : y, y}^{(2)} (n + 1) + P_{z : z, z}^{(2)} (n + 1) \\ + 2 P_{z : y, z}^{(2)} (n + 1) + 2 P_{z : x, z}^{(2)} (n + 1) + 2 P_{z : x, y}^{(2)} (n + 1), \end{array}

whereas the discretization of

J_{h}^{(2)}

[Eqs. (12-14)] for the time iteration of n + 1/2 provides,

J_{x}^{(2)} (n + 1 / 2) = \frac{P_{x}^{(2)} (n + 1) - P_{x}^{(2)} (n)}{Δ t},

J_{y}^{(2)} (n + 1 / 2) = \frac{P_{y}^{(2)} (n + 1) - P_{y}^{(2)} (n)}{Δ t},

J_{z}^{(2)} (n + 1 / 2) = \frac{P_{z}^{(2)} (n + 1) - P_{z}^{(2)} (n)}{Δ t} .

To include second-order nonlinear effects in the FDTD formalism, these

J h^{(2)}

terms are incorporated in the update equations,

\begin{array}{l} E_{x}^{p - 1 / 2, r, s} (n + 1) = E_{x}^{p - 1 / 2, r, s} (n) - \frac{Δ t}{ε_{0} ε_{r x}} J'_{x}^{(1) p - 1 / 2, r, s} (n + \frac{1}{2}) - \frac{Δ t}{ε_{0} ε_{r x}} J_{x}^{(2) p - 1 / 2, r, s} (n + \frac{1}{2}) \\ + \frac{Δ t}{ε_{0} ε_{r x} Δ y} [H_{z}^{p - 1 / 2, r + 1 / 2, s} (n + \frac{1}{2}) - H_{z}^{p - 1 / 2, r - 1 / 2, s} (n + \frac{1}{2})] \\ - \frac{Δ t}{ε_{0} ε_{r x} Δ z} [H_{y}^{p - 1 / 2, r, s + 1 / 2} (n + \frac{1}{2}) - H_{y}^{p - 1 / 2, r, s - 1 / 2} (n + \frac{1}{2})], \end{array}

\begin{array}{l} E_{y}^{p, r - 1 / 2, s} (n + 1) = E_{y}^{p, r - 1 / 2, s} (n) - \frac{Δ t}{ε_{0} ε_{r y}} J'_{y}^{(1) p, r - 1 / 2, s} (n + \frac{1}{2}) - \frac{Δ t}{ε_{0} ε_{r y}} J_{y}^{(2) p, r - 1 / 2, s} (n + \frac{1}{2}) \\ + \frac{Δ t}{ε_{0} ε_{r y} Δ z} [H_{x}^{p, r - 1 / 2, s + 1 / 2} (n + \frac{1}{2}) - H_{x}^{p, r - 1 / 2, s - 1 / 2} (n + \frac{1}{2})] \\ - \frac{Δ t}{ε_{0} ε_{r y} Δ x} [H_{z}^{p + 1 / 2, r - 1 / 2, s} (n + \frac{1}{2}) - H_{z}^{p - 1 / 2, r - 1 / 2, s} (n + \frac{1}{2})], \end{array}

\begin{array}{l} E_{z}^{p, r, s - 1 / 2} (n + 1) = E_{z}^{p, r, s - 1 / 2} (n) - \frac{Δ t}{ε_{0} ε_{r z}} J'_{z}^{(1) p, r, s - 1 / 2} (n + \frac{1}{2}) - \frac{Δ t}{ε_{0} ε_{r z}} J_{z}^{(2) p, r, s - 1 / 2} (n + \frac{1}{2}) \\ + \frac{Δ t}{ε_{0} ε_{r z} Δ x} [H_{y}^{p + 1 / 2, r, s - 1 / 2} (n + \frac{1}{2}) - H_{y}^{p - 1 / 2, r, s - 1 / 2} (n + \frac{1}{2})] \\ - \frac{Δ t}{ε_{0} ε_{r z} Δ y} [H_{x}^{p, r + 1 / 2, s - 1 / 2} (n + \frac{1}{2}) - H_{x}^{p, r - 1 / 2, s - 1 / 2} (n + \frac{1}{2})], \end{array}

\begin{array}{l} H_{x}^{p, r - 1 / 2, s - 1 / 2} (n + \frac{1}{2}) = H_{x}^{p, r - 1 / 2, s - 1 / 2} (n - \frac{1}{2}) \\ + \frac{Δ t}{μ_{0} Δ z} [E_{y}^{p, r - 1 / 2, s} (n) - E_{y}^{p, r - 1 / 2, s - 1} (n)] - \frac{Δ t}{μ_{0} Δ y} [E_{z}^{p, r, s - 1 / 2} (n) - E_{z}^{p, r - 1, s - 1 / 2} (n)], \end{array}

\begin{array}{l} H_{y}^{p - 1 / 2, r, s - 1 / 2} (n + \frac{1}{2}) = H_{y}^{p - 1 / 2, r, s - 1 / 2} (n - \frac{1}{2}) \\ + \frac{Δ t}{μ_{0} Δ x} [E_{z}^{p, r, s - 1 / 2} (n) - E_{z}^{p - 1, r, s - 1 / 2} (n)] - \frac{Δ t}{μ_{0} Δ z} [E_{x}^{p - 1 / 2, r, s} (n) - E_{x}^{p - 1 / 2, r, s - 1} (n)], \end{array}

\begin{array}{l} H_{z}^{p - 1 / 2, r - 1 / 2, s} (n + \frac{1}{2}) = H_{z}^{p - 1 / 2, r - 1 / 2, s} (n - \frac{1}{2}) \\ + \frac{Δ t}{μ_{0} Δ y} [E_{x}^{p - 1 / 2, r, s} (n) - E_{x}^{p - 1 / 2, r - 1, s} (n)] - \frac{Δ t}{μ_{0} Δ x} [E_{y}^{p, r - 1 / 2, s} (n) - E_{y}^{p - 1, r - 1 / 2, s} (n)], \end{array}

where

μ_{0}

is the permeability free space, p, r, and s are the spatial indices of the field components along the x, y, and z-directions, respectively, Δx, Δy, and Δz are the mesh steps along the x, y, and z-directions, respectively, H_x, H_y, and H_z are the magnetic field components along the x, y, and z-directions, respectively,

ε_{r x}

,

ε_{r y}

, and

ε_{r z}

are the relative permittivities along the x, y, and z-directions, respectively, and

J'_{x}^{(1)}

,

J'_{y}^{(1)}

, and

J'_{z}^{(1)}

are the dispersive parts of the linear components of the current densities along the x, y, and z-directions, respectively.

4. Frequency-conversion in a LiNbO₃ crystal

The derived FDTD method is evaluated by simulating the representative effects of optical rectification and DFG. These nonlinear effects are investigated using a LiNbO₃ crystal layer having a length l, which is illustrated in Fig. 1. LiNbO₃ 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 LiNbO₃ 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 LiNbO₃ layer, respectively. The c-axis of the crystal is oriented along the z-axis, the [100] 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.

Fig. 1 An illustration of the LiNbO₃ crystal having a length l and an excitation electric field polarization at the angle θ with respect to the crystal c-axis.

Download Full Size | PDF

It is critical to determine the frequency dependence of the LiNbO₃ refractive indices and extinction coefficients, along with its second-order nonlinear susceptibilities. The extraordinary, n_e, and ordinary, n_o, 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. n_e and к_e show a phonon resonance at 7.6 THz, corresponding to the lowest-frequency A₁ mode of LiNbO₃, whereas three E mode phonon resonances are observed in n_o and к_o, which are located at 4.6, 7.9, and 9.7 THz [15]. Since nonlinear dispersion is directly related to linear dispersion via Miller’s rule, the $χ^{(2)}$ elements of the LiNbO₃ crystal cannot be taken as constants. The LiNbO₃ crystal belongs to the 3m point group symmetry class, such that its second-order nonlinear susceptibly tensor is written as [11],

\bar{\bar{χ^{(2)}}} (Ω, ω, Ω - ω) = [\begin{matrix} \begin{matrix} 0 & 0 & 0 \\ - χ_{22}^{(2)} (Ω, ω, Ω - ω) & χ_{22}^{(2)} (Ω, ω, Ω - ω) & 0 \\ χ_{31}^{(2)} (Ω, ω, Ω - ω) & χ_{31}^{(2)} (Ω, ω, Ω - ω) & χ_{33}^{(2)} (Ω, ω, Ω - ω) \end{matrix} & \begin{matrix} 0 & χ_{15}^{(2)} (Ω, ω, Ω - ω) & - χ_{22}^{(2)} (Ω, ω, Ω - ω) \\ χ_{15}^{(2)} (Ω, ω, Ω - ω) & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \end{matrix}],

where

χ_{15}^{(2)} (Ω, ω, Ω - ω)

,

χ_{22}^{(2)} (Ω, ω, Ω - ω)

,

χ_{31}^{(2)} (Ω, ω, Ω - ω) c

, and

χ_{33}^{(2)} (Ω, ω, Ω - ω)

are the non-vanishing coefficients. Notably, since the THz region exhibits significant dispersion, the Kleinman symmetry condition is invalid and the

χ_{15}^{(2)} (Ω, ω, Ω - ω)

element differs from the

χ_{31}^{(2)} (Ω, ω, Ω - ω)

element. From 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,

P_{x}^{(2)} (Ω) = 2 ε_{0} \int_{- \infty}^{\infty} χ_{15}^{(2)} (Ω, ω, Ω - ω) E_{x} (ω) E_{z} (Ω - ω) d ω,

and,

\begin{array}{l} P_{z}^{(2)} = ε_{0} \int_{- \infty}^{\infty} χ_{31}^{(2)} (Ω, ω, Ω - ω) E_{x} (ω) E_{x} (Ω - ω) d ω \\ + ε_{0} \int_{- \infty}^{\infty} χ_{33}^{(2)} (Ω, ω, Ω - ω) E_{z} (ω) E_{z} (Ω - ω) d ω, \end{array}

such that it is necessary to consider dispersion for the

χ_{15}^{(2)} (Ω, ω, Ω - ω)

,

χ_{31}^{(2)} (Ω, ω, Ω - ω)

, and

χ_{33}^{(2)} (Ω, ω, Ω - ω)

tensor elements. Figure 3(a) illustrates the highly-dispersive

χ_{33}^{(2)} (Ω, ω, Ω - ω)

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

χ_{31}^{(2)} (Ω, ω, Ω - ω)

element [see Fig. 3(b)] is obtained by using a proportionality constant of δ = 0.21 pm/V and the experimental data in [19], whereas

χ_{15}^{(2)} (Ω, ω, Ω - ω)

[see Fig. 3(c)] is described by δ = 0.71 pm/V and the experimental data from [19]. 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.

Fig. 2 (a) Refractive index and (b) extinction coefficient for the extraordinary LiNbO₃ crystal direction in the THz frequency regime. (c) Refractive index and (d) extinction coefficient for the ordinary LiNbO₃ crystal direction in the THz frequency regime. (e) Extraordinary and ordinary refractive indices in the optical frequency regime. The experimental data is obtained from [15,16] and the LiNbO₃ crystal is taken to be lossless in the optical regime.

Download Full Size | PDF

Fig. 3 Magnitude and phase of the second-order nonlinear susceptibility elements of LiNbO₃ for (a) $χ_{33}^{(2)} (Ω, ω, Ω - ω)$ , (b) $χ_{31}^{(2)} (Ω, ω, Ω - ω)$ , and (c) $χ_{15}^{(2)} (Ω, ω, Ω - ω)$ as a function of frequency. The curve fits are calculated using Miller’s rule with the experimental data from [17–19].

Download Full Size | PDF

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 $χ_{33}^{(2)} (Ω, ω, Ω - ω)$ . The simulations are conducted using the fully-dispersive $χ_{33}^{(2)} (Ω, ω, Ω - ω)$ element [see Fig. 3(a)], as well as the frequency-independent $χ_{33}^{(2)}$ = 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,

f_{d} = \frac{c q}{l | n_{g}^{o p t} - n_{T H z} (f_{d}) |},

where c is the speed of light, q is a positive integer,

n_{g}^{o p t}

is the extraordinary group refractive index at the wavelength of 1550 nm, and n_THz is the extraordinary phase refractive index in the THz regime. When comparing ξ_z obtained using the dispersive

χ_{33}^{(2)} (Ω, ω, Ω - ω)

and the frequency-independent

χ_{33}^{(2)}

= 348 pm/V, significant disagreement is evident, especially at frequencies near the phonon resonance of LiNbO₃ at 7.6 THz. Figure 4(b) shows the relative spectral power, which is defined as ξ_z obtained using the dispersive

χ_{33}^{(2)} (Ω, ω, Ω - ω)

divided by ξ_z obtained using the frequency independent

χ_{33}^{(2)}

= 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.

Fig. 4 (a) The z-component of the spectral power obtained using the dispersive $χ_{33}^{(2)} (Ω, ω, Ω - ω)$ and a frequency-independent $χ_{33}^{(2)}$ = 348 pm/V. (b) Relative spectral power calculated when implementing the dispersive $χ_{33}^{(2)} (Ω, ω, Ω - ω)$ and the frequency-independent $χ_{33}^{(2)}$ = 348 pm/V.

Download Full Size | PDF

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, f_p, the signal frequency, f_s, and either the forward, $f_{i}^{+}$ , or backward, $f_{i}^{-}$ , propagating idler frequencies. The forward propagation DFG coherence length, $L_{c}^{+}$ , is,

L_{c}^{+} = \frac{c}{2} {| n_{p} f_{p} - n_{s} f_{s} - n_{i} f_{i}^{+} |}^{- 1},

and the backward propagation DFG coherence length,

L_{c}^{-}

, is,

L_{c}^{-} = \frac{c}{2} {| n_{p} f_{p} - n_{s} f_{s} + n_{i} f_{i}^{-} |}^{- 1},

where n_p, n_s, and n_i are the pump, signal, and idler refractive indices, respectively, and

f_{i} = | f_{p} - f_{s} |

. Since the excitation angle is at θ = 45°, the electric fields at f_p, f_s, and

f_{i}^{+}

(or

f_{i}^{-}

) are polarized along the ordinary, extraordinary, and ordinary crystal directions, respectively. When considering THz radiation generation in the forward direction [i.e. Equation (53)], phase-matching occurs at

f_{i}^{+}

= 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

f_{i}^{+}

= 3 THz, as seen in Fig. 5(b). Similarly,

L_{c}^{-}

from Eq. (54) is displayed in Fig. 5(c), which shows that phase-matching occurs at the idler frequency of

f_{i}^{-}

= 1.8 THz. By performing FDTD simulations and recording ξ_x near the input facet of the LiNbO₃ layer, we confirm that

f_{i}^{-}

= 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.

Fig. 5 (a) Coherence length for THz radiation generation in the forward direction. (b) The x-component of the spectral power recorded in free space. (c) Coherence length for THz radiation generation in the backward direction. (d) The x-component of the spectral power recorded near the input face of the LiNbO₃ crystal.

Download Full Size | PDF

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 $P_{z}^{(2)}$ is maximum and the only contribution is from the $χ_{33}^{(2)} (Ω, ω, Ω - ω)$ term [see Eq. (51)]. THz radiation generation decreases as θ increases, since the contribution from $χ_{33}^{(2)} (Ω, ω, Ω - ω)$ is reduced; nonetheless, the weaker $χ_{31}^{(2)} (Ω, ω, Ω - ω)$ coefficient now contributes to $P_{z}^{(2)}$ . At θ = 90°, THz radiation generation is lowest, since $P_{z}^{(2)}$ is only influenced by $χ_{31}^{(2)} (Ω, ω, Ω - ω)$ . Figure 6(b) depicts ξ_x, where the THz radiation generation vanishes at both θ = 0° and 90°, since $P_{x}^{(2)}$ = 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.

Fig. 6 (a) The z-component of the spectral power and (b) the x-component of the spectral power, at various electric field polarization angles.

Download Full Size | PDF

5. Summary

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 LiNbO₃ 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 LiNbO₃ 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 [20], or the generation of radiation in bulk materials that are excited by a focused laser beam.

Funding

Natural Sciences and Engineering Research Council of Canada.

Acknowledgments

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.

References

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]

4. A. Mayer and F. Keilmann, “Far-infrared nonlinear optics. I. χ ⁽²⁾ near ionic resonance,” Phys. Rev. B Condens. Matter 33(10), 6954–6961 (1986). [CrossRef] [PubMed]

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 d₁₁ 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, “LiNbO₃: 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 LiNbO₃ 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 LiNbO₃,” 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).

20. A. Cala’ Lesina, L. Ramunno, and P. Berini, “Dual-polarization plasmonic metasurface for nonlinear optics,” Opt. Lett. 40(12), 2874–2877 (2015). [CrossRef] [PubMed]

A modeling of dispersive tensorial second-order nonlinear effects for the finite-difference time-domain method

Abstract

1. Introduction

2. Derivation of the second-order nonlinear current density

3. Discretization of the second-order nonlinear current density

4. Frequency-conversion in a LiNbO₃ crystal

5. Summary

Funding

Acknowledgments

References

Cited By

Figures (6)

Equations (54)

Optics Express

Abstract

1. Introduction

2. Derivation of the second-order nonlinear current density

3. Discretization of the second-order nonlinear current density

4. Frequency-conversion in a LiNbO3 crystal

5. Summary

Funding

Acknowledgments

References

Cited By

Figures (6)

Equations (54)

Optics Express

4. Frequency-conversion in a LiNbO₃ crystal