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

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, high-sensitivity 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.

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,2,3]. Nonlinear optic (NLO) crystals—particularly with second-order 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 (LiNbO₃), lithium tantalate (LiTaO₃), and potassium titanyl phosphate (KTiOPO₄) [6,7,8,9]. A high peak-power mid-IR optical amplifier using a potassium titanyl arsenate (KTiOAsO₄) 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 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 (AgGaS₂, AGS), silver gallium selenide (AgGaSe₂, AGSe), cadmium silicon phosphide (CdSiP₂, CSP), and zinc germanium phosphide (ZnGeP₂, ZGP)), defect chalcopyrite crystals (e.g., mercury thiogallate (HgGa₂S₄, HGS) and cadmium selenide (CdSe)), and orthorhombic ternary chalcogenides (e.g., lithium thioindate (LiInS₂, LIS), lithium thiogallate (LiGaS₂, LGS), and lithium gallium selenide (LiGaSe₂, LGSe)) have been reported continuously in recent years [13,14,15,16,17,18,19,20,21,22,23,24,25,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,28,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 (CdGeAs₂, CGA), AGS, AGSe, CSP, and ZGP); (2) defect chalcopyrite showing uniaxial birefringence (e.g., thallium arsenic selenide (Tl₃AsSe₃, TASe), gallium selenide (GaSe), CdSe, and HGS); and (3) orthorhombic ternary chalcogenides showing biaxial birefringence (e.g., lithium indium selenide (LiInSe₂, 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. Non-oxide nonlinear optic (NLO) crystals considered for broadband mid-IR second-harmonic generation (SHG).

Crystals	Birefringence	Crystal System	Point Group	Transparency
AgGaS2 (AGS)	Negative uniaxial	Tetragonal	$\overline{4}2m$	0.45–13 µm [31]
AgGaSe2 (AGSe)				0.71–19 µm [32]
CdSiP2 (CSP)				0.66–6.5 µm [33]
HgGa2S4 (HGS)			$\overline{4}$	0.55–11 µm [34]
Tl3AsSe3 (TASe)		Trigonal	$3m$	1.26–17 µm [35]
GaSe		Hexagonal	$\overline{6}2m$	0.65–18 µm [36]
CdGeAs2 (CGA)	Positive uniaxial	Tetragonal	$\overline{4}2m$	2.4–18 µm [37]
ZnGeP2 (ZGP)				0.74–12 µm [36]
CdSe		Hexagonal	$6mm$	0.8–20 µm [38]
LiInS2 (LIS)	Negative biaxial	Orthorhombic	$mm2$	0.4 μm–12 µm [39]
LiInSe2 (LISe)				0.5 μm–12 µm [40]
LiGaS2 (LGS)				0.32 μm–11.6 µm [41]
LiGaSe2 (LGSe)				0.37 μm–13.2 µm [41]

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

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 n_e < n_o (for negative uniaxial crystals: AGS, AGSe, CSP, HGS, TASe, and GaSe) or n_e > n_o (for positive uniaxial crystals: CGA, ZGP, and CdSe), where n_o and n_e 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:

$$\Delta k=\left|k_e(2\omega ,\theta )-2k_o(\omega )\right|=0(\mathrm{for}\mathrm{type}\mathrm{I},\mathrm{negative}\mathrm{uniaxial}),$$

(1)

$$\Delta k=\left|k_e(2\omega ,\theta )-k_o(\omega )-k_e(\omega ,\theta )\right|=0(\mathrm{for}\mathrm{type}\mathrm{II},\mathrm{negative}\mathrm{uniaxial}),$$

(2)

$$\Delta k=\left|k_o(2\omega )-2k_e(\omega ,\theta )\right|=0(\mathrm{for}\mathrm{type}\mathrm{I},\mathrm{positive}\mathrm{uniaxial}),$$

(3)

$$\Delta k=\left|k_o(2\omega )-k_o(\omega )-k_e(\omega ,\theta )\right|=0(\mathrm{for}\mathrm{type}\mathrm{II},\mathrm{positive}\mathrm{uniaxial}),$$

(4)

where each k represents the wave number of the interacting wave and is defined as k_o(jω) = (jω/c)n_o(jω) or k_e(jω, θ) = (jω/c)n_e(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, n_e(θ), can be derived as follows using the definition in [42]:

$$n_e(j\omega ,\theta )=\frac{n_o(j\omega )n_e(j\omega )}{\sqrt{n_o^2(j\omega ){\sin }^2\theta +n_e^2(j\omega ){\cos }^2\theta }}.$$

(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_g}{c},$$

(6)

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

$$n_e^{(g)}(2\omega ,\theta )=n_o^{(g)}(\omega )(\mathrm{for}\mathrm{type}\mathrm{I},\mathrm{negative}\mathrm{uniaxial}),$$

(7)

$$2n_e^{(g)}(2\omega ,\theta )=n_o^{(g)}(\omega )+n_e^{(g)}(\omega ,\theta )(\mathrm{for}\mathrm{type}\mathrm{II},\mathrm{negative}\mathrm{uniaxial}),$$

(8)

$$n_o^{(g)}(2\omega )=n_e^{(g)}(\omega ,\theta )(\mathrm{for}\mathrm{type}\mathrm{I},\mathrm{positive}\mathrm{uniaxial}),\mathrm{and}$$

(9)

$$2n_o^{(g)}(2\omega )=n_o^{(g)}(\omega )+n_e^{(g)}(\omega ,\theta )(\mathrm{for}\mathrm{type}\mathrm{II},\mathrm{positive}\mathrm{uniaxial}),$$

(10)

where each superscript g in Equations (7)–(10) indicates the group index. Now, the broad-band 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. The BPM and GV matching conditions for broadband SHG in uniaxial crystals.

Type	Condition	Negative Uniaxial	Positive Uniaxial
Type I	BPM	$n_e(2\omega ,\theta )=n_o(\omega )$	$n_o(2\omega )=n_e(\omega ,\theta )$
	GV matching	$n_e^{(g)}(2\omega ,\theta )=n_o^{(g)}(\omega )$	$n_o^{(g)}(2\omega )=n_e^{(g)}(\omega ,\theta )$
Type II	BPM	$2n_e(2\omega ,\theta )=n_o(\omega )+n_e(\omega ,\theta )$	$2n_o(2\omega )=n_o(\omega )+n_e(\omega ,\theta )$
	GV matching	$2n_e^{(g)}(2\omega ,\theta )=n_o^{(g)}(\omega )+n_e^{(g)}(\omega ,\theta )$	$2n_o^{(g)}(2\omega )=n_o^{(g)}(\omega )+n_e^{(g)}(\omega ,\theta )$

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 λ_F and θ, which determine the center wavelength of broadband SHG and the propagation direction of the F-beam.

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 mm2 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,44,45]. In this assignment, the order of RI magnitudes is given by n_x < n_y < n_z. Then, the collinear BPM conditions can be expressed as:

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

(11)

$$\Delta k=\left|k^{(l)}(2\omega )-k^{(h)}(\omega )-k^{(l)}(\omega )\right|=0(\mathrm{for}\mathrm{type}\mathrm{II},\mathrm{biaxial}),$$

(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}}}.$$

