*Article* **Investigation of Mid-Infrared Broadband Second-Harmonic Generation in Non-Oxide Nonlinear Optic Crystals**

**Ilhwan Kim, Donghwa Lee and Kwang Jo Lee \***

Department of Applied Physics, Institute of Natural Science, Kyung Hee University, Yongin-si 17104, Korea; dlfghks383@gmail.com (I.K.); fairytale095@gmail.com (D.L.)

**\*** Correspondence: kjlee88@khu.ac.kr

**Abstract:** The mid-infrared (mid-IR) continuum generation based on broadband second harmonic generation (SHG) (or difference frequency generation) is of great interest in a wide range of applications such as free space communications, environmental monitoring, thermal imaging, highsensitivity metrology, gas sensing, and molecular fingerprint spectroscopy. The second-order nonlinear optic (NLO) crystals have been spotlighted as a material platform for converting the wavelengths of existing lasers into the mid-IR spectral region or for realizing tunable lasers. In particular, the spectral coverage could be extended to ~19 μm with non-oxide NLO crystals. In this paper, we theoretically and numerically investigated the broadband SHG properties of non-oxide mid-IR crystals in three categories: chalcopyrite semiconductors, defect chalcopyrite, and orthorhombic ternary chalcogenides. The technique is based on group velocity matching between interacting waves in addition to birefringent phase matching. We will describe broadband SHG characteristics in terms of beam propagation directions, spectral positions of resonance, effective nonlinearities, spatial walk-offs between interacting beams, and spectral bandwidths. The results will show that the spectral bandwidths of the fundamental wave allowed for broadband SHG to reach several hundreds of nm. The corresponding SH spectral range spans from 1758.58 to 4737.18 nm in the non-oxide crystals considered in this study. Such broadband SHG using short pulse trains can potentially be applied to frequency up-conversion imaging in the mid-IR region, in information transmission, and in nonlinear optical signal processing.

**Keywords:** mid-infrared photonics; chalcopyrite semiconductors; orthorhombic ternary chalcogenides; harmonic generation; continuum generation; group velocity matching

#### **1. Introduction**

The field of mid-infrared (IR) photonics is growing rapidly due to increasing demand for applications such as free space communications, remote sensing, environmental monitoring, thermal imaging, defense, IR countermeasure, medicine, gas sensing, and molecular fingerprint spectroscopy [1–3]. Nonlinear optic (NLO) crystals—particularly with secondorder nonlinearities—have been spotlighted as a material platform for converting the wavelengths of existing lasers into the mid-IR spectral region or for realizing tunable lasers [4,5]. The mid-IR spectral regions up to ~4 μm could be readily accessible with the oxide NLO crystals such as lithium niobate (LiNbO3), lithium tantalate (LiTaO3), and potassium titanyl phosphate (KTiOPO4) [6–9]. A high peak-power mid-IR optical amplifier using a potassium titanyl arsenate (KTiOAsO4) is also reported in [10]. These oxide crystals have been mainly used in the near-IR region but are still transparent within 4 μm. However, the spectral coverage of oxide NLO crystals is limited to less than 5 μm due to multi-phonon absorption [4]. The upper spectral limit of photons generated via parametric generation can be extended to ~19 μm using non-oxide crystals such as chalcopyrite semiconductors, orthorhombic ternary chalcogenides, and orientation-patterned (OP) semiconductors (e.g., OP-GaAs, OP-GaP, OP-ZnSe, and OP-GaN) [4,5,11,12]. The OP semiconductors with

**Citation:** Kim, I.; Lee, D.; Lee, K.J. Investigation of Mid-Infrared Broadband Second-Harmonic Generation in Non-Oxide Nonlinear Optic Crystals. *Crystals* **2021**, *11*, 921. https://doi.org/10.3390/ cryst11080921

Academic Editors: Xiaoming Duan, Renqin Dou, Linjun Li and Xiaotao Yang

Received: 12 July 2021 Accepted: 6 August 2021 Published: 8 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

periodic inversion of crystalline orientation use quasi-phase matching (QPM) for NLO interactions, whereas chalcopyrite semiconductors and orthorhombic ternary chalcogenides generally utilize birefringent phase matching (BPM). Extensive experimental studies of mid-IR laser sources using chalcopyrite semiconductors (e.g., silver thiogallate (AgGaS2, AGS), silver gallium selenide (AgGaSe2, AGSe), cadmium silicon phosphide (CdSiP2, CSP), and zinc germanium phosphide (ZnGeP2, ZGP)), defect chalcopyrite crystals (e.g., mercury thiogallate (HgGa2S4, HGS) and cadmium selenide (CdSe)), and orthorhombic ternary chalcogenides (e.g., lithium thioindate (LiInS2, LIS), lithium thiogallate (LiGaS2, LGS), and lithium gallium selenide (LiGaSe2, LGSe)) have been reported continuously in recent years [13–26].

