Because of their high mobility and thermal conductivity, and their good mechanical properties, the interest in van der Waals crystals has increased in recent decades. In addition to graphene, chalcogenide compound semiconductors with two-dimensional sublayer structures have been prepared and studied     , and the band characteristics of these compound materials make it possible for them to be applied for light-emission and detection, communication, imaging and other related functions   . Moreover, due to the special symmetry of van der Waals crystals, which arises from the sublayers, they are considered to be ideal nonlinear-optical materials, and this has attracted interest among researchers. The nonlinear optical response of materials with different thicknesses has been studied using several methods. Taking a monolayer of MoS2 as an example, a strong resonant SHG coefficient was obtained which was attributed to electronic structural changes at the layer edges  , and it has been reported that as the number of layers increases, the position and intensity of the SHG response curve shifts. Some attempts at applying van der Waals crystals as frequency modulation devices have been made and the prospect of using them, especially in communication, has been predicted      .
Other nonlinear optical effects are the intensity of the SHG response with respect to the polarization of the light incident on the crystal in a particular electric field; thus, the SHG response is considered to be a function of the incident light frequency. It is important to find the optimal nonlinear optical crystal for a particular frequency band. Although efforts have been made on investigating nonlinear optical properties, these have been limited by the high expenditure required for crystal growth and the construction of the optical system; thus, experimental investigations of a huge number of crystals in order to compare the nonlinear coefficients of each material is costly and inefficient. In the last decade, due to the development of theoretical models and computational methods, some results, such as the susceptibilities of TMDCs, have been obtained from first principles calculations. More specifically, calculations for MoSxSe2−x were reported by Hu et al. and an interpretation of the SHG disparity at 810 nm was put forward  . Though efforts have been made, the SHG properties at wavelengths above a hundred nanometers of other typical van der Waals crystals including the transition metal monochalcogenides (TMMDs), GaSe, InSe, have rarely been studied and still need to be investigated   and in this study we calculated and directly compared the susceptibilities of 4 kinds of typical transition metal chalcogenides (GaSe, InSe, MoS2, WS2). Considering the practical applications, we chose the bulk form of typical van der Waals crystals to do the calculation, which is also different from previous research .
This study theoretically predict the SHG response of crystals in the 500 to 1100 nm wavelength range by using the calculations, which are mainly based on density functional theory (DFT) and density functional perturbation theory (DFPT) with ab initio application of the 2n + 1 theorem   . GaSe, InSe, MoS2 and WS2 were chosen as representatives not only because they are typical van der Waals crystals but also because their structural similarities allow them to be divided into two groups—MX and MX2, where M represents the metal (Ga, In, Mo) and X is the chalcogen (S, Se). With our calculations, it is possible to filter out inappropriate candidates, and find the best optical crystal needed for a given frequency band.
According to the Hohenberg and Kohn theorem, all the physical properties of an electronic system can be uniquely determined by its ground-state charge density distribution . Thus, getting information about the ground-state is the very first step in calculating the susceptibility. Density functional theory states that the system energy can be written as a functional of charge density and that the ground-state energy functional has the minimum value  . The explicit form of the energy functional can be expressed in terms of Kohn Sham eigenvalues as below .
Here N is the number of electrons, −e is the electron charge, n(r) is the charge density at position r , Exc is the exchange-correlation energy, νxc is defined as and ϵn is the Kohn Sham eigenvalue corresponding to the nth state. In the formula we have divided the system energy into three parts, and the exchange-correlation term, Exc is known to be the smallest, but cannot be ignored. Since an exact expression for the exchange-correlation energy is unknowable, Kohn and Sham proposed a reasonable approximation for it. It is considered that the exchange-correlation energy for a real electron is equal to the Exc in a homogeneous electron gas with the same charge-density distribution, and the assumption is known as the local density approximation (LDA) . Although the LDA is a really successful approach in most cases, it also well known that this method underestimates the bandgaps of materials  . For this reason, scissors-type bandgap corrections are needed in our calculations and the specific values used will be discussed later.
The second-harmonic generation corresponds to the third-order derivative of the system energy, but due to the 2n + 1 theorem, the first-order derivative of the wave function is sufficient to calculate the SHG response, and the susceptibility is as follows:
where Ω is the dimensions of the unit cell and u is the periodic part of the Bloch wave function, E0 is the external electric field, and Pc is the state projection operator . The first derivative of the wavefunction can be obtained from perturbation theory.
3. Calculational Details
Norm-conserving Fritz-Haber pseudopotentials with a cutoff energy of 100Ry were used in our calculations, which was deemed to cover the states needed to calculate the SHG response. The sampling point densities of the GaSe, InSe, MoS2 and WS2 Brillouin zones were set to 10 × 10 × 10, 7 × 7 × 7, 8 × 8 × 8, 8 × 8 × 8 respectively.
Table 1. Primitive vectors and coordinates of the basis atoms of GaSe, InSe, MoS2 and WS2, in which t1, t2, t3 represent the base vectors in Cartesian coordinate system used to decribe the unit cells. a and c are lattice consistants. Xy represents the coordinate in (t1, t2, t3) of X-atom in y-sublayer.
Figure 1. Calculation steps with the methods we applied in each steps.
4.1. Band Structure Verification
To verify the accuracy of the calculated electronic structures, we compared the band diagrams obtained with ones published in previous papers. Underestimates of the bandgaps were revealed by the comparison. Here, the band structures of GaSe InSe MoS2 and WS2 are given in Figure 2(a)-(d) respectively. By comparison with results published previously by other researchers, the accuracy of the calculations was determined and the scissors corrections obtained. These are given in Table 2 and Figure 3.
4.2. Nonlinear Optical Tensor
As we know, the second-order susceptibility matrix has 33 = 27 elements. For GaSe, InSe, MoS2, and WS2, the 211 element is the one that concerns us most: most of the other elements are canceled to zero due to the symmetry of these crystals, and coaxial phase matching is non-negligible in nonlinear optical processes. The susceptibility is a complex number, in which the real part represents the light arising, and information about the energy loss is obtained from the imaginary part combined with the real part. The absolute values of the 211 element of the susceptibility matrix, which is proportional to the SHG emergent light intensity is given in Figure 4.
In Figure 4, it is can be seen that all four materials have significant SHG effects in the visible and infrared bands, and at particular wavelengths the nonlinear optical response is several orders of magnitude larger than the band. These results are consistent with previous results. Moreover, crystals with similar structures have similar nonlinear response curves, with a short-wave-directed shift as the actual bandgap increases. The phenomenon is understandable and the result is implied by the scissors corrections obtained.
In particular, the nonlinear responses of all four crystals in the green to red light region of the visible light band reach significant levels, with InSe, GaSe and MoS2 having their highest response efficiency in this region. The results predict that the nonlinear response will give rise to peaks at 500 nm, 570 nm, 660 nm and 570 nm for the four kinds of material. As the wavelength increases, the
Figure 2. Unit cell and band structures of van der Waals crystals for (a) GaSe, (b) InSe (c) MoS2 and (d) WS2, respectively.
Table 2. Bandgaps obtained by our calculation, actual size of the gaps and the scissors corrections obtained by comparing the values.
Figure 3. Scissors correction values for the bandgaps obtained by comparing actual bandgaps with the values obtained by calculation using DFT with LDA for (a) AB-structured materials like GaSe and (b) AB2-structured materials like MoS2, respectively. Besides the parameters we need for later discussion, the points for MoSe2 and the zero energy system are shown in the plots.
intensity of the SHG response decreases and is weaker for GaSe and InSe, but we can still see some small periodic peaks in the spectral distribution as shown in Figure 3. In contrast, WS2 and MoS2 are deemed to be better choices if an IR-band SHG device is needed. Specifically, we can expect to see clear SHG peaks at around 980 nm for MoS2 and 920 nm for WS2. However, in the DFT and DFPT calculations, the effects of thermal vibrations in the lattice and defects are not included, and these may give rise to discrepancies between our calculations and the results of experiments.
Figure 4. The absolute value of the spectral distribution of the 211 element in the (2nd order) susceptibility matrices calculated by DFT method of (a) GaSe, (b) InSe, (c) MoS2 and (d) WS2 in the wavelength range from 500 nm to 1100 nm, respectively.
We used first principles calculations to obtain the SHG transition efficiencies of four typical van der Waals crystals, including the TMMCs, MoS2 and WS2, and the TMDCs, GaSe and InSe. From our calculations we obtained a linear regression for the scissors correction, and similar SHG response curves were observed for the two groups of materials, which may be due to the similarities between the crystal structures and their electronic states. Our calculations can be used to predict the SHG response of several crystals in the 500 to 1100 nm wavelength range, and so enable unsuitable nonlinear optical materials to be filtered out. SHG from GaSe and InSe in the visible light region is predicted. The calculations predict that high peaks for light emerging at 500 nm (for GaSe) and 570 nm (for InSe) would be observed. Moreover, two high SHG response regions in the visible light and infrared regions for MoS2 and WS2 are predicted by the calculations. Peaks in the SHG transformation efficiency at 660 nm and 980 nm for MoS2 and at 580 nm and 920 nm for WS2 were obtained. In the future, the present calculation will be extended to other 2D materials such as graphene, to find their potentials of nonlinear optical application.
This work was partly supported by Grant-in-Aid for JSPS Research Fellow JP19J20564, Japan.
 Gouskov, A., Camassel, J. and Gouskov, L. (1982) Growth and Characterization of III-VI Layered Crystals like GaSe, GaTe, InSe, GaSe1-xTex and GaxIn1-xSe. Progress in Crystal Growth and Characterization of Materials, 5, 323.
 Shaw, J. C., Zhou, H., Chen, Y., Weiss, N. O., Liu, Y., Huang, Y. and Duan, X. (2014) Chemical Vapor Deposition Growth of Monolayer MoSe2 Nanosheets. Nano Research, 7, 511.
 Tanabe, T., Zhao, S., Sato, Y. and Oyama, Y. (2017) Effect of Adding Te to Layered GaSe Crystals to Increase the van der Waals Bonding Force. Journal of Applied Physics, 122, Article ID: 165105.
 Oyama, Y., Tanabe, T., Sato, F., Kenmochi, A., Nishizawa, J., Sasaki, T. and Suto, K. (2008) Liquid-Phase Epitaxy of GaSe and Potential Application for Wide Frequency-Tunable Coherent Terahertz-Wave Generation. Journal of Crystal Growth, 310, 1923.
 Tanabe, T., Tang, C., Sato, Y. and Oyama, Y. (2018) Direct Determination of the Interlayer van der Waals Bonding Force in 2D Indium Selenide Semiconductor Crystal. Journal of Applied Physics, 123, Article ID: 245107.
 Takasuna, S., Shiogai, J., Matsuzaka, S., Kohda, M., Oyama, Y. and Nitta, J. (2017) Effect of Optical Waveguide on Photoluminescence Polarization in Layered Material GaSe with Millimeter Scale. Physical Review B, 96, Article ID: 161303.
 Kenmochi, A., Tanabe, T., Oyama, Y., Suto, K. and Nishizawa, J. (2008) Terahertz Wave Generation from GaSe Crystals and Effects of Crystallinity. Journal of Physics and Chemistry of Solids, 69, 605.
 Coleman, J.N., Lotya, M., O’Neill, A., Bergin, S.D., King, P.J., et al. (2011) Two-Dimensional Nanosheets Produced by Liquid Exfoliation of Layered Materials. Science, 331, 568.
 Tang, C., Sato, Y., Tanabe, T. and Oyama, Y. (2018) Low Temperature Liquid Phase Growth of Crystalline InSe Grown by the Temperature Difference Method under Controlled Vapor Pressure. Journal of Crystal Growth, 495, 54.
 Li, D.Q., Ratner, M.A. and Marks, T.J. (1988) Molecular and Macromolecular Nonlinear Optical Materials. Probing Architecture/Electronic Structure/Frequency Doubling Relationships via an SCF-LCAO MECI. pi. Electron Formalism. Journal of the American Chemical Society, 110, 1707.
 Hu, L., Wei, D. and Huang, X. (2017) Second Harmonic Generation Property of Monolayer TMDCs and Its Potential Application in Producing Terahertz Radiation. The Journal of Chemical Physics, 147, Article ID: 244701.
 Wang, C.-Y. and Guo, G.-Y. (2015) Nonlinear Optical Properties of Transition-Metal Dichalcogenide MX2 (M = Mo, W; X = S, Se) Monolayers and Trilayersfrom First-Principles Calculations. The Journal of Physical Chemistry C, 119, 13268-13276.
 Cabuk, S. (2012) The Nonlinear Optical Susceptibility and Electro-Optic Tensor of Ferroelectrics: First-Principle Study. Central European Journal of Physics, 10, 239.
 Gonze, X., Amadon, B., Anglade, P.-M., Beuken, J.-M., Bottin, F., et al. (2009) ABINIT: First-Principles Approach to Material and Nanosystem Properties. Computer Physics Communications, 180, 2582.
 Gonze, X., Rignanese, G.-M., Verstraete, M., Beuken, J.-M., Pouillon, Y., Caracas, R., Jollet, F., Torrent, M., Zerah, G., Mikami, M., Ghosez, Ph., Veithen, M., Raty, J.-Y., Olevano, V., Bruneval, F., Reining, L., Godby, R., Onida, G., Hamann, D.R. and Allan, D.C. (2005) A Brief Introduction to the ABINIT Software Package. Zeitschrift für Kristallographie, 220, 558.
 Baroni, S., De Gironcoli, S., Dal Corso, A. and Giannozzi, P. (2001) Phonons and Related Crystal Properties from Density-Functional Perturbation Theory. Reviews of Modern Physics, 73, 515.
 Assadi, M.H.N. and Hanaor, D.A. (2013) Theoretical Study on Copper’s Energetics and Magnetism in TiO2 Polymorphs. Journal of Applied Physics, 113, Article ID: 233913.
 Nagel, S., Baldereschi, A. and Maschke, K. (1979) Tight-Binding Study of the Electronic States in GaSe Polytypes. Journal of Physics C: Solid State Physics, 12, 1625.
 Plucinski, L., Johnson, R., Kowalski, B., Kopalko, K., Orlowski, B., Kovalyuk, Z. and Lashkarev, G. (2003) Electronic Band Structure of GaSe (0001): Angle-Resolved Photoemission and ab Initio Theory. Physical Review B, 68, Article ID: 125304.
 Debbichi, L., Eriksson, O. and Lebègue, S. (2015) Two-Dimensional Indium Selenides Compounds: An ab Initio Study. The Journal of Physical Chemistry Letters, 6, 3098.
 Kobayashi, K. and Yamauchi, J. (1995) Electronic Structure and Scanning-Tunneling-Microscopy Image of Molybdenum Dichalcogenide Surfaces. Physical Review B, 51, Article ID: 17085.
 Albe, K. and Klein, A. (2002) Density-Functional-Theory Calculations of Electronic Band Structure of Single-Crystal and Single-Layer WS 2. Physical Review B, 66, Article ID: 073413.