(13)

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

$$B_j=-(b_j+c_j)k_x^2-(a_j+c_j)k_y^2-(b_j+a_j)k_z^2,$$

(14)

$$C_j=b_j{c}_j{k}_x^2+a_j{c}_j{k}_y^2+b_j{a}_j{k}_z^2,\mathrm{and}$$

(15)

$$a_j=n_x(j\omega ),b_j=n_y(j\omega ),c_j=n_z(j\omega ),$$

(16)

where each k_i 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 Δn_g = 0 obtained from Equation (6) are simplified as:

$$n_g^{(l)}(2\omega )=n_g^{(h)}(\omega )(\mathrm{for}\mathrm{type}\mathrm{I},\mathrm{biaxial})\mathrm{and}$$

(17)

$$2n_g^{(l)}(2\omega )=n_g^{(h)}(\omega )+n_g^{(l)}(\omega )(\mathrm{for}\mathrm{type}\mathrm{II},\mathrm{biaxial}).$$

(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. The BPM and GV matching conditions for broadband SHG in biaxial crystals.

Type	Condition	Negative Biaxial
Type I	BPM	$n^{(l)}(2\omega )=n^{(h)}(\omega )$
	GV matching	$n_g^{(l)}(2\omega )=n_g^{(h)}(\omega )$
Type II	BPM	$2n^{(l)}(2\omega )=n^{(h)}(\omega )+n^{(l)}(\omega )$
	GV matching	$2n_g^{(l)}(2\omega )=n_g^{(h)}(\omega )+n_g^{(l)}(\omega )$

. 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 λ_F 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.

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 $\overline{4}2m,\overline{4},3m,\overline{6}2m,$ and $6mm$ (see Table 1). The corresponding analytical expressions of the effective NLO coefficients (d_eff) are listed in column 5 of

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

Crystals	Point Group	dil Component (pm/V)	BPM Type	deff Expression 1
AGS	$\overline{4}2m$	d36 = 13.7 (@ λF = 2.53 µm [50])	Type I	$-d_{36}\sin (\theta +\rho )\sin 2\phi$
			Type II	$d_{36}\sin [2(\theta +\rho )]\cos 2\phi$
AGSe	$\overline{4}2m$	d36 = 33 (@ λF = 2.3 µm [51])	Type I	$-d_{36}\sin (\theta +\rho )\sin 2\phi$
			Type II	$d_{36}\sin [2(\theta +\rho )]\cos 2\phi$
CSP	$\overline{4}2m$	d36 = 84.5 (@ λF = 4.56 µm [52])	Type I	$-d_{36}\sin (\theta +\rho )\sin 2\phi$
			Type II	$d_{36}\sin [2(\theta +\rho )]\cos 2\phi$
HGS	$\overline{4}$	d36 = 31.5, d31 = 10.5 (@ λF = 1.064 µm [53])	Type I	$\sin (\theta +\rho )(d_{36}\sin 2\phi +d_{31}\cos 2\phi )$
			Type II	$\sin [2(\theta +\rho )](d_{36}\cos 2\phi -d_{31}\sin 2\phi )$
TASe	$3m$	d22 = 32, d31 = 20 (@ λF = 10.6 µm [54])	Type I	$d_{15}\sin (\theta +\rho )-d_{22}\cos (\theta +\rho )\sin 3\phi$
			Type II	$d_{22}{\cos }^2(\theta +\rho )\cos 3\phi$
GaSe	$\overline{6}2m$	d22 = 54 (@ λF = 10.6 µm [51])	Type I	$-d_{22}\cos (\theta +\rho )\sin 3\phi$
			Type II	$d_{22}{\cos }^2(\theta +\rho )\cos 3\phi$
CGA	$\overline{4}2m$	d36 = 186 (@ λF = 5.2955 µm [55])	Type I	$d_{36}\sin [2(\theta +\rho )]\cos 2\phi$
			Type II	$-d_{36}\sin (\theta +\rho )\sin 2\phi$
ZGP	$\overline{4}2m$	d36 = 70 (@ λF = 5.2955 µm [56])	Type I	$d_{36}\sin [2(\theta +\rho )]\cos 2\phi$
			Type II	$-d_{36}\sin (\theta +\rho )\sin 2\phi$
CdSe	$6mm$	d31 = 18 (@ λF = 10.6 µm [51])	Type I	0
			Type II	$d_{31}\sin (\theta +\rho )$

¹ The walk-off angle ρ is given as a function of the polar angle θ at resonance (see Equation (19)).

. In each d_eff 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_e^2(\omega )}-\frac{1}{n_o^2(\omega )}\right)n_e^2(\omega ,\theta ),$$