Continuum generation based on broadband second harmonic generation (SHG) (or difference frequency generation) is of great interest in a wide range of applications. Considering the bandwidth of the continuum is inversely proportional to the group velocity (GV) mismatch between the interacting optical waves, broadband parametric generation is possible through GV matching [12]. Such broadband SHG using short pulse trains can potentially be applied to frequency up-conversion imaging in the mid-IR region, in information transmission using optical pulse signals, and in nonlinear optical signal processing [27–29]. For example, a wave packet transmitted through a free-space communication system can be decomposed into individual pulse signals by Fourier analysis and then each optical signal containing information is modulated and controlled for by another clock pulse train [27]. Recently, the potential to create such a continuum using ultrashort-pulse lasers has been estimated, especially in oxide crystals [12,30]. However, for non-oxide NLO crystals, broadband parametric generation has not yet been intensively studied. In this paper, we theoretically and numerically investigate broadband SHGs in three categories of commercially available non-oxide NLO crystals: (1) chalcopyrite semiconductors exhibiting uniaxial birefringence (e.g., cadmium germanium arsenide (CdGeAs2, CGA), AGS, AGSe, CSP, and ZGP); (2) defect chalcopyrite showing uniaxial birefringence (e.g., thallium arsenic selenide (Tl3AsSe3, TASe), gallium selenide (GaSe), CdSe, and HGS); and (3) orthorhombic ternary chalcogenides showing biaxial birefringence (e.g., lithium indium selenide (LiInSe2, LISe), LIS, LGS, and LGSe). The technique is based on GV matching between interacting waves in addition to BPM. In this simultaneous BPM–GV matching approach, the GV mismatch (GVM) between the interacting waves is always zero at the resonances within a specific spectral range that satisfies both the BPM and zero GVM. For uniaxial crystals, the wavelengths satisfying the simultaneous BPM–GV matching scheme correspond to the "magic" wavelengths in oxide crystals as reported in [12,30]. For biaxial crystals, these magic wavelengths extend into specific magic spectral regions, as will be discussed later. The advantage of this scheme is that the acceptable bandwidth of the fundamental (F) wave becomes very broad. This broadband spectrum has potential applications in multi-channel nonlinear optic signal processing. A broader input bandwidth also means that a train of pulses with narrower temporal widths can be used as the F-wave because the temporal width of a pulse has a Fourier transform relationship with its spectral width. Another critical advantage of the simultaneous BPM–GV matching scheme is that it allows for the use of long crystal lengths without considering the temporal walk-off between the interacting waves. The NLO efficiency is proportional to the square of the crystal length. Of course, there is still a limit on the crystal length due to the spatial walk-off between the interacting waves but this is significantly smaller than the limit due to the GVM, as will be described in Section 3. For each kind of non-oxide crystal, Type I and Type II NLO interactions will be considered, in which the polarizations of two fundamental photons to produce a second harmonic (SH) photon are either parallel to each other (for Type I) or perpendicular to each other (for Type II). We will describe the broadband SHG characteristics of each crystal in terms of beam propagation directions, spectral positions of resonance, effective nonlinearities, spatial walk-offs between interacting beams, and spectral bandwidths. The results will show that the spectral bandwidths of the F-wave allowed for the broadband SHG to reach several hundreds of nm. The SH spectral range satisfying

the simultaneous BPM–GV matchings span from 1758.58 to 4737.18 nm in the non-oxide crystals considered in this study. Finally, potential applications of such broadband light sources will be briefly described.

#### **2. Materials and Theories**

The non-oxide NLO crystals considered in this work are listed in Table 1. As can be seen from Table 1, chalcopyrite semiconductors and defective chalcopyrite exhibit either positive or negative uniaxial birefringence, whereas orthorhombic ternary chalcogenides have biaxial birefringence. As will be described later, the type of birefringence determines the angle tunability of the spectral position of NLO resonance. The point groups for each crystal determines the effective nonlinearity for the given direction of the input F-beam and are listed together in Table 1. The transparent range for each crystal is summarized in Table 1 along with references. In this section, we will describe the theoretical details of the BPM and GV matching properties of non-oxide NLO crystals, the effective nonlinearities, and the spatial walk-offs between the interacting beams.

**Table 1.** Non-oxide nonlinear optic (NLO) crystals considered for broadband mid-IR second-harmonic generation (SHG).


#### *2.1. BPM and GV Matching for Broadband SHG*

Figure 1 illustrates the polarization relationships of F and SH waves. H and V in Figure 1 represent horizontal and vertical polarization directions, respectively. For Type I, a pair of F photons with the same frequencies, *ω,* and polarization states produces an SH photon with a frequency, 2*ω*. In this case, the SH photon and F photons have polarization states perpendicular to each other (Figure 1a). For Type II, a pair of F photons with H and V polarization states, respectively, generates an SH photon (Figure 1b). In the experiment, a pair of photons with H and V polarization states were obtained from the input light polarized at 45◦ to the horizontal direction, as shown in Figure 1b.

In uniaxial crystals, the order of the refractive index (RI) magnitude is given by either *ne* < *no* (for negative uniaxial crystals: AGS, AGSe, CSP, HGS, TASe, and GaSe) or *ne* > *no* (for positive uniaxial crystals: CGA, ZGP, and CdSe), where *no* and *ne* represent the RIs of the ordinary (*o*) and extraordinary (*e*) waves in the principal axes of an index ellipsoid [42]. Then, the collinear BPM conditions can be expressed as:

**Figure 1.** Schematic diagrams showing the polarization relationships of fundamental (F) and second harmonic (SH) waves for (**a**) Type I and (**b**) Type II. V and H represent vertical and horizontal polarization directions, respectively.

$$\Delta k = |k\_\varepsilon(2\omega, \theta) - 2k\_\vartheta(\omega)| = 0 \text{ (for type I, negative uniaxial)},\tag{1}$$

$$\Delta k = |k\_t(2\omega, \theta) - k\_o(\omega) - k\_t(\omega, \theta)| \, = \text{ (for type II, negative uniaxial)}, \tag{2}$$

$$
\Delta k = |k\_{\vartheta}(2\omega) - 2k\_{\mathfrak{c}}(\omega, \theta)| \\
= 0 \text{ (for type I, positive unitaxial)}, \tag{3}
$$

$$
\Delta k = |k\_o(2\omega) - k\_v(\omega) - k\_\varepsilon(\omega, \theta)| \\
= 0 \text{ (for type II, positive uniaxial)}, \tag{4}
$$

where each *k* represents the wave number of the interacting wave and is defined as *ko*(*jω*)=(*jω*/*c*)*no*(*jω*) or *ke*(*jω*, *θ*)=(*jω*/*c*)*ne*(*jω*, *θ*). Here, *j* can be 1 or 2 and *c* denotes the speed of light in vacuum. For light propagating at the angle *θ* to the optic axis of uniaxial birefringence, the RI of *e*-wave, *ne*(*θ*), can be derived as follows using the definition in [42]:

$$n\_{\varepsilon}(j\omega,\theta) = \frac{n\_{\sigma}(j\omega)n\_{\varepsilon}(j\omega)}{\sqrt{n\_{\sigma}^{2}(j\omega)\sin^{2}\theta + n\_{\varepsilon}^{2}(j\omega)\cos^{2}\theta}}.\tag{5}$$

The temporal walk-off between the interacting waves due to the difference in GV can be defined as the time delay (Δ*T*) per unit crystal length as follows:

$$\frac{\Delta T}{L} = \frac{\Delta n\_{\mathcal{S}}}{c},\tag{6}$$

