Ab-initio quantum transport simulation of self-heating in single-layer 2-D materials

TL;DR

This paper uses first-principles quantum simulations to analyze self-heating effects in 2-D materials like MoS2, WS2, and black phosphorus, revealing their impact on device performance and thermal management.

Contribution

It provides the first ab-initio analysis of self-heating in 2-D materials, highlighting differences among materials and their influence on electronic device performance.

Findings

Self-heating is more significant in 2-D materials than in Si nanowires.

High current densities can severely limit 2-D device performance due to self-heating.

Black phosphorus is less affected by self-heating compared to transition metal dichalcogenides.

Abstract

Through advanced quantum mechanical simulations combining electron and phonon transport from first-principles self-heating effects are investigated in n-type transistors with a single-layer MoS2, WS2, and black phosphorus as channel materials. The selected 2-D crystals all exhibit different phonon-limited mobility values, as well as electron and phonon properties, which has a direct influence on the increase of their lattice temperature and on the power dissipated inside their channel as a function of the applied gate voltage and electrical current magnitude. This computational study reveals (i) that self-heating plays a much more important role in 2-D materials than in Si nanowires, (ii) that it could severely limit the performance of 2-D devices at high current densities, and (iii) that black phosphorus appears less sensitive to this phenomenon than transition metal dichalcogenides.

Tables

Table 1. Table 1: Selected material parameters of single-layer MoS 2 , WS 2 , and black phosphorus along the Γ Γ \Gamma - X 𝑋 X and Γ Γ \Gamma - Y 𝑌 Y symmetry lines: band gap ( E g subscript 𝐸 𝑔 E_{g} ), electron ( m e subscript 𝑚 𝑒 m_{e} ) and hole ( m h subscript 𝑚 ℎ m_{h} ) effective masses at the band edges, sound velocity of the transverse ( v s , T ​ A subscript 𝑣 𝑠 𝑇 𝐴 v_{s,TA} ) and longitudinal ( v s , L ​ A subscript 𝑣 𝑠 𝐿 𝐴 v_{s,LA} ) acoustic phonon branches, as well as optical phonon energy ( ℏ ​ ω o ​ p ​ t Planck-constant-over-2-pi subscript 𝜔 𝑜 𝑝 𝑡 \hbar\omega_{opt} ). All data were extracted from DFT vasp and DFPT dfpt calculations.

	$E_g$ (eV)	$m_e^{\ast }$	$m_h^{\ast }$	$v_{s,TA}$ (km/s)	$v_{s,LA}$ (km/s)	$\hslash {\omega }_{opt}$ (meV)
MoS2	1.67	0.46	0.56	4.7	7.1	58
WS2	1.8	0.3	0.38	3.9	6.2	53
BPΓ-X	1.61	0.18	0.15	4.6	5.0	57
BPΓ-Y	1.61	1.2	2.5	3.9	8.5	57

Equations