(19)

where n_e(ω,θ) is given in Equation (5). To derive this expression, we used the definition in [49]. Note that in each d_eff 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 d_eff 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 d_eff. As can be appreciated from Equation (19), ρ is also a function of λ_F and θ. Therefore, for uniaxial crystals, d_eff can be obtained by substituting the solution sets of λ_F and θ that satisfy broadband SHG conditions (i.e., Table 2). The NLO efficiency is proportional to the square of d_eff for a given direction of beam propagation [42]. The d_il 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 d_il components as follows:

$$d_{eff}^{(t)}={\xi }_1^{(t)}{d}_{15}+{\xi }_2^{(t)}{d}_{24}+{\xi }_3^{(t)}{d}_{31}+{\xi }_4^{(t)}{d}_{32}+{\xi }_5^{(t)}{d}_{33},$$

(20)

where the superscript t can be either I or II and represents Type I and Type II NLO interactions, respectively. The d_il 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. The d_il components of orthorhombic ternary chalcogenides considered in this study.

Crystals	Point Group	dil Components (pm/V) 1
LIS	$mm2$	d31 = 7.2 ± 0.4, d32 = 5.7 ± 0.6, d33 = −16 ± 4 (@ λF = 2.53 µm [57])
LISe		d31 = 11.78 ± 5%, d32 = 8.17 ± 10%, d33 = −16 ± 25% (@ λF = 2.53 µm [40])
LGS		d31 =5.8, d32 = 5.1, d33 = −10.7 (@ λF = 2.53 µm [58])
LGSe		d31 = 9.9, d32 = 7.7, d33 = −18.2 (@ λF = 2.53 µm [58])