where *L* and Δ*ng* represent the crystal length and the group index difference between interacting photons, respectively. When Δ*ng* = 0 in Equation (6), GV matching is achieved, which can be simplified as:

$$n\_{\mathfrak{e}}^{\langle \mathfrak{g} \rangle}(2\omega, \theta) = n\_{\mathfrak{o}}^{\langle \mathfrak{g} \rangle}(\omega) \text{ (for type I, negative uniaxial)},\tag{7}$$

$$2n\_{\varepsilon}^{(\mathbf{g})}(2\omega,\theta) = n\_{\mathbf{o}}^{(\mathbf{g})}(\omega) + n\_{\varepsilon}^{(\mathbf{g})}(\omega,\theta) \text{ (for type II, negative uniaxial)},\tag{8}$$

$$n\_o^{(\mathfrak{g})}(2\omega) = n\_c^{(\mathfrak{g})}(\omega, \theta) \text{ (for type I, positive uniaxial), and} \tag{9}$$

$$2n\_o^{(g)}(2\omega) = n\_o^{(g)}(\omega) + n\_\varepsilon^{(g)}(\omega, \theta) \text{ (for type II, positive uniaxial)},\tag{10}$$

where each superscript *g* in Equations (7)–(10) indicates the group index. Now, the broadband SHG in a negative uniaxial crystal is defined as Equations (1) and (7) (for Type I) or Equations (2) and (8) (for Type II), being satisfied simultaneously. For a positive uniaxial crystal, the broadband SHG condition is defined as Equations (3) and (9) (for Type I) or Equations (4) and (10) (for Type II). Table 2 summarizes the BPM and GV matching conditions for broadband SHG in uniaxial crystals. All equations in Table 2 are expressed as two-variable functions for *θ* and *ω* (or for the F-wavelength, *λF*). Solving the system of equations for BPM and GV matching yields a solution set of *λ<sup>F</sup>* and *θ*, which determine the center wavelength of broadband SHG and the propagation direction of the F-beam.


**Table 2.** The BPM and GV matching conditions for broadband SHG in uniaxial crystals.

The orthorhombic ternary chalcogenides considered in this study (i.e., LIS, LISe, LGS, and LGSe) are all negative biaxial crystals belonging to the point group of orthorhombic *mm*2 at room temperature (see Table 1). Their crystallographic axes (*a*, *b*, and *c*) are all perpendicular to each other and have a relationship with the optical axes, where (*y*, *x*, *z*)=(*a*, *b*, *c*), within the spectral region considered in this study [39,43–45]. In this assignment, the order of RI magnitudes is given by *nx* < *ny* < *nz*. Then, the collinear BPM conditions can be expressed as:

$$
\Delta k = \left| k^{(l)}(2\omega) - 2k^{(h)}(\omega) \right| \\
= 0 \text{ (for type I, binomial) and} \\
\tag{11}
$$

$$
\Delta k = \left| k^{(l)}(2\omega) - k^{(h)}(\omega) - k^{(l)}(\omega) \right| \\
= \left. 0 \left( \text{for type II, binomial} \right) , \tag{12}
$$