$$\displaystyle\left\{\begin{array}[]{l}\sum_{l}\left(\mathbf{E}\delta_{li}-\mathbf{H}_{il}(k_{z})-\mathbf{\Sigma}^{RB}_{il}(E,k_{z})-\mathbf{\Sigma}^{RS}_{il}(E,k_{z})\right)\cdot\mathbf{G}^{R}_{lj}(E,k_{z})=\delta_{ij}\\ \mathbf{G}^{\gtrless}_{ij}(E,k_{z})=\sum_{lm}\mathbf{G}^{R}_{il}(E,k_{z})\cdot\left(\mathbf{\Sigma}^{\gtrless B}_{lm}(E,k_{z})+\mathbf{\Sigma}^{\gtrless S}_{lm}(E,k_{z})\right)\cdot\mathbf{G}^{A}_{mj}(E,k_{z}).\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{l}\sum_{l}\left(\mathbf{\omega}^{2}\delta_{li}-\mathbf{\Phi}_{il}(q_{z})-\mathbf{\Pi}^{RB}_{il}(\omega,q_{z})-\mathbf{\Pi}^{RS}_{il}(\omega,q_{z})\right)\cdot\mathbf{D}^{R}_{lj}(\omega,q_{z}z)=\delta_{ij}\\ \mathbf{D}^{\gtrless}_{ij}(\omega,q_{z})=\sum_{lm}\mathbf{D}^{R}_{il}(\omega,q_{z})\cdot\left(\mathbf{\Pi}^{\gtrless B}_{lm}(\omega,q_{z})+\mathbf{\Pi}^{\gtrless S}_{lm}(\omega,q_{z})\right)\cdot\mathbf{D}^{A}_{mj}(\omega,q_{z}),\end{array}\right.$$

$${\mathbf{\Sigma }}_{nn}^{\gtrless S}(E,k_z)$$

$${\Pi }_{nn}^{\gtrless ij}(\omega ,q_z)$$

$${\Pi }_{nl}^{\gtrless ij}(\omega ,q_z)$$

$${\mathbf{\Sigma }}^R(E,k_z)$$

$${\mathbf{\Pi }}^R(\omega ,q_z)$$

$$n({\mathbf{R}}_{\mathbf{i}})$$

$$I_{d,s\to s+1}$$

$$I_{E-el,s\to s+1}$$

$$I_{E-ph,s\to s+1}$$

$$N_{ph}^{tot}(s)$$

$$=$$

$$\mathrm{LDOS}({\mathbf{R}}_i,\omega ,q_z)$$

Full text

I Introduction

Only a fraction of the existing 2-D materials has been identified so far mounet , probably less than 50 components have been mechanically exfoliated or grown via chemical vapor deposition (CVD), but it already clearly appears that the 2-D technology might have a positive impact in various domains of application, one of them being the search for next-generation logic switches beyond Si FinFETs kuhn . While graphene, the first experimentally demonstrated 2-D compound novoselov , is not suitable as transistor channel due to the absence of a bandgap, monolayers of MoS2 and other transition metal dichalcogenides (TMDs) look very promising as silicon replacement in upcoming ultra-scaled transistors kis : they exhibit excellent electrostatic properties due to their 2-D nature, relatively large band gaps ($E_{g}\sim$1.8 eV in MoS2), high $I_{ON}$/$I_{OFF}$ current ratios, and decent room temperature single-layer mobility values, e.g. 214 cm2/Vs for WS2 encapsulated between two h-BN layers iqbal . Furthermore, because the atomic composition and crystallographic phase of TMDs fully determine whether they behave as semiconductors, insulators, or metals miro , active components relying solely on stacks of such materials can be easily imagined for the future.

The “current vs. voltage” characteristics and electron/hole mobilities of devices based on single-layer TMDs and other 2-D materials have been extensively investigated in the recent past, contrary to their thermal and electro-thermal properties (self-heating, the formation of local hot spots, breakdown failures), which have remained largely unexplored up to now. It is however expected that the influence of thermally induced effects becomes significant in ultra-scaled structures due to the close proximity of electrons and phonons and their potentially increased coupling pop0 . Few studies, mainly concerned with graphene balandin ; shi ; pop1 , but also molybdenum disulfide xing ; pop2 and black phosphorus cahill have discussed these issues. Still, it can be generally said that the thermal behavior of 2-D materials is not completely understood yet and that rapid progresses, both at the theoretical and experimental levels, are required to shed light on the underlying physics and to enable the emergence of 2-D devices with improved electro-thermal functionalities, e.g. nanoscale thermoelectric generators or Peltier coolers.

This paper intends to give a theoretical insight into the electro-thermal properties of transistors with a single-layer 2-D material as channel. For that purpose, advanced simulation techniques will be utilized. On the electrical side, the state-of-the-art consists in utilizing ballistic full-band solvers implementing the Non-equilibrium Green’s Function formalism (NEGF) datta . The necessary Hamiltonian matrix relies either on semi-empirical (tight-binding) hguo ; klimeck models, it is derived from first-principles via maximally localized Wannier functions register ; fiori1 , or it is constructed at the ab-initio level qw0 . With many TMDs, such approaches induce non-physical negative differential resistance (NDR) behaviors that have never been experimentally observed at room temperature register . The NDR originates from the TMD bandstructure, which exhibits several narrow energy bands that cannot carry any current if the electrostatic potential undergoes large variations between the source and drain contacts gnani . It is therefore an artifact of the ballistic approximation. The inclusion of electron-phonon scattering helps get rid of NDR by connecting bands that would otherwise remain independent from each other szabo1 . Combining a full-band approach with electron-phonon scattering has been applied in several occasions to 2-D materials, but only to extract their electron mobility kaas ; qw ; marzari ; fiori2 , not to obtain their current characteristics as a function of the externally applied voltages.

On the thermal modeling side, it is widely recognized that the tiny dimensions of 2-D materials require a quantum mechanical approach to treat heat transport, which should be done at the phonon level. Two different methods have been successfully applied by various research groups to obtain the thermal conductivity of single-layer 2-D materials: (i) molecular dynamics (MD) simulations md1 ; md2 and (ii) the coupling of density-functional theory (DFT) with the linearized Boltzmann Transport equation (LBTE), where the former provides the bandstructure and scattering parameters for the latter lbt1 ; lbt2 ; marzari2 . The advantage of MD over DFT+LBTE is that it allows for the consideration of non-homogeneous atomic structures with spatial or material variations. Nevertheless, neither one nor the other model lends itself naturally to a self-consistent coupling with electron transport. The NEGF formalism offers more flexibility for that because it can capture the energy exchanges between the electron and phonon populations via scattering self-energies, thus automatically accounting for local lattice temperature increases asai ; pecchia ; frederiksen ; jtlu1 ; jtlu2 .

Here, by combining the electro-thermal simulation framework of Ref. prb14 , originally developed for nanowire transistors, and the ab-initio modeling platform of Ref. szabo1 , specifically dedicated to the investigation of 2-D materials, the first theoretical results of electro-thermal transport through logic devices with a 2-D monolayer semiconductor as channel will be demonstrated. To the best of our knowledge such an accurate technique had not been tested in the context of 2-D materials before, mainly because its application is normally limited to small systems composed of 100 atoms or less due to the high computational burden associated with the description from first-principles of electron and phonon bandstructures. The paper is organized as follows: after a detailed presentation of the modeling approach in Section II, simulation results will be introduced in Section III for transistors with a single-layer MoS2, WS2, and black phosphorus channel. The low-field phonon-limited mobility, I-V characteristics, lattice temperature, and dissipated power of these devices will be analyzed and compared to each other. Finally, conclusions will be drawn in Section IV.

II Simulation Approach

II.1 Electron and Phonon Transport

Transistor structures similar to the one depicted in Fig. 1 have been simulated in this work. Transport occurs along the $x$-axis, $y$ is a direction of confinement, and the out-of-plane dimension, $z$, is assumed periodic and gives rise to a $k_{z}$(electron momentum along $z$)- or $q_{z}$(phonon momentum along $z$)-dependence of all computed physical quantities. The electron and phonon properties are modeled at the ab-initio level with NEGF. Their coupling is realized via scattering self-energies that ensure current and energy conservation. In this context, the following system of equations must be solved for the electron population

[TABLE]

In Eq. (3), $\mathbf{E}$ is a diagonal matrix that contains the electron energy $E$ as entry. The $\mathbf{G}_{ij}(E,k_{z})$’s represent the electron Green’s Functions at energy $E$ and momentum $k_{z}$ between atoms $i$ and $j$ situated at position $\mathbf{R}_{i}$ and $\mathbf{R}_{j}$, respectively. They are of size $N_{orb,i}\times N_{orb,j}$, where $N_{orb,i}$ is the number of orbitals (basis components) describing atom $i$. The $\mathbf{G}_{ij}(E,k_{z})$’s can be either retarded ($R$), advanced ($A$), lesser ($<$), or greater ($>$). The same conventions apply to the self-energies $\mathbf{\Sigma}_{ij}(E,k_{z})$ whose additional index $B$ ($S$) stands for boundary (scattering). Details about the $\mathbf{\Sigma}_{ij}(E,k_{z})$’s are provided below.

The Hamiltonian matrix elements $\mathbf{H}_{il}(k_{z})$ are key ingredients in Eq. (3). To simulate 2-D materials, they are usually either expressed in a semi-empirical basis such as tight-binding hguo ; klimeck or in a localized density-functional theory (DFT) basis qw0 . Here, an intermediate scheme is used: after identifying a small unit cell that is representative for the considered 2-D crystal, a DFT calculation is performed with a standard tool such as VASP vasp . The plane-wave outputs of VASP are then converted into a set of maximally localized Wannier functions (MLWFs) with the wannier90 package wannier90 . This allows for the creation of a tight-binding-like Hamiltonian matrix that is scaled up to match the desired structure dimensions. Finally, quantum transport simulations are run with the produced Hamiltonian as input. The whole process is highlighted in Ref. szabo1 . Note that MoS2 and WS2 have been treated here within the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) pbe , whereas the hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE06) hse06 has been employed for black phosphorus. A 25$\times$1$\times$25 Monkhorst-Pack k-point grid and a 500 eV plane-wave cutoff energy have been enforced in all electronic structure calculations. Spin-orbit coupling has been neglected. The accuracy of the plane-wave $\rightarrow$ MLWF conversion is demonstrated in Fig. 2(a-c) for the selected 2-D materials: 3 $p$-like orbitals have been adopted for S, 5 $d$-like orbitals for Mo and W, and 4 $sp^{3}$ hybrid for P.

