The effect of non-analytical corrections on the phononic thermal transport in InX (X = S, Se, Te) monolayers

We investigate the effect of non-analytical corrections on the phonon thermal transport properties in two-dimensional indium chalcogenide compounds. The longitudinal optical (LO) and transverse optical (TO) branches in the phonon dispersion are split near the Γ-point. The lattice thermal conductivity of monolayer InS is increased by 30.2% under non-analytical corrections because of the large LO-TO splitting at Γ-point. The predicted lattice thermal conductivities with non-analytical corrections at room temperature are 57.1 W/mK, 44.4 W/mK and 33.1 W/mK for the monolayer InS, InSe and InTe, respectively. The lattice thermal conductivity can be effectively reduced by nanostructures because the representative mean free paths are found very large in these monolayers. By quantifying the relative contribution of the phonon modes to the lattice thermal conductivity, we predict that the longitudinal acoustic branch is the main contributor to the lattice thermal conductivity. Due to the low lattice thermalconductivities of these monolayers, they can be useful in the nanoscale thermoelectric devices.

phonon spectra and lattice thermal conductivity of InX monolayers would be much stronger than their bulk and meaningful. Secondly, we expect low lattice conductivities in these monolayers because of the low elastic moduli in comparison to other two-dimensional materials 27 and these monolayers contain heavy elements such as In, Te, and Se. Wickramaratne et al. have predicted excellent electronic thermoelectric properties for these monolayers 28 , but there is a lack of study on phononic thermal transport properties which motivates us to study.
Here, we present a comprehensive study on the phonon transport properties in monolayer InX by solving the phonon Boltzmann transport equation (PBTE) based on first-principles calculations. The long-wavelength dispersion of longitudinal optical branches and lattice thermal conductivites of these monolayers are strongly affected by the non-analytical corrections to the dynamical matrix. The long-wavelength dispersion of longitudinal optical branches and lattice thermal conductivities of these monolayers are strongly affected by the non-analytical corrections to the dynamical matrix. The LO-TO splitting in monolayer InS is ten times larger than monolayer MoS 2 , which strongly affects the lattice thermal conductivity of monolayer InS. The lattice thermal conductivity trend (κ κ κ > > InS I nSe I nTe ) is explained with the help of the phonon spectra and their anharmonicities. Furthermore, the contribution of each mode toward total lattice thermal conductivity is extracted. We also discuss the dependence of the lattice thermal conductivity on temperature and size.

