A Revisit to High Thermoelectric Performance of Single-layer MoS2

Both electron and phonon transport properties of single layer MoS2 (SLMoS2) are studied. Based on first-principles calculations, the electrical conductivity of SLMoS2 is calculated by Boltzmann equations. The thermal conductivity of SLMoS2 is calculated to be as high as 116.8 Wm−1K−1 by equilibrium molecular dynamics simulations. The predicted value of ZT is as high as 0.11 at 500 K. As the thermal conductivity could be reduced largely by phonon engineering, there should be a high possibility to enhance ZT in the SLMoS2-based materials.

where κ is the thermal conductivity, B k the Boltzmann constant, V the system volume, T the temperature, the angular bracket denoting ensemble average. i v  , i  , and i r  are the velocity vector, the site energy, and the position vector of atom respectively. The distance between atom i and atom j is denote the two-body and the three-body force, respectively.
The MD simulations are carried out utilizing LAMMPS software package S3 . The Stillinger-Weber potential, whose parameters are fitted using GULP by Jiang et al. S4 , is adopted in our simulations (shown in Table S1). An appropriate empirical potential is fundamental to obtain a reliable MD calculation result. The phonon dispersion relations, which are extracted from two different empirical potentials for SLMoS2, are compared with the results from DFT calculation, as shown in Fig. S1. The potential function used in this work provides a better reproduction of the phonon dispersion, which will assure a more reliable result on thermal conductivity.
Different from previous works on MoS2 ribbons S5 , our MoS2 model is a SLMoS2 sheet. The periodic boundary conditions are applied along the two in-plane directions. There is no inter-layer interaction due to the single layer structure. Therefore, in the simulations, the inter-layer van der Waals interactions are not taken into account. The two in-plane directions studied correspond to the armchair direction and the zigzag direction.
Generally, the temperature in MD simulation, TMD, is calculated from the kinetic energy of atoms according to the Boltzmann distribution: where E is the total kinetic energy, m the atomic mass, and N the number of particles in the system, respectively.
The velocity Verlet algorithm is employed to integrate equations of motion, and the time step is 0.5 fs. Initially, the Nose-Hoover heat reservoir is used to equilibrate the system at 300 K for 5×10 5 time steps (250ps). Then, simulations run in the NPT ensemble for 500 ps (10 6 time steps) to relax the structure. After relaxation, the converged values of both the cell size and the potential energy are obtained, which makes sure that there is no stress or strain effects. Then, the structure runs another 5 ns under NVE ensemble for relaxation. Now the system is ready for heat flux recording. The heat current vector is calculated and recorded each 2.5 fs for 8 × 10 6 time steps to obtain the autocorrelation function and thermal conductivity.
In Fig. S2, the black curve shows the normalized heat current autocorrelation function (HCACF) used in Green-Kubo formula to calculate thermal conductivity, where the side length of simulation cell is 8 units and the temperature is 300 K. Additionally, the blue curve shows the thermal conductivity which is an integration of the HFACF. The thermal conductivity converges to 93.55 Wm -1 K -1 after around 220 ps due to the decay of the HFACF. The final thermal conductivity is the mean value of twelve realizations with different initial conditions.

II. Finite size effect in simulations
When using Green Kubo formula to calculate thermal conductivity, the finite size effect could arise if the simulation cell is not sufficiently large S2 . Fig. 4 shows the thermal conductivity of SLMoS2 at room temperature with different cell sizes, from 2  2  1 to 32  32  1 unit 3 . The thermal conductivity shows a strong dependence on the size when it's smaller than 8  8  1 unit 3 . However, it changes little when the size is larger. The simulation cell in our calculations is selected as large as 32  32  1 units 3 (34.7  30.0  0.616 nm 3 ) which is large enough to overcome the finite size effect.
III. The anharmonic effects of the empirical potential function In MD simulations, the inter-atomic potential parameters are fundamental to the accuracy of the thermal conductivity calculations. The inter-atomic potential parameters herein predict the thermal expansion coefficient to be 4.85 × 10 -6 K -1 at room temperature, which is lower than the results of Huang et al S6 and C. Sevik S7 that is about 6.74 × 10 -6 K -1 from the predictions of the first principles. It means that this empirical potential function somewhere underestimates the real anharmonicity and phonon-phonon scatterings. That is, the thermal conductivity would be overestimated in our work. Moreover, an overvalued thermal conductivity results in an underestimated ZT.

IV. Electron transport properties calculation details
We compared our DFT calculation in electron and phonon dispersion curves with previous results. Fig. S3 shows The Boltzmann transport equation (BTE) is utilized to predict electronic transport properties.
The assumptions, the constant relaxation time and the rigid band approximation S8 , are used in BTE calculation S9 . These strategies in transport coefficients calculation have been verified through previous works S8, S10 . With constant relaxation time assumption, transport coefficients for electrons can be obtained by where L (α) is defined as : Kim et. al predicted the mobility in MoS2 both theoretically and experimentally S11 . It shows that the experimental results are consistent with the theoretical predictions. In their theoretical model, the impurity scatterings, acoustic phonon scatterings, and optical scatterings are all taken into consideration, which produces comprehensive results. Therefore, their values of mobility are used in our calculation. However, they provide values of mobility below 300 K only. We obtain the mobility for higher temperature according to the reciprocal relationship S12,S13 . The values of mobility for 300 K, 400 K and 500 K are fitted as 180.27 cm 2 V -1 s -1 , 117.55 cm 2 V -1 s -1 and 79.92 cm 2 V -1 s -1 , respectively. Besides, based on the results from Kaasbjerg et. al S12 , the mobility is not sensitive to carrier concentration. Therefore, it is assumed that the mobility and relaxation time are independent on the concentration in our calculation.
The relaxation time of n-type SLMoS2 for 300 K, 400 K, and 500 K are fitted as 5.1710 -14 s, 3.3710 -14 s and 2.2910 -14 s, respectively. The values are in the range of the prediction based on deformation potential theory by Fan et al. S14 .

V. Thermoelectric property
With the calculated band structures, we can obtain the transport tensor of SLMoS2. In Fig.   S5, there is little difference for Seebeck coefficients along two in-plane orthogonal directions, x and y. So, average values are taken and the SLMoS2 is treated as isotropic here.
As the Fermi level shifts from band gap to valence band, we can get the transport properties under p-type carrier concentration with rigid band approximation. Fig. S6 shows the thermoelectric properties when SLMoS2 is p-type doped. Compared with the results of ZT for ntype SLMoS2 (shown in Fig. 3), the p-type has smaller ZT values, because the n-type has bigger Seebeck coefficients.