To handle phonon transport through transistors with a 2-D channel, as in Fig. 1, the following NEGF-based system of equations must be processed

[TABLE]

where the $\mathbf{D}(\omega,q_{z})$’s are the phonon Green’s functions at frequency $\omega$ and momentum $q_{z}$, $\mathbf{\omega}$ a diagonal matrix containing the phonon frequency as single entry, the $\mathbf{\Pi}^{B}(\omega,q_{z})$’s the boundary self-energies, and the $\mathbf{\Pi}^{S}(\omega,q_{z})$’s the scattering ones, while $\mathbf{\Phi}(q_{z})$ refers to the dynamical (Hessian) matrix of the studied domain. The same Green’s Function types as for electrons also exist for phonons (retarded, advanced, lesser, and greater). All blocks involved in Eq. (6) are of size 3$\times$3, which corresponds to the degrees of freedom of the crystal vibrations. The matrix $\mathbf{\Phi}(q_{z})$ is constructed in an analogous way as $\mathbf{H}(kz)$: through density functional perturbation theory (DFPT) dfpt , as available in VASP, the dynamical matrix of a representative unit cell is computed from first-principles. It is then scaled up to obtain the entries for a larger device structure, e.g. the transistor in Fig. 1. The phonon bandstructures of single-layer MoS2, WS2, and black phosphorus are displayed in Fig. 2(d-f).