Results and Discussions
Bulk indium chalcogenides exist in rhombohedral, tetragonal, cubic, orthorhombic, monoclinic, and hexagonal structures [29][30][31] . Here, we study only the energetically and dynamically stable hexagonal monolayer of indium chalcogenides with space group P m 6 2 (187) with four atoms in a primitive unit cell as shown in Fig. 1. The optimized lattice parameters of monolayer InS, InSe, and InTe are a = b = 3.919 Å, 4.093 Å, and 4.382 Å, respectively, and they agree well with previous reports 32 (See Table 1). The vertical distance between two chalcogen atoms (d X−X ) is used as a thickness in the lattice thermal conductivity calculation.  The phonon band structures are shown in Fig. 2. Two in-plane acoustic modes (longitudinal acoustic (LA) and transverse acoustic (TA) modes) are linear near at the Γ point, and one out-of-plane acoustic mode (flexural acoustic (ZA) mode) has a quadratic nature near the Γ-point. The quadratic nature of ZA mode is a common feature of the two-dimensional materials, and it studied very well for graphene 13 , hexagonal boron nitride 33 , silicene 18 , and monolayer MoS 2 34 . In two-dimension systems, the total energy of the system remain the same under a global rotation due to the invariant mechanics with respect to the orientation of their reference system. This enforces an additional linear constraints on usual acoustic sum rules of the force constants derived from translational symmetry which results the quadratic nature near the Γ-point. The ZA mode is critical in thermal transport because it contributes the major part of the lattice thermal conductivity in graphene 13 . The absence of the imaginary line in phonon band structures confirms the dynamical stability of these monolayers. The phonon dispersions of these monolayers look similar, and the band gap between low-frequency optical modes and high-frequency optical modes is 82.4 cm −1 , 41.3 cm −1 , and 27.8 cm −1 for the monolayer InS, InSe, and InTe, respectively.
The non-analytical corrections are applied to the dynamical matrix by calculating the dielectric constants and Born effective charges as summarized in Table 2. The corrections split longitudinal optical (LO) and transverse optical (TO) branches at the Γ-point in these monolayers as shown in Fig. 2. These LO-TO splitting are very strong and they are about ten times larger than monolayer MoS 2 . The polarization density produced by the atomic displacement (u LO a ) and the associated long-range electric fields are the responsible for the LO-TO splitting. The polarization density (P) Fourier transform can be written as 35 : where e is the electron charge, q p is the in-plane phonon momentum, and Z a is the Born effective charge tensor associated with atom a. The polarization charge density (q P q ( ) p p ) is zero for the TO branch because the direction of propagation and the polarization is orthogonal to each other and the LO branch produces an electric field. The restoring force on the atoms is increased due to the electric field, and additional energy is required for the displacement of the LO branch with respect to the TO branch. The relationship between the frequency squares of these branches can be expressed as 35 : is the screened Coulomb interaction (which is inversely proportional to the dielectric constant) and = | | e q q / q p p p is momentum unit vector. The LO-TO splitting depends on the screened Coulomb interaction and the momentum direction (q p ) along the Born effective charges. The large LO-TO splitting in the monolayer InS, compare to monolayer InSe and InTe because the dipole moment is strongly dependent on the electronegativity difference, the larger the difference in electronegativity results the larger dipole moment. When we go down the group (from S to Te), the electronegativity decreases which decreases the dipole moment and hence the Coluomb's interaction. The second reason of the large LO-TO splitting in the monolayer is the lower mass of the S.
The lattice thermal conductivities as a function temperature are plotted for the InS, InSe, and InTe monolayers in Fig. 3(a). The lattice thermal conductivities decrease in the temperature range from 100 K to 750 K, and they are fitted well with the T 1/ lκ relationship, which demonstrates that the dominant three-phonon scattering processes in this temperature range are the Umklapp process. The lattice thermal conductivities of the three monolayers are 57.09 W/mK (InS), 44.43 W/mK (InSe), and 33.05 W/mK (InTe) at room temperature. Our calculated value of lattice thermal conductivity for the monolayer InSe agrees well with the recently reported value 36 . They possess low lattice thermal conductivity, especially for monolayer InTe. The lattice thermal conductivity of InTe is lower as compared to a lot of other two-dimensional materials, such as silicene 18 , phosphorene 15  , and stanene 20 . The possible reasons of lower lattice conductivity in monolayer InTe are the small phonon band gap between the optical modes because the small gap causes stronger scattering between the optical modes phonon and heavy mass of In and Te.
We have also computed the lattice thermal conductivities without the non-analytical corrections as shown in Fig. 3(b) to estimate the effect of LO-TO splitting and we found that the lattice thermal conductivities are strongly affected by the dipole-dipole interactions. The lattice thermal conductivities without the non-analytical corrections are decreased by 23.17%, 21.20%, and 12.62% for monolayer InS, InSe, and InTe, respectively. The non-analytical corrections will shift the optical bands in a neighborhood of Γ-point, thus changing the amount of phase space available for three-phonon scattering. The phase space for these monolayers increases substantially when we remove the non-analytical corrections, which explains the increase in thermal conductivities.
The lattice thermal conductivities are changed drastically in the low-temperature range (100 K ~ 300 K), and this change is partially attributed to the low Debye temperature (Θ α D ) of the acoustic phonon modes as given in Table 3. The Debye temperature corresponds to the temperature at which a phonon mode starts to be excited, and it is defined , where α shows ZA, TA, and LA modes, and ω α m is the maximum frequency of the cor-   Table 3. Debye temperature ( D Θ α ), representative mean free path (rMFP), specific heat (C v ), room temperature lattice thermal conductivity with non-analytical corrections ( l κ ), and lattice thermal conductivity without nonanalytical corrections (κ′ l ) of the InX monolayers.