¹ d₃₁ = d₁₅ and d₃₂ = d₂₄ under Kleinman symmetry.

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 mm2 crystals (e.g., LIS, LISe, LGS, and LGSe) showing the relationships of (y, x, z) = (a, b, c), the ξ-coefficients are given as follows:

$$\begin{array}{l}{\xi }_1^{(I)}=2AH(BCH+EG)(BCE-GH),\\ {\xi }_2^{(I)}=2AH(BEG+CH)(BGH-CE),\\ {\xi }_3^{(I)}=AE{(BCH+EG)}^2,\\ {\xi }_4^{(I)}=AE{(BGH-CE)}^2,\\ {\xi }_5^{(I)}=A^3{H}^{2}E,\end{array}$$

(21)

$$\begin{array}{l}{\xi }_1^{(II)}=-AE(BCH+GE)(BCE-GH)-AH{(BCE-GH)}^2,\\ {\xi }_2^{(II)}=-AE(BGE+CH)(BGH-CE)-AH{(BGE+CH)}^2,\\ {\xi }_3^{(II)}=-AE(BCH+GE)(BCE-GH),\\ {\xi }_4^{(II)}=-AE(BGE+CH)(BGH-CE),\\ {\xi }_5^{(II)}=-A^3{E}^{2}H,\end{array}$$

(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},$$

(23)

where V_z represents the angle between the z-axis and the optic axis of biaxial birefringence. For biaxial crystals with a relationship of n_x < n_y < n_z, such as all four chalcogenides considered in this study, V_z is given by:

$$\sin V_z=\frac{n_z{(n_y^2-n_x^2)}^{1/2}}{n_y{(n_z^2-n_x^2)}^{1/2}}.$$

(24)

As V_z is a function of RI, d_eff is given as a function of the three variables λ_F, θ, and φ. Thus, d_eff 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 d_eff for a given direction of beam propagation [42].

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ω (for negative uniaxial crystals), where ρ_ω and ρ_2ω represent the walk-off angles of the e-waves with frequencies ω and 2ω, respectively. We note that for Type II, the relationship of ρ_2ω > ρ_ω 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 λ_F and θ, and thus can be obtained by substituting a set of solutions for λ_F and θ that satisfy the broadband SHG conditions in Table 2.

For biaxial crystals with the point symmetry mm2, 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}$$

(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)}.$$