With the knowledge of $\mathbf{H}(k_{z})$ and $\mathbf{\Phi}(q_{z})$, electron and phonon quantum transport simulations can be performed, first in the ballistic limit of transport, i.e. $\mathbf{\Sigma}^{S}(E,k_{z})$=$\mathbf{\Pi}^{S}(\omega,q_{z})$=0. Eqs. (3) and (6) can be solved with a recursive Green’s Function (RGF) algorithm rgf that produces only the desired entries of the Green’s Functions. The open boundary conditions are derived from the solution of eigenvalue problems with a shift-and-invert technique prb06 . As an illustration, the ballistic electron and phonon transmission functions through 40 nm long monolayers of MoS2, WS2, and black phosphorus (in the armchair configuration) have been computed and the results are plotted in Fig. 3 (a) and (b) for $k_{z}$=0 and $q_{z}$=0, respectively. A flat band potential is assumed in all cases so that a typical step-like behavior of the transmission as a function of the electron/phonon energy can be observed. Attention should be paid to the fact that the mini-gaps (regions with no bands) present in the phonon bandstructures in Fig. 2(d-f) are clearly visible in the transmission functions in Fig. 3(b).

To accurately compute the current that flows through a 2-D transistor, it is not sufficient to keep one single momentum point, e.g. $k_{z}$=$q_{z}$=0. Here, it has been found that a total of 11 $k_{z}$ and $q_{z}$ momentum points is necessary to reliably model the periodicity of the out-of-plane direction $z$ (no more than 1% current variations as compared to simulations with 21 momentum points). It should finally be underlined that electrons and phonons can only enter and leave the simulation domain at the source and drain contact extremities, no escape through the top or bottom oxide of the transistor in Fig. 1 is possible. This restriction certainly causes an overestimation of the lattice temperature apl17 when the electron and phonon populations fully interact with each other, as described in the next Section.

II.2 Electron and Phonon Coupling

To complete the picture it remains to define the scattering self-energies that couple electron and phonon transport and drive both populations out-of-equilibrium. To reduce the computational intensity, only the diagonal blocks of the electron-phonon components $\mathbf{\Sigma}^{\gtrless S}(E,k_{z})$ are retained. It can be analytically demonstrated, as in Ref. jtlu1 , that these scattering self-energies ensure total energy conservation. Each individual block is defined as

[TABLE]

