High-Performance Solid-State Thermionic Energy Conversion Based on 2D van der Waals Heterostructures: A First-Principles Study

Two-dimensional (2D) van der Waals heterostructures (vdWHs) have shown multiple functionalities with great potential in electronics and photovoltaics. Here, we show their potential for solid-state thermionic energy conversion and demonstrate a designing strategy towards high-performance devices. We propose two promising thermionic devices, namely, the p-type Pt-G-WSe2-G-Pt and n-type Sc-WSe2-MoSe2-WSe2-Sc. We characterize the thermionic energy conversion performance of the latter using first-principles GW calculations combined with real space Green’s function (GF) formalism. The optimal barrier height and high thermal resistance lead to an excellent performance. The proposed device is found to have a room temperature equivalent figure of merit of 1.2 which increases to 3 above 600 K. A high performance with cooling efficiency over 30% of the Carnot efficiency above 450 K is achieved. Our designing and characterization method can be used to pursue other potential thermionic devices based on vdWHs.

Most of theoretical work on thermionic energy conversion uses the Richardson's law for the thermionic cur- 2 3 is the Richardson constant, e is the electron charge, m * is the electron effective mass, τ is the averaged electron transmission denoting the fraction of electrons transmitted from the metal to the semiconductor, k B is the Boltzmann constant,  is the reduced Planck constant, T is the absolute temperature, and E b is the thermionic barrier height which is usually taken as a parameter to optimize the thermionic energy conversion efficiency 1,2,12 . A small E b of several k B T is found to be optimal for a single barrier thermionic energy converter 1 . In the case of 2D vdWHs, one can benefit from the large database of 2D semiconductors with different electron affinity 15 . In addition, due to quantum confinement effects, 2D materials have layer-dependent band alignment 16 . Therefore, one can tune E b by changing the number of layers. However, evaluating E b of the vdWH-based thermionic device is nontrivial due to the Fermi level pinning effect at the metal-2D material interface 17 . The interfacial property also affects τ and C th , both of which are interface dependent and cannot be estimated easily. We note that all the above mentioned difficulties can be remedied by the parameter-free density functional theory (DFT) based first-principles calculations. First-principles study of solid-state thermionic energy conversion, though rarely reported 10 , holds great potential in the field due to the strong predictive power of DFT.
A great challenge for DFT to calculate the thermionic transport in vdWHs is to accurately evaluate E b , which is directly related to the bandgap of the semiconductor layer. DFT usually underestimates the semiconductor bandgap due to the self-interaction error. In addition, the bandgap of 2D materials shows large renormalization due to the substrate or dielectric screening [18][19][20] , a dynamic polarization effect not captured by DFT 21 . This would overestimate the bandgap. The two effects seem to cancel each other to some extent but not completely. The accurate calculation of E b is of great importance since the thermionic current changes exponentially with E b and therefore the device performance is quite sensitive to E b . To this end, we use the GW approximation 22 to calculate the quasiparticle band structures and the band alignment of the vdWHs.
We have scanned a series of vdWHs using first-principles calculations and found two promising structures to have the potential for high performance. In particular, we characterize the thermionic energy conversion performance of one of the devices, namely, the Sc-WSe 2 -MoSe 2 -WSe 2 -Sc, as shown in Fig. 1(a). We predict an equivalent figure of merit, ZT of 3 at temperatures above 600 K for the proposed device and a high cooling performance with the cooling efficiency over 30% of the Carnot efficiency above 450 K.