where each *k* represents the wave number of the interacting wave and is defined as *k*(*m*) (*jω*)=(*jω*/*c*)*n*(*m*) . The unit *k*-vector is defined in the spherical coordinate as, (sin*θ*cos*ϕ*, sin*θ*sin*ϕ*, and cos*θ*). Here, *θ* and *ϕ* are the polar and azimuthal angles, respectively. The RIs of the two eigen-polarization modes of light traveling inside a biaxial crystal are expressed as follows by solving the Fresnel equation of the wave normal [46]:

$$n^{(m)}(j\omega) \, = \sqrt{\frac{2}{-B\_j \pm \sqrt{B\_j^2 - 4C\_j}}} \,. \tag{13}$$

The parameters in Equation (13) are defined as follows:

$$B\_{\dot{j}} = -(b\_{\dot{j}} + c\_{\dot{j}})k\_x^2 - (a\_{\dot{j}} + c\_{\dot{j}})k\_y^2 - (b\_{\dot{j}} + a\_{\dot{j}})k\_{z\prime}^2 \tag{14}$$

$$\mathbf{C}\_{\rangle} = \,^b \mathbf{c}\_{\rangle} \mathbf{c}\_{\rangle}^2 + a\_{\rangle} \mathbf{c}\_{\rangle} \mathbf{k}\_{y}^2 + b\_{\rangle} a\_{\rangle} \mathbf{k}\_{z\prime}^2 \text{ and} \tag{15}$$

$$a\_{\rangle} = \ n\_x(j\omega), \ b\_{\rangle} = \ n\_y(j\omega), \ c\_{\rangle} = \ n\_z(j\omega), \tag{16}$$

where each *ki* represents the *x*, *y*, and *z*-axis components of the wave vector. *j* can be 1 or 2 and then *n*(*m*) (*ω*) and *n*(*m*) (2*ω*) denote the RIs of the F-wave and the SH wave with frequencies *ω* and 2*ω*, respectively. *m* in Equation (13) can be either *l* or *h*, representing low or high RI. *l* and *h* are obtained by taking the plus and minus from the ± sign of the denominator in Equation (13), respectively. The GV matching conditions Δ*ng* = 0 obtained from Equation (6) are simplified as:

$$n\_{\mathcal{J}}^{(l)}(2\omega) = n\_{\mathcal{J}}^{(h)}(\omega) \text{ (for type I, biaxial) and} \tag{17}$$

$$2n\_{\mathcal{S}}^{(l)}(2\omega) \;=\; n\_{\mathcal{S}}^{(h)}(\omega) + n\_{\mathcal{S}}^{(l)}(\omega) \; (\text{for type II, binomial}).\tag{18}$$

Each subscript *g* in Equations (17) and (18) indicates the group index. Now, the broadband SHG in a biaxial crystal is defined as Equations (11) and (17) (for Type I) or Equations (12) and (18) (for Type II), being satisfied simultaneously. The BPM and GV matching conditions for broadband SHG in biaxial crystals are listed in Table 3. Each of these equations in Table 3 are a function of three variables: F-wavelength (*λF*), *θ*, and *ϕ*. Therefore, by solving the system of equations for BPM and GV matching while changing *λF*, we can get a set of solutions for *θ* and *ϕ*, i.e., the direction of the F-wave vector for

broadband SHG at the given *λF*. Thus, the resonant *λ<sup>F</sup>* can be continuously tuned by sweeping the F-wave vector along the direction characterized by the solution set of *θ* and *ϕ*. In other words, for the orthorhombic ternary chalcogenide considered in this study, the spectral position of the SH wave can be selectively determined or tuned within the range of the solution sets. In contrast, for chalcopyrite crystals exhibiting uniaxial birefringence, SH resonance can only be achieved at a single wavelength, as discussed earlier.

**Type Condition Negative Biaxial** Type I BPM *<sup>n</sup>*(*l*)(2*ω*) = *<sup>n</sup>*(*h*)(*ω*) GV matching *<sup>n</sup>*(*l*) *<sup>g</sup>* (2*ω*) = *<sup>n</sup>*(*h*) *<sup>g</sup>* (*ω*) Type II BPM <sup>2</sup>*n*(*l*)(2*ω*) = *<sup>n</sup>*(*h*)(*ω*) + *<sup>n</sup>*(*l*)(*ω*) GV matching <sup>2</sup>*n*(*l*) *<sup>g</sup>* (2*ω*) = *<sup>n</sup>*(*h*) *<sup>g</sup>* (*ω*) + *<sup>n</sup>*(*l*) *<sup>g</sup>* (*ω*)

**Table 3.** The BPM and GV matching conditions for broadband SHG in biaxial crystals.

#### *2.2. Effective Nonlinearities*

The effective nonlinearity of a crystal depends on its point group and the NLO interaction type [47,48]. The point groups of chalcopyrite and defect chalcopyrite crystals considered in this study are 42*m*, 4, 3*m*, 62*m*, and 6*mm* (see Table 1). The corresponding analytical expressions of the effective NLO coefficients (*deff*) are listed in column 5 of Table 4. In each *deff* expression, *θ* and *ϕ* represent the polar and azimuthal angles, respectively. The spatial walk-off angle (*ρ*) between the wave vector and the Poynting vector within a uniaxial crystal is derived as:

$$\tan \rho = \frac{\sin 2\theta}{2} \left( \frac{1}{n\_c^2(\omega)} - \frac{1}{n\_o^2(\omega)} \right) n\_c^2(\omega, \theta), \tag{19}$$

where *ne*(*ω*,*θ*) is given in Equation (5). To derive this expression, we used the definition in [49]. Note that in each *deff* expression (column 5 of Table 4), *θ* is corrected as much as *ρ* (i.e., *θ* → *θ+ρ*). As the interacting beams do not usually propagate along the crystallographic axis due to birefringence, there is a spatial walk-off between the beams even in the case of collinear BPM. Thus, we need to correct the angle from *θ* to *θ+ρ* to obtain *deff* accordingly. Considering the azimuthal angle has no effect on both BPM and GV matching conditions as shown in Table 2, *ϕ* can be chosen as an arbitrary value to maximize *deff*. As can be appreciated from Equation (19), *ρ* is also a function of *λ<sup>F</sup>* and *θ*. Therefore, for uniaxial crystals, *deff* can be obtained by substituting the solution sets of *λ<sup>F</sup>* and *θ* that satisfy broadband SHG conditions (i.e., Table 2). The NLO efficiency is proportional to the square of *deff* for a given direction of beam propagation [42]. The *dil* components on the crystallographic axes were chosen as the measurements closest to the spectral region to be considered in this study and their references are also given in column 3 of Table 4.

For biaxial birefringent crystals such as orthorhombic ternary chalcogenides, the effective NLO coefficients for a given direction of the F-wave can generally be expressed as a linear combination of *dil* components as follows:

$$d\_{eff}^{(t)} = \mathfrak{J}\_1^{\mathfrak{x}(t)} d\_{15} + \mathfrak{J}\_2^{\mathfrak{x}(t)} d\_{24} + \mathfrak{J}\_3^{\mathfrak{x}(t)} d\_{31} + \mathfrak{J}\_4^{\mathfrak{x}(t)} d\_{32} + \mathfrak{J}\_5^{\mathfrak{x}(t)} d\_{33} \tag{20}$$

where the superscript *t* can be either I or II and represents Type I and Type II NLO interactions, respectively. The *dil* components on the crystallographic axes were chosen as the measurements closest to the spectral region to be considered in this study and the values are given in Table 5 with references. Each *ξ* coefficient in Equation (20) is determined by the relationship between the optical axes and the crystallographic axes of a biaxial crystal [48]. For *mm*2 crystals (e.g., LIS, LISe, LGS, and LGSe) showing the relationships of (*y*, *x*, *z*)=(*a*, *b*, *c*), the *ξ*-coefficients are given as follows:


**Table 4.** Effective NLO coefficients (*deff*) of chalcopyrite crystals showing uniaxial birefringence.

<sup>1</sup> The walk-off angle *ρ* is given as a function of the polar angle *θ* at resonance (see Equation (19)).

$$\begin{array}{l} \mathbb{Z}\_{1}^{(I)} = 2AH(BCH + EG)(BCE - GH),\\ \mathbb{Z}\_{2}^{(I)} = 2AH(BEG + CH)(BGH - CE),\\ \mathbb{Z}\_{3}^{(I)} = AE(BCH + EG)^{2},\\ \mathbb{Z}\_{4}^{(I)} = AE(BGH - CE)^{2},\\ \mathbb{Z}\_{5}^{(I)} = A^{3}H^{2}E,\end{array} \tag{21}$$

$$\begin{array}{lcl} \mathcal{I}\_{5}^{(II)} &=& -AE(BGH+GE)(BCE-GH)-AH(BCE-GH)^{2}, \\ \mathcal{I}\_{5}^{(II)} &=& -AE(BGE+CH)(BGH-CE)-AH(BGE+CH)^{2}, \\ & \mathcal{I}\_{5}^{(II)} &=& -AE(BGH+GE)(BCE-GH), \\ & \mathcal{I}\_{4}^{(II)} &=& -AE(BGE+CH)(BGH-CE), \\ & & \mathcal{I}\_{5}^{(II)} &=& -AE(BGE+CH)(BGH-CE), \\ & & \mathcal{I}\_{5}^{(II)} &=& -A^{3}E^{2}H, \end{array} \tag{22}$$

where the angle-dependent parameters of *A*, *B*, *C*, *G*, *E*, and *H* represent sin*θ*, cos *θ*, sin*ϕ*, cos*ϕ*, sin*δ*, and cos*δ*, respectively. The angle *δ* introduced for convenience only is defined as:

$$\tan \delta \equiv \frac{2BGC}{A^2 \cot^2 V\_z - B^2 G^2 + C^2} \tag{23}$$

where *Vz* represents the angle between the *z*-axis and the optic axis of biaxial birefringence. For biaxial crystals with a relationship of *nx* < *ny* < *nz*, such as all four chalcogenides considered in this study, *Vz* is given by:

$$\sin V\_z = \frac{n\_z \left(n\_y^2 - n\_x^2\right)^{1/2}}{n\_y \left(n\_z^2 - n\_x^2\right)^{1/2}}.\tag{24}$$

As *Vz* is a function of RI, *deff* is given as a function of the three variables *λF*, *θ*, and *ϕ*. Thus, *deff* can be obtained by substituting the set of solutions for *λF*, *θ*, and *ϕ*, satisfying the broadband SHG conditions (i.e., Table 3). The SHG efficiency is proportional to the square of *deff* for a given direction of beam propagation [42].


**Table 5.** The *dil* components of orthorhombic ternary chalcogenides considered in this study.

<sup>1</sup> *d*<sup>31</sup> = *d*<sup>15</sup> and *d*<sup>32</sup> = *d*<sup>24</sup> under Kleinman symmetry.

#### *2.3. Spatial Walk-Off*

If the wave vector of the input light is not parallel to one of the crystallographic axes, a spatial walk-off occurs between the wave vector and the Poynting vector inside the crystal, as described in Equation (19). For uniaxial crystals, *ρ* corresponds to the angle between the *o*-wave and *e*-wave [42]. Then, the largest walk-off angle (*w*) between the SH-beam and one of the two F-beams can be defined as the maximum angle between the interacting *o* and *e*-waves. As can be appreciated from the BPM conditions in Table 2, *w* is equal to *ρω* (for positive uniaxial crystals) or *ρ*2*<sup>ω</sup>* (for negative uniaxial crystals), where *ρω* and *ρ*2*<sup>ω</sup>* represent the walk-off angles of the *e*-waves with frequencies *ω* and 2*ω*, respectively. We note that for Type II, the relationship of *ρ*2*<sup>ω</sup>* > *ρω* is valid in the uniaxial crystals considered in this study [59]. As expected from Equation (19), *w* is also given as a two-variable function for *λ<sup>F</sup>* and *θ*, and thus can be obtained by substituting a set of solutions for *λ<sup>F</sup>* and *θ* that satisfy the broadband SHG conditions in Table 2.

For biaxial crystals with the point symmetry *mm*2, the walk-off angle is expressed as:

$$\tan \rho\_j^{(m)} \equiv \left\{ n^{(m)}(j\omega) \right\}^2 \left[ \left( \frac{k\_x}{\left\{ n^{(m)}(j\omega) \right\}^{-2} - a\_j} \right)^2 + \left( \frac{k\_y}{\left\{ n^{(m)}(j\omega) \right\}^{-2} - b\_j} \right)^2 + \left( \frac{k\_z}{\left\{ n^{(m)}(j\omega) \right\}^{-2} - c\_j} \right)^2 \right]^{-1/2} \tag{25}$$

where all parameters of Equation (25) are defined in Equations (13)–(16). For Type I, *w* is simply determined by the angle between the low-RI SH-beam and the high-RI F-beam, as can be appreciated from Table 3. However, for Type II, the low-RI F-beam is always placed between the high-RI F-beam and the low-RI SH-beam [60,61]. In this case, *w* is defined as the largest angle formed by the high-RI F-beam and the low-RI SH-beam:

$$\cos w = \cos \rho\_1^{(h)} \cos \rho\_2^{(l)}.\tag{26}$$

Now, *w* is also given as a function of the three variables *λF*, *θ*, and *ϕ*, which can be obtained by substituting the solution sets of *λF*, *θ*, and *ϕ*, satisfying the broadband SHG into Equation (26). The maximum deviation between the interacting beams after passing through an NLO crystal of length *L* can be expressed as:

$$
\Delta = L \tan w.\tag{27}
$$

#### **3. Simulations and Discussion**

#### *3.1. Broadband SHG in Uniaxial Chalcopyrite and Defect Chalcopyrite Crystals*

Figure 2 shows the BPM properties of negative uniaxial crystals among chalcopyrite semiconductors and defect chalcopyrite, namely AGS, AGSe, CSP, HGS, TASe, and GaSe (see Table 1). The red and blue surfaces represent *ne*(*λF*/2,*θ*) and *no*(*λF*) (for Type I) or 2*ne*(*λF*/2,*θ*) and *no*(*λF*) + *ne*(*λF*,*θ*) (for Type II), as expected from Table 2. The base of the coordinates is the plane formed by the angles *θ* and *λF*, indicating the direction of the F-wave vector and the wavelength, respectively. Each intersection of the two surfaces in Figure 2 indicates that the BPM is possible over a specific range of *λF*. The results show that both Type I and Type II NLO interactions are possible in all considered cases of AGS, AGSe, CSP, HGS, TASe, and GaSe. Figure 3 shows the BPM properties of chalcopyrite semiconductors and defect chalcopyrite exhibiting positive uniaxial birefringence (i.e., CGA, ZGP, and CdSe, as listed in Table 1). In this case, the red and blue surfaces represent *no*(*λF*/2) and *ne*(*λF*,*θ*) (for Type I) or 2*no*(*λF*/2) and *no*(*λF*) + *ne*(*λF*,*θ*) (for Type II). As shown in Figure 3d,f, the two surfaces do not intersect and thus Type II interactions are not possible for ZGP and CdSe. Therefore, for chalcopyrite crystals with positive uniaxial birefringence (see Table 1), two types of BPM are possible in CGA, whereas only Type I interactions are possible in ZGP and CdSe. However, as the effective nonlinearity of CdSe for Type I is zero (i.e., *deff* = 0 in Table 4), consequently, SHG is not possible for both Types I and II in CdSe.