In Eq. (7) $\nabla_{i}\mathbf{H}_{nm}$=$\delta\mathbf{H}_{nm}/\delta(R_{m}^{i}-R_{n}^{i})$ is the derivative of the Hamiltonian matrix block $\mathbf{H}_{nm}$ with respect to variations along the $i$=$x$, $y$, or $z$ axis of the bond $\mathbf{R}_{m}-\mathbf{R}_{n}$ connecting atoms $n$ and $m$. Since the Hamiltonian matrix in the MLWF basis includes connections with the 30 (36) nearest-neighbors (NNs) of each atom in MoS2 and WS2 (black phosphorus), the sum over $l$ in Eq. (7) has to cover the same range of interactions. Practically, summing over the 12 (MoS2 and WS2) or 13 (black phosphorus) NNs is sufficient. Still, although the scattering self-energy $\mathbf{\Sigma}^{\gtrless S}(E,k_{z})$ is block diagonal only, its entries couple one atom at $\mathbf{R_{i}}$ with its 12-13 NNs at $\mathbf{R_{l}}$, thus indirectly accounting for non-diagonal effects. The fact that each momentum $k_{z}$ is connected to all possible $k_{z}-q_{z}$ points also contributes to the presence of non-diagonal effects.

The phonon-electron scattering self-energy matrix $\mathbf{\Pi}^{S}(\omega,q_{z})$ cannot be assumed diagonal as $\mathbf{\Sigma}^{\gtrless S}(E,k_{z})$ because this would violate the energy conservation rule between electrons and phonons. It has to take the following form

[TABLE]

for the diagonal block entries and

[TABLE]

for the non-diagonal ones. In Eqs. (8) and (9), “$\mathrm{tr}\left\{\right\}$” refers to the trace operator and atoms $l$ and $i$ must be coupled with each other via a non-zero matrix element $\mathbf{H}_{li}$. The last missing components are the retarded electron and phonon scattering self-energies, which can be approximated as

[TABLE]

It is obvious from Eqs. (7-9) that the scattering self-energies bridge the electron and phonon populations because $\mathbf{\Sigma}^{S}(E,k_{z})$ depends on $\mathbf{D}(\omega,q_{z})$ and $\mathbf{\Pi}^{S}(\omega,q_{z})$ on $\mathbf{G}(E,k_{z})$. This is also the reason why Eqs. (3-11) must be self-consistently solved till convergence is reached in the so-called Born approximation. The resulting numerical problem is particularly challenging since all energies $E$, frequencies $\omega$, and momentum $k_{z}$/$q_{z}$ are connected altogether. Such systems of equations cannot be tackled on a single computational core and must therefore be treated in parallel. To manage the distribution of the numerical tasks and the gathering of the needed data to calculate the scattering self-energies, the multi-level parallelization scheme introduced in Ref. sc10 has been adapted and enhanced to enable electro-thermal transport simulations of 2-D materials from first-principles.

II.3 Calculation of Observables

A second self-consistent loop must be established between the solution of the Schrödinger and Poisson equations to properly take the electrostatics of the examined devices into account. This requires evaluating the electron concentration $n(\mathbf{R}_{i})$ (here only $n$-type transistors are investigated) for each atomic position $\mathbf{R}_{i}$ according to

[TABLE]

and plugging the outcome into Poisson’s equation expressed on a finite element method (FEM) grid. After convergence the electrical current flowing between two adjacent unit cells labeled $s$ and $s+1$ of the 2-D device structure

[TABLE]

the electron component of the energy current between $s$ and $s+1$

[TABLE]

as well as its phonon part

[TABLE]

can be computed from the electron and phonon Green’s functions. In all these equations, $\hbar$ is Planck’s reduced constant, $q$ the elementary charge, and it is implied that atom $i$ ($j$) belongs to the unit cell $s$ ($s+1$). It should be pointed out in Eqs. (13) to (15) that (energy) current conservation does not only require an integration over energy (or frequency), but also a summation over the momentum ($k_{z}$ or $q_{z}$).

Another quantity of interest is the effective lattice temperature of unit cell $s$, $T_{eff}(s)$. It can be obtained by assuming that the total phonon population within one unit cell $s$, $N_{ph}^{tot}(s)$, which is first calculated with NEGF, locally obeys a Bose-Einstein distribution function $N_{Bose}(T_{eff},\omega)=1/(\mathrm{exp}(\hbar\omega/k_{B}T_{eff})-1)$ with $T_{eff}$ as the governing temperature and $k_{B}$ as Boltzmann’s constant

[TABLE]