Scientific RepoRtS |
(2020) 10:1093 | https://doi.org/10.1038/s41598-020-57644-0 www.nature.com/scientificreports www.nature.com/scientificreports/ responding phonon mode. More phonon modes are activated in this temperature range, and the population of phonons is increased. Enhancement in the phonon population leads to an increase in phonon scattering rates and hence, the lattice thermal conductivity is dramatically decreased.
The contribution of the ZA, TA, LA, and optical branches to the lattice thermal conductivity at room temperature is calculated as given in Table 4. The main contributor to the lattice thermal conductivity is the LA branch because of the large LA branch phonon group velocity and long phonon lifetime. In the case of graphene, the main contributor is ZA branch where ZA contributes by 76% 14 . In these monolayers, the acoustic branches are granted by approximately 85%, and optical branches are contributed approximately 15%. Optical branches contribute more significant as compared to graphene, stanene and monolayer MoS 2 .
The phonon properties are investigated to understand the underlying phenomena of lower lattice thermal conductivity in these monolayers and the trend of lattice thermal conductivity (InS > InSe > InTe). The solution of the PBTE within single mode relaxation time approximation (SMRTA), the lattice thermal conductivity of a two-dimensional material can be written as is the specific heat. The phonon heat capacities for the monolayers InX are calculated using the relation: and the values are given in Table 3. The specific heat of the monolayer InTe is lower than those of the monolayers InS and InSe. The lower specific heat of the monolayer InTe is due to low vibrational frequency, and it is partially responsible for lower lattice thermal conductivity of InTe among monolayers InX. Phonon group velocity is an important factor that affects the lattice thermal conductivity. Phonon group velocities of monolayer InX are calculated along the Γ-M, and Γ-K directions as shown in Fig. 4 and they are determined from the slope of the phonon dispersion. Phonon group velocities of the monolayer InTe are found lower as compared to the monolayer InS and InSe and the acoustic phonon group velocities for monolayer InTe at the Γ-point are 1793 m/s and 2747 m/s for TA, and LA branches, respectively. The group velocities of the optical modes are very low as compared to the acoustic branches, and this low group velocities of the optical branches cause lower contribution to lattice thermal conductivity.
The phonon lifetimes are extracted for each phonon mode in order to get more physical insight as shown in Fig. 5. The phonon-phonon scattering rates are dominated by isotopic and boundary scattering rates in the finite sample. In the monolayer InS, ZA mode is contributed 27.40% to the lattice thermal conductivity (larger than LA, TA and optical modes) because of longer phonon lifetime. However, LA mode is contributed more considerable in the monolayer InSe and InTe due to longer phonon lifetimes and large group velocity. The optical phonon lifetimes are very short that why they contribute very little to the lattice thermal conductivity.
Grüneisen parameter measures the anharmonicity in the chemical bonding, which drives the normal and umklapp phonon-phonon scattering processes. It is calculated from the change in phonon frequency with respect to change in lattice constant, and it can be expressed as: where γ α is the Grüneisen parameter of the α branch and a 0 is the equilibrium lattice constant. The γ for each phonon branch is computed as shown in Fig. 5(a-c) to clarify the origin of low lattice thermal conductivities in these monolayers. The γ for the ZA mode is inversely proportional to the wave vector squared (1/q 2 ), It can be easily explained from the definition of Grüneisen parameter for two-dimensional system and from the quadratic phonon dispersion of the ZA mode. The γ ZA is proportional to 1/q 2 under small positive and negative strain, since the second term in γ ZA (q) will not depend on q. Similar behavior is found graphene and BN-dope graphene 39 . The γ are anomalously large for these monolayer InX, which lead to the low lattice thermal conductivity. Large Grüneisen parameters are the consequence of weak bonding in these monolayers. Strong anharmonicity leads to short phonon lifetime because phonon-phonon scattering rates also depend on anharmonicity of the material. The effect of size on the lattice thermal conductivity is significant in the nanoscale devices because when the sample size decreases from maximal phonon mean free path (MFP), the phonon-boundary scattering is increased and thus lattice thermal conductivity is decreased. To investigate the size-dependence, the cumulative lattice thermal conductivity as a function of phonon MFP is calculated as illustrated in Fig. 6 Table 4. Percentage contribution of the ZA, TA, LA, and optical phonon branches to lattice thermal conductivity at room temperature for the monolayer InS, InSe, and InTe. (2020) 10:1093 | https://doi.org/10.1038/s41598-020-57644-0 www.nature.com/scientificreports www.nature.com/scientificreports/ evaluate the representative mean free path (rMPF, L 0 ) and the fitted curves are shown in Fig. 6. The uniparametric function is given as 40 : where κ l is the cumulative lattice thermal conductivity, and l max κ is the maximal lattice thermal conductivity. The rMFP values are tabulated in Table 3, which are larger than those of phosphorene, monolayer SnS 2 and SnSe 2 , and smaller than that of stanene 20,38 . The rMFP is very important in the designing of nanostructure because the  www.nature.com/scientificreports www.nature.com/scientificreports/ phonon-boundary scattering dominates over the three-phonon scattering when the size of the sample below rMFP.