**Figure 2.** BPM properties of negative uniaxial crystals among chalcopyrite semiconductors and defect chalcopyrite exhibiting negative uniaxial birefringence: (**a**) Type I and (**b**) Type II in AGS; (**c**) Type I and (**d**) Type II in AGSe; (**e**) Type I and (**f**) Type II in CSP; (**g**) Type I and (**h**) Type II in HGS; (**i**) Type I and (**j**) Type II in TASe; and (**k**) Type I and (**l**) Type II in GaSe. The red and blue surfaces represent *ne*(*λF*/2,*θ*) and *no*(*λF*) (for Type I) or 2*ne*(*λF*/2,*θ*) and *no*(*λF*) + *ne*(*λF*,*θ*) (for Type II).

**Figure 3.** BPM properties of chalcopyrite semiconductors and defect chalcopyrite exhibiting positive uniaxial birefringence: (**a**) Type I and (**b**) Type II in CGA; (**c**) Type I and (**d**) Type II in ZGP; and (**e**) Type I and (**f**) Type II in CdSe. The red and blue surfaces represent *no*(*λF*/2) and *ne*(*λF*,*θ*) (for Type I) or 2*no*(*λF*/2) and *no*(*λF*) + *ne*(*λF*,*θ*) (for Type II).

Figure 4 shows the numerical simulation results of the BPM and GV matching for the broadband SHG: Figure 4a–h correspond to AGS, AGSe, CSP, HGS, TASe, GaSe, CGA, and ZGP, respectively. The solid red and blue lines in each graph represent the BPM curves for Type I and Type II, respectively, corresponding to the intersection lines in Figures 2 and 3. The dashed magenta and cyan lines in each graph correspond to the GV matching curves for Type I and Type II, respectively, which are plotted using Equations (7)–(10). The intersection of the red and magenta (or blue and cyan) curves in Figure 4 indicates the specific direction of the F-wave vector (i.e., *θBPM*) and the *λ<sup>F</sup>* value corresponding to the SH resonance for Type I (or Type II). As discussed earlier, for ZGP, the BPM and GV matching curves intersect only for Type I (see Figure 4g). The F-wave resonances for the broadband SHG, the corresponding SH wavelengths, and the directions of the F-wave vector for the BPM obtained at the intersections in Figure 4 are summarized in Table 6. The Sellmeier equations of each crystal used in the calculations are listed together in Table 6. The F-wave resonances are interspersed in the spectral range of 4.37–9.47 μm, which corresponds to the spectral range of various mid-IR lasers such as high-power quantum cascade lasers (QCLs) based on buried-ridge or strain-balanced waveguides (WGs); distributed feedback (DFB) lasers based on plasmon-enhanced WGs or corrugated surface gratings; external cavity lasers in various bound-to-continuum designs; solid state lasers based on chalcogenide crystals doped with Fe2+; and gas lasers [62–75].