where the local density-of-states $\mathrm{LDOS}(\mathbf{R}_{i},\omega,q_{z})$ is defined as

[TABLE]

The effective lattice temperature of unit cell $s$ is then retrieved by matching Eqs. (16) and (17). The resulting non-linear system of equations can be solved, for example, with a Newton-Raphson iterative method. If necessary, an atomic resolution of $T_{eff}$ is also possible.

III Device Results

To illustrate the influence of electro-thermal effects on the performance of $n$-type transistors made of a 2-D material, the single-gate structure in Fig. 1 has been simulated with the models described in Section II together with single-layer MoS2, WS2, and armchair-oriented black phosphorus (AC BP) channels. The logic switch specifications include a gate contact of length $L_{g}$=15 nm, 12.5 nm long source and drain extensions with a donor doping concentration $N_{D}$=6e13 cm*-2*, perfectly ohmic contacts, a top (bottom) oxide layer of thickness $t_{ox}$=3 (20) nm and relative permittivity $\epsilon_{R}$=20 (3.9), and a supply voltage $V_{DD}$=0.67 V.

The $I_{d}$-$V_{gs}$ transfer characteristics of the MoS2 device at a drain-to-source voltage $V_{ds}$=$V_{DD}$ and OFF-current $I_{OFF}$=0.1 $\mu$A/$\mu$m are reported in Fig. 4 for three different situations: (i) ballistic transport ($\mathbf{\Sigma}^{S}(E,k_{z})$=$\mathbf{\Pi}^{S}(\omega,q_{z})$=0 in Eqs. (3-8)), (ii) dissipative transport in the presence of electron-phonon scattering, but with an equilibrium phonon population (only $\mathbf{\Pi}^{S}(\omega,q_{z})$=0), and (iii) in the same configuration as in (ii), but with both the electron and phonon populations driven out-of-equilibrium. The latter case incorporates electro-thermal effects, contrary to the others where the lattice temperature does not vary as a function of the electrical current and remains equal to the room temperature $T_{0}$=300 K. It might look surprising in Fig. 4 that the ballistic current is not the largest one. This peculiar feature has already been discussed in details in Ref. szabo1 : in typical transistor structures, the ballistic transmission function $T(E)$ from left (source) to right (drain) only slightly depends on the applied $V_{ds}$ at a given $V_{gs}$. With strong drain-induced barrier lowering (DIBL), a shift towards lower energy of the point where $T(E)$ turns on is expected, but no shape changes. However, with MoS2 and most other TMDs, the transmission function exhibits a strong dependence on $V_{ds}$, as can be seen by comparing $T(E)$ in Fig. 3 (flat-band conditions, zero electrostatic potential everywhere) and in Fig. 5 (linear potential drop between both device ends). Obviously, at energies close to the conduction band edge, one of the two available channels in the source (sub-bands with a positive velocity) of MoS2 stops conducting at high $V_{ds}$ so that $T(0)$=1 with a linear potential drop, while $T(0)$=2 with flat bands. As a consequence, the ballistic current is reduced. This “anomaly” is due to the specificity of the TMD bandstructures, where sub-bands with a very narrow energy width $\Delta E$ smaller than $qV_{ds}$ can be found, as indicated in Fig. 2(a). These states are not conductive in the ballistic limit of transport, but become active as soon as electron-phonon scattering is switched on szabo1 . In this case an electron occupying a non-conductive band can be transferred to a conductive one by absorbing or emitting a phonon. This explains the rather counterintuitive reinforcement of the dissipative current with respect to the ballistic one in Fig. 4.

When the influence of the phonon-electron scattering self-energy $\mathbf{\Pi}^{S}(\omega,q_{z})$ is added to the picture, phonon emission and absorption processes induce local variations of the lattice temperature. This causes a current decrease by about 10%, as compared to the case with electron-phonon interactions only, which is marked by a double arrow in Fig. 4 and labeled self-heating. To better understand what happens inside the MoS2 transistor, the spectral electron and phonon current distributions in Fig. 6(a-b) are investigated as next step. It can be observed in sub-plot (a) that electrons loose a substantial portion of their energy through phonon emission after passing the maximum of the electrostatic potential energy barrier, when they are accelerated by the strong electric field on the drain side. From sub-plot (b) it can be inferred that at the location with the highest phonon generation rate (and therefore the highest lattice temperature too, as explained later) at around $x$=31 nm, the phonon current is equal to 0 and the created lattice vibrations either propagate towards the drain ($x$=40 nm, positive current) or the source ($x$=0 nm, negative current) contact. The phonons that reach the source can interact there with electrons moving towards the drain and increase their backscattering rate prb14 , thus leading to the current reduction visible in Fig. 4. To complement this analysis, it is demonstrated in Fig. 6(c) that energy is conserved in the developed simulation approach. The sum of the electron and phonon energy currents, as given in Eqs. (14) and (15), is indeed constant from source to drain. This indicates that the energy lost by electrons is correctly transferred to the phonon bath and vice versa. Slight fluctuations of the total energy current still exist (less than 1%): they are due to the slow convergence of Eqs. (3-10). It has been verified that more self-consistent iterations between the Green’s Functions and the scattering self-energies improve the situation, at the expense of the computational time.

