Anisotropic thermoelectric behavior in armchair and zigzag mono- and fewlayer MoS2 in thermoelectric generator applications

In this work, we have studied thermoelectric properties of monolayer and fewlayer MoS2 in both armchair and zigzag orientations. Density functional theory (DFT) using non-equilibrium Green’s function (NEGF) method has been implemented to calculate the transmission spectra of mono- and fewlayer MoS2 in armchair and zigzag directions. Phonon transmission spectra are calculated based on parameterization of Stillinger-Weber potential. Thermoelectric figure of merit, ZT, is calculated using these electronic and phonon transmission spectra. In general, a thermoelectric generator is composed of thermocouples made of both n-type and p-type legs. Based on our calculations, monolayer MoS2 in armchair orientation is found to have the highest ZT value for both p-type and n-type legs compared to all other armchair and zigzag structures. We have proposed a thermoelectric generator based on monolayer MoS2 in armchair orientation. Moreover, we have studied the effect of various dopant species on thermoelectric current of our proposed generator. Further, we have compared output current of our proposed generator with those of Silicon thin films. Results indicate that thermoelectric current of MoS2 armchair monolayer is several orders of magnitude higher than that of Silicon thin films.

Compared to the research progress in its electronic and mechanical characteristics [22][23][24] , thermoelectric (TE) properties of MoS 2 have not been widely studied. Thermoelectrics provide a way of converting thermal energy into electricity 25 . Thermoelectric generator is expected to play an important role in increasing demand for clean energy in future 33 . In general, a TE generator module is made of an array of thermocouples. As illustrated in Fig. 2, each thermocouple, the basic unit of a TE generator, is made of a p-type and an n-type semiconductor, named as legs, connected thermally in parallel and electrically in series. Temperature gradient across thermocouple is the driving force inducing electrical current.
The research on thermoelectric materials has been one of the major topics since 1950 s when basic science of thermoelectrics was well founded 26 . Bi 2 Te 3 and the similar alloys have played a main role in the application of thermoelectric devices. It is well-known that efficiency of thermoelectric conversion can be evaluated by a dimensionless figure of merit κ κ = /( + ) ZT GS T e p h 2 , in which G, S, k e , k ph and T are electrical conductance, Seebeck's coefficient, electronic contribution to thermal conductance, phonon contribution to thermal conductance and absolute temperature, respectively 26 . In order to have a high ZT, it is desirable to have a high electrical conductance and large Seebeck's coefficient and low thermal conductance. These parameters mainly depend on the intrinsic properties of materials and they are generally coupled with each other. Enhancement to one of them may degrade the other and the overall effect will not change. In three decades after 1950 s, only incremental progress was made due to the difficulty in fine-tuning of these parameters 27 .
Recently, new wave of research on thermoelectric field has been initiated because nanoscale structures may enhance thermoelectric efficiency. It was shown that quantum confinement of charge carriers in quantum-well super-lattices 28 , quantum-wires 29 as well as bulk samples containing nanostructured constituents 27 will enhance thermoelectric conversion. It is known that Density of States (DOS) of low-dimensional materials exhibits sharp changes around Fermi level [27][28][29] . As a result, Seebeck's coefficient, which depends on logarithmic derivative of DOS, is significantly enhanced, and hence, the ZT increases. In addition to an increase in Seebeck's coefficient, low dimensional materials are known to benefit from higher phonon scattering and consequently lower phonon thermal conductance 27 . Low phonon thermal conductivity (k ph ) has been already reported for TMDs: κ ph of MoS 2 thin films and disordered layered WSe 2 are about 0.1 W/mK to 1 W/mK 30 and 0.05 W/mK 31 , respectively. In addition, it has been reported that MoS 2 has anisotropic thermal properties 32 , which provides another degree of freedom to optimize TE conversion performance. The advantage of nano-scale structures with respect to their bulk counterparts motivates us to study thermoelectric properties of 2D MoS 2 in both armchair and zigzag orientations.
In this work, thermoelectric properties of mono-, bi-, tri-and quadlayer MoS 2 in armchair and zigzag directions have been studied for electricity generation. ZT of bulk MoS 2 has already been reported to be 0.1 at 700 K 33 . In a later study, effect of pressure on thermoelectric properties of MoS 2 has been investigated 34 . ZT has been reported to increase up to 0.65 in a wide range of pressure and temperature. Thermoelectric performance of monolayer MoS 2 has been studied and ZT is reported to reach 0.58 in room temperature 35 . In this study, ZT values up to 1.2 in armchair direction has been achieved which is higher than ZT values reported for omnidirectional MoS 2 structures. Well-established thermoelectric materials include PbTe 36,37 and Bi 2 Te 3 38 based alloys, from which higher ZT values around 2.4 have been already achieved at 900 K. However, their substitution with abundant materials is favorable due to scarcity of Te element. This study aims to investigate the possibility of forming high performance thermoelectric generator based on highly available MoS 2 . In addition to abundance of MoS 2 , mono-and fewlayer structures have benefits of possibility of forming high density thermoelectric modules, due to their nano-scale size. In this study, it is found that as the number of layers increases from monolayer to quadlayer, both transmission spectrum and phonon thermal conductance increase. In addition, It is composed of a p-type and an n-type semiconductor, known as legs. Temperature gradient across thermocouples will induce an electrical current through thermocouple based on thermoelectric phenomenon. strong electronic and thermal transport anisotropy is found between zigzag and armchair orientations. Transmission coefficient and phonon thermal conductance of zigzag orientation is higher than those of armchair with the same number of layers. Their effect on ZT has been studied in this work. In addition, TE conversion of Si thin film TE generator with the same thickness as MoS 2 armchair mono-and fewlayer TE generator has been studied by using Synopsys TCAD software. The comparison indicates that proposed MoS 2 generator exhibits superior TE conversion efficiency.