**Figure 4.** BPM and GV matching curves of (**a**) AGS, (**b**) AGSe, (**c**) CSP, (**d**) HGS, (**e**) TASe, (**f**) GaSe, (**g**) CGA, and (**h**) ZGP. The intersection of the red and magenta (or blue and cyan) curves indicates the direction of the F-wave vector (i.e., *θ*) and the *λ<sup>F</sup>* value corresponding to the SH resonance for Type I (or Type II).

° ° ° ° **Table 6.** Broadband SHG conditions in mid-IR uniaxial crystals: the resonant F-wavelengths (*λF*) and the corresponding SH wavelengths, (*λSH*); the BPM direction (*θBPM*); the effective NLO coefficients (*deff*); the maximum walk-off angles between the interacting waves (*w*); and the acceptable F-bandwidths for the broadband SHG (Δ*λF*).


<sup>1</sup> The tensor components (*dil*) used for the calculations are listed in Table 4.

Considering the NLO efficiency is proportional to the square of *deff*, it is critical to estimate the effective nonlinearity for the given direction of the F-wave vector that satisfies the broadband SHG. Figure 5 shows the polar-plots of the *deff* values as a function of *θ*, which were numerically calculated using the analytical equations for chalcopyrite uniaxial crystals in Table 4. The walk-off effects in Equation (19) were considered in the calculation results as discussed in Table 4. As shown in Figure 5, each polar-plot of *deff* exhibits two-fold or four-fold rotational symmetry with respect to *θ*, which is determined by the point symmetry of the crystal. The *deff* values calculated for the directions (*θBPM*) of the F-wave satisfying the broadband SHG conditions are listed in column 6 of Table 6. For all considered cases except CdSe, the *deff* values are very large, from 7.5 pm/V to 163.6 pm/V, which are much larger than the maximum *deff* for typical 5-mol% MgO-doped periodically poled LiNbO3 (MgO:PPLN) using the first-order QPM: (2/π)*d*<sup>31</sup> = ~2.80 pm/V for Type I and II [59]. The walk-off angle *w*, calculated using the solution set of *λ<sup>F</sup>* and *θ* for the broadband SHG, are listed in column 7 of Table 6. The calculated values of *w* are either positive (for negative uniaxial crystals) or negative (for positive uniaxial crystals) because the deviation of the *e*-ray with respect to the *o*-ray is opposite depending on the type of uniaxial crystal. The estimated walk-offs can be sufficiently overcome by a large beam window in thick crystals. The calculated values of *deff* and *w* are also summarized in Table 6 along with other conditions for the broadband SHG.