With Eqs. (16) and (17), the non-equilibrium phonon population can be converted into a position-resolved effective lattice temperature $T_{eff}$ that is resolved at the unit cell level. The results are provided in Fig. 7 for the MoS2 transistor. What is shown is the lattice temperature increase $\Delta T$=$T_{eff}$-$T_{0}$ as a function of the location along the $x$-axis and as a function of the applied gate-to-source voltage $V_{gs}$. Two points should be emphasized: first, $T_{eff}$ peaks at the same $x$-coordinate as where the phonon current vanishes in Fig. 6(b) and where the phonon generation rate reaches its maximum. The formed hot spot is in fact situated close to the channel-drain interface, in the region with the sharpest drop of the electrostatic potential energy and therefore the highest phonon emission probability. Secondly, the $\Delta T$ values are extremely high $>$400 K at $V_{gs}$=0.7 V. Before approaching such temperatures, the MoS2 single-layer crystal would already have oxidized (oxidation temperature: 675 K) zhou and the whole device would have broken down. As a comparison, in Ref. pop2 , the $\Delta T$ of a micrometer-scale monolayer MoS2 transistor has been estimated through Raman thermometry measurements to be around 250 K for a much larger $V_{ds}$=30 V vs. 0.67 V here, but a smaller current $I_{d}\approx$200 $\mu$A/$\mu$m vs. $\sim$1000 $\mu$A/$\mu$m here and a higher contact resistance $R_{c}$=5k$\Omega\times\mu$m vs. perfectly ohmic contacts here. This confirms that self-heating effects are very important in monolayer MoS2, but that the simulated temperature increases are overestimated, mainly because no phonon escape through the surrounding oxide layers is allowed. Such mechanisms would definitively lead to a decrease of the effective lattice temperature, but their exact contribution is difficult to quantify.

Keeping in mind that $\Delta T$ is overestimated, it is still relevant to compare the electro-thermal properties of MoS2 with those of other 2-D materials, namely WS2 and AC black phosphorus, and with those of more conventional ultra-scaled Si nanowire transistors since the same set of approximations is applied everywhere, i.e. phonons escape at the source and drain only, not through oxide layers. The goal of this study is to determine what material/design is the least sensitive to self-heating and why. In Fig. 8(a) the electrical power dissipated inside the channel of all switches, $P_{diss}$, is plotted as a function of the drive current $I_{d}$. For the Si nanowire transistor, a $<$100$>$-oriented structure with a diameter $d$=3nm, a gate length $L_{g}$=15 nm, and a gate-all-around configuration has been selected and the data from Ref. prb14 have been recalled. Here, $P_{diss}$ is defined as the difference in the electron energy current (see Eq. (14)) between the source and the drain. Normally, it is expected that the total dissipated power is equal to $V_{ds}\times I_{d}$ and that it has the same value for all structures at a given $I_{d}$. However, due to the short length of the considered devices, only a fraction of the total power dissipation takes place in the channel region, the rest in the contacts. This explains why $P_{diss}<V_{ds}\times I_{d}$ in our simulations and why the curves are different in Fig. 8(a).