Method
The computational model used in this paper is based on self-consistent density functional theory (DFT) using non-equilibrium Green's function (NEGF) method 39,40 implemented in QuantumWise ATK software package. Prior to the calculations of carrier transport, the structure has been relaxed to optimized force and stress of 0.05 eV/Å and 0.05 eV/ Å 3 , respectively. The relaxation calculations is implemented by using Generalized Gradient Approximation (GGA) exchange correlation with a Double Zeta Polarized (DZP) basis set and a mesh cut-off energy of 75 Ha.
Top view of structures studied in this paper is illustrated in Fig. 3. They can be divided into three regions: left, right and central. Left and right regions are called electrodes, treated with semi-infinite boundary conditions. Their properties can be described by solving for the bulk material. The voltage and temperature bias are applied on electrode regions. Central region includes a repetition of each electrode region in order to screen out perturbations introduced in the scattering regions. In order to have an insight on the effect of lattice orientation and thickness on the intrinsic TE properties of MoS 2 , no perturbation is introduced in the scattering region in calculating ZT.
Central region shown in Fig. 3, should be large enough to accommodate both the voltage and temperature drop within itself. Due to computational constraints, we used 149, 299, 449 and 599 atoms supercell as central region in mono-, bi-, tri-and quadlayer structures, respectively. Using infinitesimal voltage and temperature drop, i.e. working in linear regime, makes our approximation valid. In addition, a vacuum spacing of 20 Å is added to each side of the structure super cell to suppress any interaction caused by periodic boundary condition at out-of-plane direction.
In order to calculate linear transport properties of the system, Landauer-Buttiker 41 formula is used, in which transport coefficients are calculated from Green's function. This formulism is correct in absence of inelastic scattering and phase-changing mechanisms. DFT-NEGF method is chosen since it is proven to be a fast and computationally efficient method for a systems containing a large number of atoms 42,43 . For DFT calculations, Monkhorst-Pack k-grid 44 of 1 × 1 × 100 with a density mesh cut-off of 10 Ha is used for supercell within Localized Density Approximation (LDA) 45 with DZP basis set.
Electrical current I in the linear transport regime is given by: where factor 2 counts for spin degeneracy, q is electrical charge of carrier, h is Planck's constant, ( ) is the Fermi distribution of left (right) electrode. In linear response regime, it is assumed that As a result, equation (1) will be reduced to: Electronic contribution to TE properties, which is including electrical conductance (G), Seebeck's coefficient (S) and electronic thermal conductance (κ e ), can be calculated by using the followings: where L n is expressed as: Phonon transmission spectrum is calculated based on parameterization of Stillinger-Weber 46 potential for MoS 2 47 as implemented in Quantum Wise ATK package. Phonon thermal conductance (κ ph ) can be calculated as: and equation (7) becomes: It is worth mentioning that the phonon thermal conductance calculations in this paper are performed in the absence of any phonon decaying mechanisms. Hence, the calculations set the upper limit for phonon thermal conductance of pure MoS 2 . In reality however, there would be a few mechanisms which tend to suppress phonon conduction such as rough surface, edge imperfectness, scattering centers, etc. ZT values calculated in this paper, therefore is the minimum of what actually can be achieved by these materials. TE figure of merit is calculated by using the above information:

Results and Discussion
Transmission spectrum characterizes the electrical behavior of the proposed devices. Electrical factors that affect TE figure of merit include electrical conductance (G), electronic thermal conductance (κ e ) and Seebeck's coefficient (S). These factors can be derived from transmission spectrum as described in the previous section. Transmission spectrums for monolayer and fewlayer MoS 2 in armchair and zigzag orientations are illustrated in Fig. 4. Fermi level for is located at − = E E 0 f . Further study of Fig. 4 indicates that as the number of layers increases from one to four layers, the band gap decreases from ≈ . E 1 8 g eV to ≈ . E 1 1 g eV which is in good agreement with previously reported values 4-6 . In addition, the amplitude of transmission spectrum coefficient increases as number of layers increases from one to four, indicating that each layer provides an independent channel to conduct carriers 48 . Furthermore, zigzag orientation is found to have higher transmission coefficients in comparison with armchair. It is expected to be more conductive than armchair consequently.
In semiconducting materials, phonon thermal conductance (κ ph ) is several times larger than κ e and outplays the impact of κ e on TE figure of merit. κ ph of monolayer and fewlayer for armchair and zigzag orientations vs. temperature are illustrated in Fig. 5. As shown in Fig. 5, κ ph is almost independent of temperature. It is closely a constant in a wide range of temperatures (from 200 K to 500 K). In addition, zigzag orientation shows larger κ ph than armchair as was pointed out by Jiang 32 due to the alignment of one vibrational mode in transport direction along zigzag orientation. These results also suggest that κ ph of both zigzag and armchair orientations increases as the number of layers increases. The rate of increase in κ ph is more in zigzag than in armchair orientation. Our results of κ ph for monolayer MoS 2 is in a good agreement with findings by Huang 35 .
From factors playing role in TE figure of merit, G and κ e follow the profile of transmission spectrum, i.e. as the Fermi level moves into valence or conduction bands, transmission increases, and hence, there are more carriers to be conducted both thermally and electrically. In contrast to G and κ e , it is typical for semiconductor materials that Seebeck's coefficient (S) decreases as Fermi level moves into valence and conduction bands. Therefore G and S are competing with each other and their product in the form of S G 2 , known as power factor, reaches its maximum at an optimum position of Fermi energy 27,35 . ZT values of monolayer and fewlayer MoS 2 in armchair and zigzag orientations vs. Fermi level position at four temperatures are illustrated in Fig. 6. There are two main peaks in ZT, separated by a bandgap, corresponding to valence band and conduction band. Valence band maximum (VBM) and conduction band minimum (CBM) are specified in each plot by vertical dashed lines. In this study, TE figure of merit is referred to as ZT of n-doped or ZT of p-doped as Fermi level is approaching conduction band or valence band, respectively. It is shown in Fig. 6 that for all monolayer and fewlayer structures, ZT values of n-doped are higher than those of p-doped.
As temperature increases, amplitude of ZT also increases since ZT is proportional to the temperature. In addition, rising temperature broadens Fermi distribution. This broadening will populate states in energies higher than Fermi level, which were unpopulated in lower temperatures. These newly occupied states contribute to both electrical and thermal conduction. It means that electrical conductance increases in energies which has insignificant contribution to conduction in lower temperatures, resulting in broadening of ZT peaks vs. energy. Further study of Fig. 6 shows that profile of ZT broadens as number of layers increases for both armchair and zigzag orientations. This behavior can be attributed to the broadening of transmission spectra of both armchair and zigzag orientations as number of layer increases as illustrated in Fig. 4. Despite the increase in transmission coefficients from monolayer to quadlayer, ZT values tend to decrease as number of layers increases. This may be contrary to what is expected. One may expect that higher transmission coefficients is equivalent to more conductivity and hence higher ZT value. The reason for this behavior is due to suppression of out-of-plane vibrational mode in monolayer structures. As it can be seen from Fig. 6, ZT values undergo a sharp drop as structure changes from monolayer to bilayer for both orientations. This drop in ZT value is less pronounced when structure changes from bilayer to quadlayer. In addition, Fig. 6 suggests that ZT value of p-doped structures are smaller than those of n-doped for both orientations. This characteristic can be attributed to lower growth rate in transmission modes as Fermi level moves into valence band compared to when it moves into conduction band as illustrated in Fig. 4.
Peak values of ZT for monolayer and fewlayer armchair and zigzag orientations vs. temperature are shown in Fig. 7. As it was expected from equation (9), ZT is quite linear with temperature. ZT value is monotonously decreasing as number of layers increases. ZT value is larger than unity in n-doped armchair-oriented monolayer at = T K 500 . This structure has also the highest p-doped ZT value. Therefore, for both n-type and p-type legs in thermocouple, armchair-oriented monolayer MoS 2 is the best choice among all structures studied in this paper.
As illustrated in Fig. 6, in order to take advantage of the highest ZT value possible, MoS 2 should be doped in order to shift Fermi level to the optimum energy of peak value of ZT profile. Substitutional doping of TMD samples has been observed experimentally under exposure to 80 keV electron beam irradiation 49 . Also, a first principal study of effect of this doping approach for transition metal dopants as well as non-metal dopants is reported 50 . In order to examine the TE current of MoS 2 , we have simulated a monolayer MoS 2 in armchair orientation doped with various dopant species. We followed the same substitutional approach for doping our structure. Transition metal atoms (Re, Ru and Ta) are used as the replacing dopants for Molybdenum, and non-transition metal atoms (As, Br, Cl and P) are used for Sulphur 51 . In order to screen out the perturbation caused by doping properly, only one dopant atom was inserted in central region of device. A temperature gradient has been set across the nanoribbon by fixing the temperature of right electrode to = T K 300 and changing temperature of left electrode from = T K 300 to = T K 350 (for device configuration, see Fig. 3). Results are shown in Fig. 8. TE current of monolayer armchair MoS 2 shows strong dependence on the type of dopant atom. While Arsenic does not show any effect on thermoelectric current, P and Ta showed a similar boost to current. For n-type   at ∆ = T K 50 , more than two orders of magnitude larger. Decreasing thickness of Si film makes them more resistive and TE current decreases consequently, as shown in Fig.8. Superiority of MoS 2 -based thermocouples will be more dramatic if we compare its TE performance with those of thinner Si films, especially nm 1 -thick Si films which is almost the same thickness of monolayer MoS 2 .
Thermocouples, as was mentioned in previous section, are made of both p-type and n-type semiconductors. In order to compare the performance of monolayer MoS 2 -based and Si-based thermocouples, TE current of both of these materials is illustrated in Fig. 9. For Si-based thermocouples, p-doped (B) and n-doped (As) films with thickness of = t nm 5 and with doping concentration of = × , − N cm 1 10 A D 16 3 is used. For monolayer MoS 2 TE conversion, Ru-doped and P-doped are the best n-type and p-type structures, respectively. These two structures are chosen to construct the thermocouple based on monolayer MoS 2 . Fig. 9 shows that thermocouples based on monolayer MoS 2 are far more superior to thermocouples based on Si thin films.

Conclusion
In summary, we proposed a TE generator based on monolayer armchair-oriented MoS 2 . In order to find the optimum structure for the proposed thermocouple, first-principle simulation has been performed to calculate TE figure of merit, ZT, for monolayer and fewlayer MoS 2 in armchair and zigzag orientations. Results indicate that by increasing number of layers, ZT value tends to decrease. This behavior is in contrast to the fact that fewlayer MoS 2 is more conductive than monolayer in both directions and can be explained by suppression of out-of-plane vibrational modes in monolayer structures. Among all structures studied, monolayer armchair-oriented MoS 2 is shown to have the highest ZT value as both n-type and p-type semiconducting legs. Moreover, compared to Si thin films, TE current of monolayer MoS 2 is two orders of magnitude higher.