(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.$$

(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 n_e(λ_F/2,θ) and n_o(λ_F) (for Type I) or 2n_e(λ_F/2,θ) and n_o(λ_F) + n_e(λ_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 n_o(λ_F/2) and n_e(λ_F,θ) (for Type I) or 2n_o(λ_F/2) and n_o(λ_F) + n_e(λ_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., d_eff = 0 in Table 4), consequently, SHG is not possible for both Types I and II in CdSe.

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 Figure 2 and Figure 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 λ_F 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. 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 (d_eff); the maximum walk-off angles between the interacting waves (w); and the acceptable F-bandwidths for the broadband SHG (Δλ_F).

Crystals	BPM Type	λF (nm)	λSH (nm)	θBPM (°)	|deff| (pm/V) 1	w (°)	ΔλF (nm)
AGS [76]	Type I	4386.04	2193.02	31.92	7.5	1.16	488.69
	Type II	4373.83	2186.92	48.87	13.5	1.27	484.71
AGSe [77]	Type I	6224.24	3112.12	39.93	21.5	0.70	870.19
	Type II	6271.35	3135.67	67.21	23.2	0.50	868.39
CSP [52]	Type I	5315.10	2657.55	42.36	57.5	0.53	697.04
	Type II	5270.76	2635.38	74.35	43.2	0.27	669.20
HGS [78]	Type I	4438.55	2219.28	35.16	27.2	0.94	506.49
	Type II	4412.72	2206.36	54.18	31.2	0.93	506.54
TASe [35]	Type I	9076.20	4538.10	17.95	18.6	2.12	1122.97
	Type II	9140.37	4570.18	25.64	24.7	2.79	1145.78
GaSe [79]	Type I	5625.28	2812.64	10.57	52.5	2.99	760.91
	Type II	5620.10	2810.05	15.08	48.2	4.12	760.70
CGA [55]	Type I	9474.36	4737.18	32.13	163.6	−1.33	1164.03
	Type II	9474.36	4737.18	48.44	135.9	−1.49	1164.24
ZGP [80]	Type I	6006.92	3003.46	45.75	70.0	−0.70	644.54
	Type II	-	-	-	-	-	-
CdSe [81]	Type I	6380.33	3190.17	53.04	0	−0.44	-
	Type II	-	-	-	-	-	-

¹ The tensor components (d_il) used for the calculations are listed in Table 4.

. 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 Fe²⁺; and gas lasers [62,63,64,65,66,67,68,69,70,71,72,73,74,75].

Considering the NLO efficiency is proportional to the square of d_eff, 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 d_eff 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 d_eff exhibits two-fold or four-fold rotational symmetry with respect to θ, which is determined by the point symmetry of the crystal. The d_eff 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 d_eff values are very large, from 7.5 pm/V to 163.6 pm/V, which are much larger than the maximum d_eff for typical 5-mol% MgO-doped periodically poled LiNbO₃ (MgO:PPLN) using the first-order QPM: (2/π)d₃₁ = ~2.80 pm/V for Type I and II [59]. The walk-off angle w, calculated using the solution set of λ_F 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 d_eff 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 up-conversion 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,28,29].

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 λ_F 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. Broadband SHG conditions of mid-IR biaxial crystals: the ranges of λ_F and λ_SH, and the corresponding range of the F-wave vector direction (θ and φ).

Crystals	BPM Type	λF Range (nm)	λSH Range (nm)	φ Range (°)	θ Range (°)
LIS [39]	Type I	3769.29—3879.91	1884.64—1939.95	26.10—0.00	90.00—45.57
	Type II	3809.60—3856.09	1904.80—1928.04	49.39—0.00	90.00—25.75
LISe [43]	Type I	4926.06—5095.49	2463.03—2547.75	22.31—0.00	90.00—51.91
	Type II	4981.24—5039.20	2490.62—2519.60	40.71—0.00	90.00—37.53
LGS [44]	Type I	3517.15—3601.95	1758.58—1800.98	31.69—0.00	90.00—54.48
	Type II	3535.79—3567.88	1767.90—1783.94	49.69—0.00	90.00—37.22
LGSe [45]	Type I	4259.98—4288.26	2129.99—2144.13	27.21—0.00	90.00—56.60
	Type II	4258.53—4269.03	2129.26—2134.52	43.32—0.00	90.00—42.01

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,71,72,73,74,82,83,84,85].

The magnitudes of the effective NLO coefficients (d_eff) 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 d_eff values for the given conditions because the efficiency of SHG is proportional to the square of d_eff. The d_eff 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 λ_F 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 d_eff values for Type I and II, respectively. For each interaction type, the change in the d_eff value with the increasing λ_F 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 d_eff 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 λ_F values obtained at the peak points of each graph are summarized in

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

Crystals	BPM Type	λF (nm)	λSH (nm)	φ (°)	θ (°)	|deff| (pm/V) 1	w (°)	ΔλF (nm)
LIS	Type I	3781.02	1890.51	26.58	76.84	9.4	0.86	402.05
	Type II	3809.60	1904.80	49.39	90.00	6.3	0.94	410.30
LISe	Type I	4940.40	2470.20	22.55	78.97	12.3	0.82	629.79
	Type II	4981.24	2490.62	40.71	90.00	10.2	1.01	648.12
LGS	Type I	3524.58	1762.29	30.67	80.07	7.1	1.03	369.30
	Type II	3535.79	1767.90	49.69	90.00	5.4	1.11	375.76
LGSe	Type I	4262.46	2131.23	26.41	80.65	12.2	1.06	488.26
	Type II	4258.53	2129.26	43.32	90.00	8.9	1.25	491.01

¹ The tensor components (d_il) used for the calculations are listed in Table 5.

, in addition to the corresponding λ_SH, the direction of the F-wave vector (θ and φ), and the d_eff values. As can be seen from Table 8, the d_eff 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).

The maximum walk-off angles (w) between the interacting waves calculated at the maximum d_eff 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].

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 d_eff 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,28,29].

We note that the spectra plotted in Figure 6 and Figure 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.