Now, we will discuss the acceptable bandwidth of the F-wave for the broadband SHG. Figure 6 shows the spectra of F-waves that are acceptable for the broadband SHG in uniaxial mid-IR crystals, in which the solid red and dashed blue lines indicate Type I and II, respectively. The spectra of F-waves were plotted using the well-known spectral equation in the coupled mode theory, in which the NLO interaction length used in the calculation is 10 mm [42]. In the proposed simultaneous BPM–GV matching scheme, the crystal length is not limited by the temporal walk-off between the interacting wave due to the zero GVM. A crystal length of 10 mm was chosen as a reference value for comparison. In the F-wave spectra of GaSe and CGA shown in Figure 6f,g, respectively, the curves for Types I and II almost overlap because the spectral positions of the resonances and their bandwidths in the two types are almost the same. In the case of ZGP, as described earlier, Type II does not exist, thus only the red curve corresponding to Type I is plotted in Figure 6h. The center wavelength (*λF*) of each graph in Figure 6 and the calculated bandwidths (Δ*λF*) are summarized Table 6. The calculated bandwidths span from 488.69 nm to 1164.24 nm in full-width-half-maximum (FWHM) as listed in Table 6. Considering such broad spectral bandwidth can cover the full spectral width of sub-picosecond pulses, as shown in Figure 6, the SHG process can be used for continuum generation, as expected in [12,30]. For example, assuming transform-limited Gaussian pulses, the calculated input bandwidths correspond to the temporal widths ranging from 56.4 fs to 113.2 fs. In particular, parametric upconversion using short pulse trains can potentially be applied to optical imaging and microscopy in the mid-IR region, in information transmission, and nonlinear optical signal processing. [27–29].

**Figure 5.** The polar-plots of the *deff* values shown as a function of *θ*, which were numerically calculated using the analytical equations for chalcopyrite uniaxial crystals in Table 4. (**a**) Type I and (**b**) Type II in AGS; (**c**) Type I and (**d**) Type II in AGSe; (**e**) Type I and (**f**) Type II in CSP; (**g**) Type I and (**h**) Type II in HGS; (**i**) Type I and (**j**) Type II in TASe; (**k**) Type I and (**l**) Type II in GaSe; (**m**) Type I and (**n**) Type II in CGA; and (**o**) Type I in ZGP.

**Figure 6.** The spectra of F-waves that are acceptable for broadband SHG in uniaxial mid-IR crystals: (**a**) AGS, (**b**) AGSe, (**c**) CSP, (**d**) HGS, (**e**) TASe, (**f**) GaSe, (**g**) CGA, and (**h**) ZGP. Solid red line: Type I; dashed blue line: Type II.

#### *3.2. Broadband SHG in Biaxial Orthorhombic Ternary Chalcogenides*

Figure 7 shows the BPM and GV matching properties of biaxial orthorhombic ternary chalcogenides, namely LIS, LISe, LGS, and LGSe in Table 1. The vertical axis of *λ<sup>F</sup>* covers the mid-IR spectral range and the base of the coordinates is the plane formed by the angles *θ* and *ϕ* representing the direction of the F-wave vector. The red and magenta (or blue and cyan) surfaces in each graph indicate the BPM and GV matching surfaces for Type I (or Type II, respectively), which were calculated using the equations in Table 3. In each graph in Figure 7, the intersection line of the two surfaces spans a specific range of *λF*, *θ*, and *ϕ*, where the BPM and GV matching is satisfied simultaneously (i.e., the broadband SHG condition). This means that the spectral position of the F-wave resonance can be selectively determined or tuned within that range of *λF*, *θ*, and *ϕ*, satisfying the broadband SHG conditions. In contrast, for uniaxial mid-IR crystals such as chalcopyrite semiconductors, the F-wave resonance can only be achieved at a single wavelength, as described in Section 3.1.

Figure 8 plots the directions of the F-wave vector (i.e., *θ* and *ϕ*) as functions of the resonant *λF*, corresponding to the intersection lines in Figure 7. As can be seen from Figure 8, for both types, the F-wave vector satisfying the broadband SHG conditions deviates more from the optical *x* and *z*-axis at shorter wavelengths. The F-wave resonance (*λF*) ranges, corresponding SH wavelength (*λSH*) ranges, and F-wave vector directions (*θ* and *ϕ*) for the broadband SHG are listed in Table 7 by the BPM type of each crystal. The Sellmeier equations for the mid-IR biaxial crystal used in the calculations are listed together in Table 7. The F-wave resonances span over the spectral range of 3.5–5.1 μm, which corresponds to the spectral range of mid-IR lasers such as high-power QCLs, DFB lasers, optical parametric oscillator lasers, solid state crystalline lasers, and gas lasers [65,66,70–74,82–85].

**Figure 7.** Numerical simulation results of BPM and GV matching in biaxial orthorhombic ternary chalcogenides: (**a**) Type I and (**b**) Type II in LIS; (**c**) Type I and (**d**) Type II in LISe; (**e**) Type I and (**f**) Type II in LGS; and (**g**) Type I and (**h**) Type II in LGSe.

**Table 7.** Broadband SHG conditions of mid-IR biaxial crystals: the ranges of *λ<sup>F</sup>* and *λSH*, and the corresponding range of the F-wave vector direction (*θ* and *ϕ*).


The magnitudes of the effective NLO coefficients (*deff*) of biaxial crystals are determined by the set of solutions for *λF*, *θ*, and *ϕ*, satisfying the broadband SHG conditions, as described in the paragraphs with Equations (20)–(24). We note, again, that it is important to estimate the *deff* values for the given conditions because the efficiency of SHG is proportional to the square of *deff*. The *deff* values were calculated numerically using Equations (20)–(24) and the results are plotted in Figure 9 as functions of the resonant *λF*, corresponding to the horizontal axis in each graph in Figure 8. This region of *λ<sup>F</sup>* corresponds to the spectral region of the F-wave resonance that satisfies the broadband SHG conditions (i.e., the intersection lines shown in Figure 7). The solid red and blue curves in Figure 9 represent the calculated *deff* values for Type I and II, respectively. For each interaction type, the change in the *deff* value with the increasing *λ<sup>F</sup>* shows a similar trend for the four kinds of biaxial crystals (LIS, LISe, LGS, and LGSe). It is interesting to note that for Type II, the *deff* value is greatest when the direction of the F-wave vector is in the *x*-*y* plane (i.e., *θ* = 90◦), as shown Figure 9 and in column 6 of Table 7. In Figure 9, the *λ<sup>F</sup>* values

obtained at the peak points of each graph are summarized in Table 8, in addition to the corresponding *λSH*, the direction of the F-wave vector (*θ* and *ϕ*), and the *deff* values. As can be seen from Table 8, the *deff* values calculated for the orthorhombic ternary chalcogenides are still larger than that of MgO:PPLN (~2.8 pm/V), as discussed in Section 3.1, but overall are slightly smaller than those of the chalcopyrite crystals (see Table 6).