The amount of dissipated power in a given component can be put in relation with the phonon-limited mobility $\mu_{ph}$ of its underlying channel material nl16 : the higher $\mu_{ph}$ the less phonons are emitted and the less electrical power is converted into heat. To verify whether this statement is valid or not here, the phonon-limited electron mobility of single-layer MoS2, WS2, and AC black phosphorus, as computed with the “dR/dL” method rim , is presented in Fig. 9 as a function of the carrier concentration. Due to the heavier atomic mass of tungsten as compared to molybdenum, the amplitude of the crystal vibrations and in turn the electron-phonon coupling strength is weaker in WS2 than in MoS2, which gives a higher mobility and a lower dissipated power at a fixed current magnitude, as postulated above. This demonstration fails however when AC black phosphorus comes into play: its mobility is smaller than that of MoS2 and WS2 due to more significant crystal vibrations, but it still dissipates less power. The reason behind this apparent inconsistency can be traced back to the distinctive bandstructure of TMDs, where sub-bands with a narrow energy width cannot carry current at high $V_{ds}$, as shown in Fig. 5. As a result, electron-phonon interactions are needed to make those sub-bands conductive and to reach high current densities. The number of phonon emission processes occurring in TMDs under the application of a high $V_{ds}$ is therefore much more important than predicted by the mobility, which is computed under flat-band conditions. While electron-phonon interactions contribute to a current increase in TMDs, they simultaneously consume electrical power and lead to the higher heat dissipation of MoS2 and WS2 seen in Fig. 8(a). Although the Si nanowire transistor has about the same phonon-limited electron mobility as AC black phosphorus (300 cm2/Vs at an electron concentration $n_{D}\approx$1e13 cm*-2*), it dissipates more power. Since the crystal and device structures are quite different, it is not possible to explain this behavior with the same line of arguments as before.

Figure 8(b) then reports the maximum effective lattice temperature increase extracted from the same transistors as before versus the electrical current $I_{d}$. What matters in this sub-plot is no more the number of phonons that is generated (heat dissipation), but once they have been created how efficiently they can leave the active region of the logic switch. If the emitted phonons have a high group velocity, they more rapidly diffuse away from the region where they originate such that the effective lattice temperature at this location cannot rise as much as it would have if the phonon population would keep accumulating. Hence, the thermal transport properties determine to a large extend the magnitude of the self-heating effects. For example, in WS2, less power is dissipated at a given $I_{d}$ than in MoS2, but the lattice temperature increase $\Delta T$ is the same in both 2-D materials: less crystal vibrations are generated in WS2 due to the weaker electron-phonon coupling and lower electron effective mass, but they propagate slower due to the higher atomic mass of tungsten. The poorer thermal transport properties of single-layer WS2 are indirectly illustrated in Table 1, which gives the sound velocity of all studied 2-D crystals, as derived from ab-initio calculations. It can also be seen in Fig. 8(b) that the temperature increase is lower in AC black phosphorus and in the Si nanowire than in both TMDs.

Finally, the maximum lattice temperature increase can be represented as a function of the power dissipated inside the transistor channel. The results are shown in Fig. 8(c). While the curves for the 2-D materials are not identical, they all exhibit a much higher lattice temperature increase than the Si nanowire at a given $P_{diss}$. The strong confinement of phonons in 2-D crystals leads to rather poor thermal transport properties and prevents the rapid evacuation of heat from the active region of 2-D logic switches. As mentioned earlier, allowing phonons to escape through the surrounding oxide layers would certainly reduce the $\Delta T$, but the same would happen in the Si nanowire where phonons are also confined in the semiconductor apl17 . A recent experimental study on single-layer MoS2 pop2 has clearly established that self-heating is much more significant in this 2-D material than in Si, which qualitatively agrees with the finding of this paper.

IV Conclusion

An ab-initio electro-thermal transport modeling approach has been introduced in this paper and applied to the simulation of single-gate transistors with a single-layer 2-D material as channel. After demonstrating the functionality of the implemented electron and phonon quantum transport solver, the effective lattice temperature and the power dissipation of logic switches made of MoS2, WS2, AC black phosphorus monolayers as well as of a Si nanowire have been computed and compared to each other. It turns out that the considered 2-D semiconductors are much more affected by self-heating than Silicon, especially at high current densities. This might become a serious problem in the future when the contact series resistances of devices with a 2-D channel material will shrink and their drive current increase. Enhancing heat dissipation through the bottom oxide layer could be the solution to mitigate this issue, for example by carefully chosen the oxide material or by making it thinner. Further model developments will be needed to account for these effects and support the on-going experimental activity.

Acknowledgment

This work was supported by the Swiss National Science Foundation Grant No. PP00P2_159314, by ETH Research Grant ETH-32 15-1, by the European Research Council under Grant Agreement No 335684-E-MOBILE, and by CSCS under Projects s662 and s714.