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

: 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 ﬁngerprint 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. broadband of propagation directions, bandwidths. The results will show that the spectral bandwidths of the F-wave for the SHG to several hundreds of nm. The SH spectral range


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 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 (LiNbO 3 ), lithium tantalate (LiTaO 3 ), and potassium titanyl phosphate (KTiOPO 4 ) [6][7][8][9]. A high peak-power mid-IR optical amplifier using a potassium titanyl arsenate (KTiOAsO 4 ) 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

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. LiInSe 2 (LISe) 0.5 µm-12 µm [40] LiGaS 2 (LGS) 0.32 µm-11.6 µm [41] LiGaSe 2 (LGSe) 0.37 µm-13.2 µm [41] 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: 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: ω ωθ 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]: 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: 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: 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]: 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: 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: e (ω, θ) (for type I, positive uniaxial), and 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 λ 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: 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]: The parameters in Equation (13) are defined as follows: C j = b j c j k 2 x + a j c j k 2 y + b j a j k 2 z , and 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: g (ω) (for type I, biaxial) and (17) 2n (l) 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 λ 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. Table 3. The BPM and GV matching conditions for broadband SHG in biaxial crystals.

Type
Condition Negative Biaxial

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 42m, 4, 3m, 62m, and 6mm (see Table 1). The corresponding analytical expressions of the effective NLO coefficients (d eff ) are listed in column 5 of Table 4. 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: 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: 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 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:  1 The walk-off angle ρ is given as a function of the polar angle θ at resonance (see Equation (19)).
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: 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: 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].

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: 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: 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: 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 fo broadband SHG: Figure 4a-h correspond to AGS, AGSe, CSP, HGS, TASe, GaSe, C and ZGP, respectively. The solid red and blue lines in each graph represent the curves for Type I and Type II, respectively, corresponding to the intersection lines in ures 2 and 3. The dashed magenta and cyan lines in each graph correspond to th matching curves for Type I and Type II, respectively, which are plotted using Equa (7)- (10). The intersection of the red and magenta (or blue and cyan) curves in Fig indicates the specific direction of the F-wave vector (i.e., θBPM) and the λF value c sponding to the SH resonance for Type I (or Type II). As discussed earlier, for ZGP BPM and GV matching curves intersect only for Type I (see Figure 4g). The F-wave nances for the broadband SHG, the corresponding SH wavelengths, and the directio the F-wave vector for the BPM obtained at the intersections in Figure 4 are summa in Table 6. The Sellmeier equations of each crystal used in the calculations are liste gether in Table 6. The F-wave resonances are interspersed in the spectral range of 9.47 µm, which corresponds to the spectral range of various mid-IR lasers such as power quantum cascade lasers (QCLs) based on buried-ridge or strain-balanced w guides (WGs); distributed feedback (DFB) lasers based on plasmon-enhanced WGs o rugated surface gratings; external cavity lasers in various bound-to-continuum des solid state lasers based on chalcogenide crystals doped with Fe 2+ ; and gas lasers [62-  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 λ 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. 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 2+ ; and gas lasers [62][63][64][65][66][67][68][69][70][71][72][73][74][75].  The intersection of the red and magenta (or blue and cyan) curves indicates the direction of the F-wave vector (i.e., θ) and the λF 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).  The intersection of the red and magenta (or blue and cyan) curves indicates the direction of the F-wave vector (i.e., θ) and the λ F 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 (d eff ); the maximum walk-off angles between the interacting waves (w); and the acceptable F-bandwidths for the broadband SHG (∆λ F ). 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 3 (MgO:PPLN) using the first-order QPM: (2/π)d 31 =~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 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][28][29]. 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 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].   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 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].

Broadband SHG in Biaxial Orthorhombic Ternary Chalcogenides
(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 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][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, 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 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 λF corresponds to the spectral region of the F-wave resonance that satisfies the broadband SHG conditions 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]. 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.   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 ).  Table 5.

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

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.