**Figure 8.** The directions (*θ* and *ϕ*) of F-wave vectors plotted as functions of *λ<sup>F</sup>* satisfying the broadband SHG conditions. (**a**) Type I and (**b**) Type II in LIS; (**c**) Type I and (**d**) Type II in LISe; (**e**) Type I and (**f**) Type II in LGS; and (**g**) Type I and (**h**) Type II in LGSe.

° ° The maximum walk-off angles (*w*) between the interacting waves calculated at the maximum *deff* value points (in Figure 9) using Equation (26) are listed in column 8 of Table 8. The calculated *w* values are less than 1.25◦ for all considered cases. The corresponding maximum beam deviation (Δ) calculated using Equations (27) is 22.8 μm/mm, which can be sufficiently overcome by using a larger-size F-beam in thick crystals. Considering these values are smaller than those of other widely used NLO crystals (e.g., *w* = 3.5◦ and Δ = 61.16 μm/mm for BBO), longer NLO interaction lengths within crystals can be used for higher SHG efficiency [42]. In particular, the efficient SHG with these small spatial-walk-offs is very desirable for terahertz generation from mid-infrared two-color laser filaments [86,87].

**Figure 9.** The effective NLO coefficient (*deff*) plotted as functions of *λ<sup>F</sup>* satisfying the broadband SHG conditions. (**a**) Type I and (**b**) Type II in LIS; (**c**) Type I and (**d**) Type II in LISe; (**e**) Type I and (**f**) Type II in LGS; and (**g**) Type I and (**h**) Type II in LGSe.

**Table 8.** Broadband SHG conditions for the maximum *deff*: the resonant F-wavelength (*λF*), the direction of the F-wave vector (i.e., *θ* and *ϕ*), the effective NLO coefficient (*deff*), the walk-off angle (*w*), and the acceptable bandwidth of the F-wave (Δ*λF*).


<sup>1</sup> The tensor components (*dil*) used for the calculations are listed in Table 5.

Finally, we describe the spectral bandwidths of F-waves allowed for broadband SHG in LIS, LISe, LGS, and LGSe. Figure 10 shows the F-wave spectra acceptable for the broadband SHG, in which the solid red and dashed blue curves indicate the spectra for Type I and II, respectively. The graphs were plotted using the well-known NLO coupled mode theory, in which the NLO interaction length of 10 mm was used in the calculations [42]. Each

graph in Figure 10 is plotted for the case where *deff* has a maximum value in each crystal (i.e., the cases listed in Table 8). The spectral positions of F-wave resonances (*λF*) in each graph of Figure 10 and the calculated bandwidths (Δ*λF*) are summarized in Table 8. The calculated bandwidths span from 369.30 nm to 648.12 nm in FWHM as listed in Table 8. In contrast to uniaxial crystals, the spectral position of the resonance in biaxial crystals can be tuned by sweeping the F-wave vector within the ranges of *θ* and *ϕ* as in Figure 8, while maintaining the broad bandwidth. As in the cases of uniaxial chalcopyrite crystals, such broad spectral bandwidth can cover the full spectral width of sub-picosecond pulses, allowing the SHG process to be used for continuum generation [12,30]. In particular, parametric up-conversion using short pulse trains can potentially be applied to optical imaging and microscopy in the mid-IR region, in information transmission using optical pulse signals, and in nonlinear optical signal processing [27–29].

**Figure 10.** The spectra of F-waves that are acceptable for broadband SHG in biaxial mid-IR crystals: (**a**) LIS, (**b**) LISe, (**c**) LGS, and (**d**) LGSe. Solid red line: Type I; dashed blue line: Type II.

We note that the spectra plotted in Figures 6 and 10 show the spectral shapes of the input F-waves that allow for the generation of SH waves but not the spectral shape of the output SH-beam (or pulse). In order to accurately estimate the pulse broadening (or distortion) after passing through the crystal, the parameters for dispersion and other third-order nonlinearities at the resonances that we specified in this study (satisfying the BPM and zero GVM, simultaneously) should be first measured, but this has not been done enough yet in the previous publications. In other words, sufficient information has not yet been accumulated for pulse analysis after passing the crystal at resonances located in the spectral region where GVM is zero. We therefore believe that our study is very timely and can provide initial conditions for researchers to conduct potential experiments with light pulses in a simultaneous BPM–GV matching scheme. We hope that this study will inspire our readers and help them in their future experimental research using optical pulses.

#### **4. Conclusions**

We have theoretically and numerically investigated the broadband SHG properties of non-oxide NLO crystals including uniaxial chalcopyrite semiconductors and defect chalcopyrite (i.e., AGS, AGSe, CSP, HGS, TASe, GaSe, CGA, ZGP, and CdSe) as well as biaxial orthorhombic ternary chalcogenides (i.e., LIS, LISe, LGS, and LGSe). The BPM and GV matching properties (including the spectral position of resonance and the corresponding range of the F-vector direction), effective nonlinearities, spatial walk-offs, and acceptable F-wave spectra and bandwidths were discussed both theoretically and numerically. The effective nonlinearities were generally high in the uniaxial birefringent crystals, whereas the biaxial birefringent crystals had the advantage of exhibiting angular tunability of the spectral resonance. The calculated acceptable bandwidths of the F-wave range from 488.69 nm to 1164.24 nm (for uniaxial crystals) or from 369.30 nm to 648.12 nm (for biaxial crystals) in FWHM. Such broadband SHG using short pulse trains can potentially be applied to frequency up-conversion imaging in the mid-IR region, in information transmission using optical pulse signals, and in nonlinear optical signal processing.

**Author Contributions:** Conceptualization, K.J.L.; simulation, D.L. and K.J.L.; writing—original draft preparation, I.K. and D.L.; writing—review and editing, K.J.L.; supervision, K.J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Research Foundation of Korea (grant NRF-2019R1F1A1063937); Korea Institute of Science and Technology (grant 2E29580-19-147); and Institute for Information and Communications Technology Promotion (IITP), grant funded by the Korean government (MSIT) (number 2020-0-00947).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