Conclusions
In conclusions, phonon thermal transport properties, and temperature-and size-dependent lattice thermal conductivities of the monolayer InX investigated by employing first-principles calculations coupled with an iterative solution of the phonon Boltzmann transport equation The lattice thermal conductivity of these monolayers decreased with increasing temperature and perfectly follows the relation κ~T 1/ l . The predicted values of the lattice thermal conductivity at room temperature are low as compared to lots of other two-dimensional materials. The low lattice thermal conductivities originated from the strong anharmonicity, low phonon group velocity, low Debye temperature, and short phonon lifetimes. The lattice thermal conductivity can be effectively reduced by nanostructuring due to the large phonon MFP. Our work proposes that these materials can be considered for thermoelectric applications.

Methods
When a material is placed in the temperature gradient (∇T), the phonon heat flux (J) is induced and it can be written as 41 : where α is the phonon branch index, V is the volume of the unit cell, N is the number of q points in the first Brillouin zone, and ω αq , v αq , and n αq are the phonon frequency, the phonon group velocity, and the phonon distribution function, respectively. In the steady state, the phonon distribution function obeys the PBTE 41,42 : The first term denotes the diffusion of the n αq due to the temperature gradient and the second term depends on the different scattering processes, such as the phonon-phonon scattering, the phonon-boundary scattering, and the phonon-isotope scattering. For a small temperature gradient, n αq changes from its equilibrium state: = + α α α n n n q q q 0 1 , where α n q 0 is the Bose-Einstein distribution function and n q 1 α is the fluctuation in the distribution function. For small temperature gradient, the n αq can be determined by linearization of the Eq. 5 41 : The term F αq is expressed as 40 : www.nature.com/scientificreports www.nature.com/scientificreports/ where τ αq is the phonon lifetime and ∆ → αq is the correction term for the iterative solution of the PBTE. The correction term q ∆ → α is zero in the case of the single-mode relaxation time approximation (SMRTA). According to Fourier's law, the heat flux is directly proportional to the heat flux where κ l is the lattice thermal conductivity tensor. It can be estimated by solving the Boltzmann transport equation and comparing Eqs. 8 to 4 as 40 : where k B is the Boltzmann constant, i and j represent the Cartesian coordinates. The α v i is the phonon group velocity of the α branch along the direction i.
The phonon lifetime is calculated using the Matthiessen's rule, which is given as: τ α is the three-phonon scattering rate, τ α 1/ q iso is the phonon-isotope scattering rate, and 1/ q b τ α is the phonon-boundary scattering rate. The three-phonon scattering rate is calculated as 40 : where α′ and α″ represent the second and third phonon branch scattering with the phonon branch α, q q q Γ ″ ″ α α α ′ ′ + and q q q Γ ″ ″ α α α ′ ′ − are the three-phonon scattering rates of the absorption process (α + α′ → α″ and the emission process (α → α′ + α″), respectively. The phonon-isotope scattering rate can be obtained using the Tamura's formula 43 where p is the specularity parameter which describes the roughness of the boundary and L is the system size. The optimized lattice parameters and interatomic force constants (IFCs) are obtained from the total energy calculations by using plane augmented wave method 45 based on density functional theory with Vienna ab initio simulation package (VASP) 46 . We use generalized gradient approximation (GGA) parameterized by the Perdew-Burke-Ernzerhof functional 47 as exchange-correlation potential with a plane wave energy cutoff of 500 eV. All atoms in the unit cell are allowed to relax until the maximum force on each atom is smaller than 10 −4 eV/Å with a k-point mesh of 25 × 25 × 1. A vacuum thickness of 25 Å is used in order to avoid interactions between the periodic images.
Harmonic force constants are determined using the finite displacement method as implemented in the Phonopy package 48 . A 6 × 6 × 1 supercell is used for the calculations of the phonon spectra, phonon group velocity, and harmonic force constants. A 5 × 5 × 1 supercell is employed for the anharmonic force constants including the fifth nearest-neighbor interaction. ShengBTE code 40 is used to calculate the lattice thermal conductivity with a q-point mesh of 120 × 120 × 1.