Results
Design principles towards high-performance solid-state thermonics based on vdWHs. There are two crucial factors which can lead to good solid-state thermionic devices: low thermal conductance and suitable thermionic barrier height within several k B T. The former is satisfied in vdWHs due to their weak bonding. The thermal conductance can be further reduced by selecting 2D materials with heavier atomic masses. Therefore, WSe 2 , with record low cross-plane thermal conductivity 23 , is a good candidate. The second factor, the barrier height, can be tuned by either using different metal electrodes with different work functions or changing the number of WSe 2 layers to tune the electron affinity due to the quantum confinement effect. Note that, we should avoid using too thin WSe 2 layers, namely, mono-or bilayer WSe 2 , to avoid the quantum tunneling effect, which would degrade the device performance. On the other hand, too thick WSe 2 layers are also not preferred due to the requirement of the ballistic transport. Anderson's rule, which is a simple and approximate method to estimate the band alignment, does not always predict the correct barrier height 10 . The strong Fermi level pinning effect in 2D materials 17 rules out this simple method. Therefore, we evaluate the band alignment between WSe 2 and several metals in the metal-WSe 2 -metal configuration from DFT calculations. During the calculations, the in-plane lattice parameters of the metals (<111> surfaces are used.) are adapted to that of the WSe 2 layer while the crossplane lattices are relaxed using the van der Waals functional (methods section). Both the n-type and p-type barrier heights, defined as where CBM and VBM stands for the conduction band minimum and valence band maximum of the WSe 2 layers, are summarized in Table 1. As can be seen, the barrier heights with all the metals considered in the metal-WSe 2 -metal configurations are large both for electron and hole transport. The work functions of the metals span a range from 4.2 eV of Al to 6.1 eV of Pt. The ionization potential (IP) of WSe 2 is calculated to be 4.9 eV. Therefore, Anderson's rule which would predict ohmic contact for Pt, would not be true according to our calculations. Increasing the WSe 2 layer thickness would reduce the barrier height which is however still too high for 5 layers. Further increase of the WSe 2 layer thickness could reduce the barrier height. But thick barriers will violate the ballistic transport condition. Another method which disrupts the direct metal-2D material interactions by inserting a hexagonal boron nitride or graphene layer in between, was found to effectively reduce the barrier height 24,25 . We check this method by inserting a single graphene layer between the Au or Pt and WSe 2 layers and the corresponding configuration is Au/Pt-G-WSe 2 -G-Au/Pt. Another benefit for this configuration is that introducing graphene could lead to clean interfaces while reducing the thermal conduction further due to the phonon interface scattering. The barrier heights for both Au and Pt covered by graphene show significant reduction compared to those of the pure metal. In particular, the p-type barrier height is found to be only 0.04 eV and 0.02 eV for 3 and 5 layers of WSe 2 , respectively, sandwiched within the Pt/G electrodes. The barrier height is on the order of k B T, so Pt-G-WSe 2 -G-Pt is a promising p-type candidate for high-performance solid-state thermionic devices.
For the n-type device, first we evaluate the configuration of Sc-WSe 2 -Sc, since Sc has a very low work function of 3.5 eV. The calculated n-type barrier height for the 4 layer WSe 2 is 0.17 eV which is too high for thermionic energy conversion. Inserting a graphene layer between the Sc and WSe 2 didn't show any reduction of the barrier height. A usual way to reduce the Schottky barrier height for conventional metal-semiconductor junctions is through doping. Instead of applying ionic substitutional doping, we find the 2D transition metal dichalcogenide (TMD) family have a staggered band alignment 26 , so charge transfer doping could be obtained by stacking two or more of the TMD layers. MoSe 2 is used for this purpose since the lattice matches that of WSe 2 , which is favored for the calculations due to the small in-plane supercell. In particular, we evaluate the configuration of Sc-WSe 2 -MoSe 2 -WSe 2 -Sc. For 4 layer MoSe 2 , the calculated barrier height is 0.03 eV which is as promising as that of the p-type Pt-G-WSe 2 -G-Pt. The above discussions are based on the DFT results. In what follows, we show the more accurate GW calculations and the thermionic energy conversion performance of one of the two promising configurations, namely, the Sc-WSe 2 -MoSe 2 -WSe 2 -Sc, due to its smaller supercell size.
Quasiparticle band structure. The DFT and GW band structures of the scattering region of the proposed Sc-WSe 2 -MoSe 2 -WSe 2 -Sc device are shown in Fig. 1(b) and (c) respectively. The fatbands of the WSe 2 layers (red lines) are significantly distorted from the ideal isolated bands 27 , indicating strong hybridization of the wave functions of WSe 2 and Sc. This is the result of the short distance of 2.0 Å between Sc and the contacting Se. The band structure of the MoSe 2 layers (blue lines) are well reproduced 27 , resembling the van der Waals bonding nature. The indirect bandgap from Γ to Q, of the central quadlayer MoSe 2 was calculated to be 0.90 eV and 1.11 eV and the values for the direct bandgap at K are 1.38 eV and 1.65 eV using DFT and GW calculations, respectively. Normally, GW predicts much larger bandgap than DFT for 2D materials, e.g., around 1 eV, 0.7 eV, and 0.5 eV larger for monolayer, bilayer, and bulk MoS 2 , respectively [27][28][29][30] . Such small correction for the present structure, i.e., 0.2-0.3 eV, reminds us of the significant bandgap renormalization due to the substrate screening effect 19,20 . Note that the renormalization is substrate dependent, e.g., the bandgap of MoS 2 is reduced by almost 1 eV on Au substrate 20 , while the reduction is only 0.13 eV for MoSe 2 on bilayer graphene 19 . In the proposed device configuration, the dynamic screening effect is more effective and significant 31 , due to both higher screening of the sandwiched structure compared to MoS 2 on a substrate (single sided) and the presence of metallic electrodes with substantially higher density of electrons compared to the graphene electrodes.
The Fermi level E F is located near the conduction bands of MoSe 2 , as shown in Fig. 1(b) and (c), which means MoSe 2 is n-type doped. We find that the GW correction on the transport barrier height is valley dependent. For  Table 1. PBE barrier heights of the metal-WSe 2 -metal configurations. G denotes graphene which covers the metal surface, so the corresponding configuration is metal-G-WSe 2 -G-metal. All the numbers are in eV.
the K valley, the barrier height is changed from 0.20 eV of DFT to 0.33 eV of GW. However, for the Q valley, the correction is only 0.03 eV, i.e., from 0.03 eV to 0.06 eV. Since the large energy difference between the K and Q valleys compared to k T B , the transport property is dominated by the Q valley. Thermionic transport properties. Figure 2(a) displays the electron transmission function τ el of the proposed vdWH device. There is a sharp slope of τ el with respect to energy at 0.06 eV above E F . This slope is due to the Q valley states of the MoSe 2 layers, consistent with the band structures. In the case of thermoelectric materials, it is known that sharp slopes of the differential conductivity with respect to energy results in asymmetry in electron-hole transport and therefore enhances the Seebeck coefficient. Similarly, sharp slopes of τ el curve favors the thermionic energy conversion by increasing the asymmetry between low and high energy electrons or in other words by acting as an effective barrier to filter out the low energy electrons. On the hole side, τ el arises at 1.3 eV below E F . From Fig. 1(c), we know that the valence band maximum (VBM) of the MoSe 2 layers is located at 1 eV and 1.3 eV below E F at Γ and K, respectively. Therefore, the holes at Γ do not contribute to the transmission, which can be understood from the small LDOS shown in Fig. 1(d). The τ el is the electrical conductance G at zero temperature. At higher temperatures, electrons have higher kinetic energies and more electrons in Sc can be emitted to the MoSe 2 conduction bands overcoming the barrier E b , contributing to the increase of G, as shown in Fig. 2(c). The n-type doping of the vdWH can be further verified by the negative sign of the Seebeck coefficient S. We obtained a maximum S of 375 μVK −1 at 160 K, while the room temperature value is 320 μVK −1 .
The phonon transmission function τ ph is shown in Fig. 2(b). The cutoff frequency ω c of τ ph is determined by the smallest ω c of the components since three-phonon scattering was not considered. The ω c of the Sc, WSe 2 , and MoSe 2 layers are 7.5 THz, 8.2 THz, and 10.4 THz, respectively, as determined by the phonon projected DOS (see Figure S1), in agreement with the literature [32][33][34] . Hence, in the present structure, the highest phonon transport channels are determined by Sc. The phonon thermal conductance C ph is saturated at 16 MW m −2 K −1 above 200 K, as shown in Fig. 2(d). This ideal value is four times larger than that of similar configuration of graphene-phosphorene-graphene sandwiched by gold electrodes 10 . This is because gold has lower ω c of 4.7 THz 35 and there is larger acoustic mismatch in that structure. The electron thermal conductance C el increases following the trend of G as the temperature increased. At room temperature, phonons dominate the thermal transport, until above 500 K, when C el becomes larger. = + , to compare with thermoelectric materials. The obtained room temperature ZT is 1.2 which is competitive to the commercial thermoelectric materials with ZT around unity. The ZT increases above 3 at 600 K and reaches a maximum value of 3.4 at 800 K. The experimental record ZT is 2.6 obtained for SnSe at 923 K 36 . Hence, our proposed device is very promising for energy conversion. Since the transport barrier E b can be tuned by changing the number of MoSe 2 layers, we plot ZT at varied E b assuming the same C ph and shifted τ el , as shown in Fig. 3(b). This approximation was verified by explicitly calculating the transport properties of the WSe 2 -2MoSe 2 -WSe 2 structure sandwiched by Sc. The E b of the bilayer MoSe 2 from GW calculation is 0.11 eV, also 0.03 eV larger than that of DFT. The calculated C ph is 15 MW m −2 K −1 which is only 1 MW m −2 K −1 smaller than that of the quadlayer case (see Figure S2). The room temperature ZT is 0.5 which agrees with the plot in Fig. 3(b). For the hexalayer MoSe 2 , the E b is 0.04 eV after applying a GW correction of 0.03 eV (see Figure S3). We label these structures of Sc-WSe2-nMoSe2-WSe2-Sc with the number (n) of the MoSe2 layers in Fig. 3(b). For solid-state thermionic transport, the general constraint on the barrier width  is where λ is the electron mean free path and  t is the minimum thickness to prevent the electron from tunneling through the barrier 1 . For the present structure with 6 layers of MoSe 2 , we have  = 4.8 nm, smaller than typical λ of 5-10 nm for most semiconductors. For the bilayer MoSe 2 with  = 2.4 nm, however, we found electron tunneling (see Figure S4) will degrade the device performance. As shown by the dashed line in Fig. 3(b), the optimum barrier height for the maximum ZT increases as the temperature increased. Label X is for Sc-4MoSe2-Sc which has a E b of 0.2 eV and a C ph of 19 MW m −2 K −1 (see Figures S5 and  S2). This structure has a large ZT at very high temperature but nearly zero at room temperature. Note that our calculations are based on ballistic transport regime which may not be valid at very high temperatures.
Thermionic refrigeration. We also propose the present vdWH as a thermionic refrigerator. The working principle is shown in Fig. 4(a). A bias V is applied with the chemical potentials µ µ > R L to pump heat J Q from right to left. As a result of this heat flow a temperature difference ΔT is established between left and right, causing a heat back flow of C ph ΔT. The cooling coefficient of performance (COP) is defined as the heat coming out of the cold side divided by the input electrical work: where J is the electrical current. We have shown that COP can be evaluated from first-principles 10 . We define the cooling efficiency normalized by the Carnot COP for cooling ( The optimum efficiency η max obtained for the optimum E b at different working temperatures is shown in Fig. 4(b). η max follows the same trend as ZT shown in Fig. 3(a). Over 30% of the Carnot efficiency is achieved with the temperature above 450 K. The efficiency changes slightly as ΔT increases. At room temperature, η max increases from 19.45% to 19.89% as ΔT increases from 5 K to 40 K.
For thermoelectric materials, the relation between ZT and the efficiency of thermoelectric cooling 1 ( 1)] so that the cooling efficiency with respect to Carnot becomes η . Take the room temperature ZT of 1.2 of the present structure and ΔT of 1 K, we obtain a η TE of 19.3%, almost equal to the thermionic η TI of 19.4% obtained from the definition of COP under the same condition. We use this argument to show that the "ballistic" ZT of the thermionic device, which includes the contacts to leads, using the standard definition of the efficiency 1 leads to almost the same efficiency directly derived from the COP introduced a few lines above. In this sense, the "ballistic" ZT leads to a consistent comparison with thermoelectrics ZT. We also compared our first-principles results with those obtained based on Richardson's law 1,2 taking the calculated E b of 0.06 eV and C ph of 16 MW m −2 K −1 as input. We assumed the averaged electron transmission τ = 1 and the electron effective mass m * = 1 in the Richardson constant A. As can be seen from the Fig. 4(b), the overall agreement is satisfactory. We note that this agreement is fortuitous since the perfect electron transmission is assumed in the Richardson equation. In reality, τ should be smaller than 1 due to the back-scattering of the electron waves at the interface and the existence of the barrier; therefore the effective mass must be greater than 1 to reproduce the first-principles results. We note that the device could be further optimized by engineering the phonon thermal conductance, e.g., by introducing lattice mismatch or disorder, since the room temperature thermal conductance is dominated by phonons as shown in Fig. 3(d).
In summary, we have shown how simple design principles and use of accurate first principles calculations, could lead to high-performance for solid-state vdWH thermionic energy conversion devices. we proposed two promising devices, namely, the p-type Pt-G-WSe 2 -G-Pt and n-type Sc-WSe 2 -MoSe 2 -WSe 2 -Sc. The performance of the latter is characterized by more accurate GW and real space Green' function calculations. We find even though the GW correction due to the large dynamic screening effect from the Sc electrodes to the barrier is small, in absolute value (0.03 eV), it still is on the order of the energy barrier itself, and thus strongly affects the thermionic current. The proposed structure has a room temperature ZT of 1.2 which increases to 3 above 600 K. A high performance with cooling efficiency over 30% of the Carnot efficiency above 450 K is predicted. Our findings show that vdWHs with appropriate electrodes have great potential when used in thermionic energy conversion devices. norm-conserving pseudopotentials 38 , PBE 39 exchange correlation functional, kinetic energy cutoff of 60 Ry, and a k mesh of 12 × 12 × 1 were used. Lattice was relaxed with the force on each atom less than 1.0e-4 Ry/Bohr. We first optimized separately the lattices of Sc, WSe 2 , and MoSe 2 . The optimized in-plane lattice constants are 3.319 Å, 3.320 Å, and 3.320 Å, respectively. We then constructed the vdWH including the electrodes with AB stacking by fixing the in-plane lattice constant of 3.319 Å across all the layers. In this configuration, there are almost no in-plane strains on any layer. With the in-plane lattice constant fixed, we further optimized the whole structure using optB86 functional 40 to include the vdW correlations. The optimized interlayer distances of MoSe 2 -MoSe 2 and MoSe 2 -WSe 2 are 6.5 Å. The relaxed coordinates can be found in supplementary files. DFT band structure calculations were performed on the relaxed structure but with PBE functional. Spin-orbit coupling (SOC) was not included since the effect is less important for Mo than W, given that the thermionic barrier is determined by the MoSe 2 layer. Moreover, SOC is less important for CBM than for VBM. Even for the CBM, the Q valley is less affected than the K valley (see Figure S6). Fatbands were calculated using the maximally localized Wannier functions (MLWFs) as implemented in the wannier90 code 41 . Initial projections for Sc, Se, W, and Mo atoms were chosen as (d; sp 3 −1), p, d, d, respectively. GW calculation. Single-shot GW or G 0 W 0 calculations based on PBE wave functions were performed using the ABINIT code 42,43 . SG15 norm-conserving pseudopotentials, wave function cutoff of 60 Ry, dielectric matrix cutoff  C of 10 Ry, exchange self-energy cutoff of 240 Ry, 1500 bands for both the dielectric matrix and self-energy calculations, and a k mesh of 9 × 9 × 1 were used. We included 4 layers of Sc at each side (total 8 layers) in the GW calculation. The transport barrier did not change with 8 more Sc layers added. The contour deformation method with 10 grids on the imaginary axis, 20 grids on the real axis, and the cutoff frequency of 1 Ry was used to calculate the dielectric matrix and self-energy. Further increasing the imaginary axis grid to 20 and real axis grid to 50 changed the results by less than 1%. GW calculations highly depend on the convergence of the parameters, of which the most important are k mesh,  C , and the number of bands (nbands). The latter two are interdependent. We checked the convergence of k mesh with  C = 10 Ry, nbands = 500, and plasmon-pole model used. Increasing the k mesh from 9 × 9 × 1 to 12 × 12 × 1 changed the barrier height by only 3 meV. The convergence of  C and nbands were checked with k meshes of 3 × 3 × 1 and 6 × 6 × 1 (see Figure S7). The accuracy with error within 0.01 eV can be obtained with  C = 20 Ry and nbands = 3000. However, using these values for the 9 × 9 × 1 k mesh is computationally prohibitive. We found that the convergence curves at different k meshes with the same C  are quite similar. Therefore, we first evaluated the barriers with C  = 10 Ry, nbands = 1500, and  C = 20 Ry, nbands = 3000 for k mesh of 6 × 6 × 1. The corresponding barrier heights are b 1 and b 2 . We then calculated the barrier height b 3 with  C = 10 Ry and nbands = 1500 for k mesh of 9 × 9 × 1. The final correct barrier height was estimated using b 3 + b 2 -b 1 . For our present calculations, we found that the correction b 2 -b 1 was very small, i.e., −0.03 eV. Thus we used the result from C  = 10 Ry and nbands = 1500 to plot the band structure, but the corrected bands were used in the transport calculations. The GW band structure was generated using the MLWFs.
Electron transport. We used the real space Green's function method 44,45 with localized basis Hamiltonian constructed using MLWFs to calculate the ballistic electron transport. The retarded Green's function of the scattering region reads: where I the identity matrix, ε the electron energy, iδ is a small imaginary number. Σ L R , denotes the self-energy of the left (L) or right (R) lead. H is the scattering region Hamiltonian, on which we imposed a minor scissor correction of 0.03 eV to take the GW correction into account. The electron transmission function is / . Note that the above method is used for 1D transport. For the case of 3D transport, one samples the 2D Brillouin zone perpendicular to the transport direction with a fine k mesh. For each transverse k point, the k-dependent transmission τ ε k ( , ) el should be calculated and the total transmission would be τ ε τ ε = ∑ k w ( ) ( , ) el k el k , where w k is the k point weight. In the present calculation, we used a k grid of 150 × 150 × 1 to obtain a smooth transmission function. The electron transmission was calculated using the WanT code 46 . With τ el , one can obtain the coherent transport coefficients under the linear response approximation 47 :  ) and obtain the remaining two from the above two equations.
Phonon transport. The phonon Green's function method 48,49 was used to calculate the ballistic phonon transport. The equations to calculate the phonon transmission function τ ph are similar to those of the electron transport. One only needs to substitute ε by ω 2 and H by Φ, where Φ is the interatomic force constant matrix divided by atomic masses. The phonon thermal conductance was calculated as: where n is the Bose-Einstein distribution function. We used the finite difference method as implemented in the SIESTA code 50 to calculate the interatomic force constant matrix. A 3 × 3 × 1 supercell was constructed and the displacement was 0.01 Å for each degree of freedom. Troullier-Martins pseudopotentials 51 , double zeta plus polarization basis set, PBE exchange correlation functional, a 6 × 6 × 1 k mesh, energy shift of 50 meV, and real space grid cutoff of 300 Ry were used for the supercell calculations. With the calculated interatomic force constant matrix as input, the phonon transmission was calculated using the transport module of the WanT code. The phonon PDOS was evaluated using the Phonopy code 52 .