Mobility of 2D materials from first principles in an accurate and automated framework

# Mobility of 2D materials from first principles in an accurate and automated framework

## Abstract

We present a first-principles approach to compute the transport properties of 2D materials in an accurate and automated framework. We use density-functional perturbation theory in the appropriate bidimensional setup with open-boundary conditions in the third direction. The materials are charged by field-effect via the presence of planar counter-charges. In this approach, we obtain electron-phonon matrix elements in which dimensionality and doping effects are inherently accounted for, without the need for post-processing corrections. The framework shows some unexpected consequences, such as an increase of electron-phonon coupling with doping in transition-metal dichalcogenides. We use symmetries and define pockets of relevant electronic states to limit the number of phonons to compute; the integrodifferential Boltzmann transport equation is then linearized and solved beyond the relaxation-time approximation. We apply the entire protocol to a set of much studied materials with diverse electronic and vibrational band structures: electron-doped MoS, WS, WSe, phosphorene and arsenene, and hole-doped phosphorene. Among these, hole-doped phosphorene is found to have the highest mobility, with a room temperature value around cmVs. We identify the factors that affect most the phonon-limited mobilities, providing a broader understanding of the driving forces behind high-mobility in two-dimensional materials.

## I Introduction

The scientific and engineering community is devoting a major effort towards the identification and fabrication of novel 2D materialsEfetov2010 (); Radisavljevic2011 (); Lebegue2013 (); Hanlon2015 (); Yasaei2015 (); Gibaja2016 (); Ares2016 (); Bandurin2016 (); Ashton2017 (); Cheon2017 (); Mounet2018 () and their use in electronic devicesWang2012 (); Chhowalla2016 (). To accelerate the discovery of the best candidates for electronic transport, systematic and accurate methods to compute phonon-limited resistivity and mobility from first principles would be very beneficial. Nevertheless, these quantities are not straightforward to compute and various degrees of approximations are often needed.

A very successful approach relies on the solution of the Boltzmann transport equationZiman (); Allen1978 () (BTE), in which the electron-phonon coupling (EPC) enters the expression for the scattering ratesGrimvall (); Giustino2017 () and can be computed from first principles within density-functional perturbation theory (DFPT)Savrasov1994 (); Liu1996 (); Mauri1996 (); Bauer1998 (); Baroni2001 () or with finite-differences schemesDacorogna1985 (); Gunst2016a (). Although both methods proved to be very successful, DFPT is more efficient in dealing with long-wavelength effects since it does not need large supercells to accommodate the perturbation, and will thus be the method of choice in the following.

Either way, the computation of EPC in 2D materials from first principles presents some challenges. Being ruled by long-range Coulomb interactions, the long-wavelength electron-phonon dynamics are highly dependent on dimensionality. The most obvious example is the Fröhlich interaction in polar materials, which diverges in the long-wavelength limit in 3D while remains finite in 2D Sohier2016 (). This fundamentally different behaviour is difficult to capture in standard electronic-structure codes based on plane-wave basis sets because of the spurious interaction with artificial periodic images, which leads to erroneous results when . Thus, special care must be adopted when computing the response of a 2D system to long-wavelength perturbationsKozinsky2006 (); Sohier2015 (); Sohier2016 (); Sohier2017nl ().

An appropriate treatment of doping is an additional challenge. While field-effect charging is omnipresent in experimental setups, it is usually not included in first-principles simulations. Often, calculations are performed for the neutral case, and the Fermi level is shifted within a frozen-band approximation when integrating the EPC matrix elements to obtain a certain scattering-related quantity. This is approximate because the EPC itself depends on doping, most notably via screening, but also via other mechanisms that will be discussed later. A more realistic description is now possible thanks to a recent development of DFPT for gated 2D materialsSohier2017 (), where the doping charge is neutralized by adding planar distributions of counter-charges, mimicking experimental conditions where 2D materials are doped electrostatically through gates.

Another source of challenge is that, in general, a very large number of EPC matrix elements must be computed to obtain a fine enough sampling of the electron-phonon dynamics over the Brillouin zone, which is computationally very expensive. To boost efficiency, one may use Wannier interpolations Lee2005 (); Yates2007 (); Giustino2007 (); Calandra2010 (); Marzari2012 (); Ponce2016 () to obtain inexpensively the EPC over a dense mesh. However, interpolations work seamlessly only with short-range interactions. Long-range effects, such as Fröhlich or piezo-electric coupling, need to be modeled and treated separately; this has been solved for three-dimensional systemsSjakste2015 (); Verdi2015 (), while in 2D, Fröhlich coupling has been modeledSohier2016 () and efficient ways to interpolate the phonon dispersions have been put forwardSohier2017nl (), but a thorough description of EPC interpolation is still missing. Nonetheless, the absence of interpolation techniques might not be, in the long run, a major constraint. Indeed, the reduced dimensionality entails one less dimension to sample, with a drastic reduction of electronic states and phonons to consider. Thus, it can be argued that for 2D materials, one should focus on addressing dimensionality and charging first, rather than interpolation.

Here, we thus choose a direct approach and compute all EPC matrix elements from DFPT in the appropriate boundary and charging conditionsSohier2017 (). The use of symmetries and energy selection rules, in addition to the momenta being restricted to two dimensions, makes the process computationally feasible.

For the solution to the Boltzmann transport equation for electrons and holes, different approaches involving various degrees of approximation have been put forward. In the case of metals, Allen derived an approximate solutionAllen1978 () by suitably modifying Eliashberg theoryGrimvall () to transport. This approach has been successfully applied in the early efforts to compute mobilities from first principlesSavrasov1996 (); Bauer1998 () and it is now customary in many available transport codes such as EPWPonce2016 (). In the general case, including semiconductors, the major challenge is that the integrodifferential BTE does not have a closed-form solution. Most approaches use then some form of the relaxation-time approximationRestrepo2009 (); Shishir2009 (); Borysenko2010 (); Jin2014 (); Restrepo2014 (); Park2014 (); Gunst2016a (); Jin2016 (); Rudenko2016 (); Trushkov2017 (); Gaddemane2018 () to make it closed.

Still, the errors associated with the various approximate solutions to the BTE in the literature are often difficult to quantify. An iterative scheme to solve the inelastic part of the BTE Kaasbjerg2012a (); Kaasbjerg2013 (); Li2015 () has been proposed recently, while a very efficient preconditioned conjugate-gradient approach has been reported in Ref. Fiorentini2016, following a recipe introduced in the context of the phonon Boltzmann equation for thermal transportFugallo2013 (); Fugallo2014 (); Cepellotti2015 (). Ref Li2015, also offers an interesting comparison of the different methods to solve the BTE and shows a broad agreement for MoS, but this might not be the case for all materials. A full numerical solution to the BTE beyond the relaxation-time approximation was introduced in Ref. Sohier2014a, for graphene. Here, we propose its generalization to any 2D material.

Relatively few 2D materials have been theoretically investigated up to now. Graphene has been studied extensivelyAndo2006 (); Hwang2007 (); Manes2007 (); Samsonidze2007 (); Stauber2007 (); Fratini2008 (); Hwang2008a (); Shishir2009 (); Hwang2009 (); Borysenko2010 (); Mariani2010 (); DasSarma2011 (); DasSarma2011 (); Kaasbjerg2012 (); Hwang (); Restrepo2014 (); Park2014 (); Sohier2014a (); Gunst2016a (); Gunst2017 (), showing excellent agreementPark2014 () with experimentsEfetov2010 () and a detailed understanding of the main processes limiting mobilityBorysenko2010 (), including the effects of dimensionality and charging by field-effectGunst2017 (); Sohier2017 (). MoS, another prototypical 2D material, has also been studied in several works Gunst2016a (); Kaasbjerg2012a (); Kaasbjerg2013 (); Li2015 (); Jin2014 (); Restrepo2014 (), as well as phosphorene Gaddemane2018 (); Qiao2014 (); Jin2016 (); Trushkov2017 (); Rudenko2016 (); Liao2015 (), arsenene Wang2015 (); Pizzi2016nc (); Wang2017 (), silicene Gunst2016a (), and other TMDs Jin2014 (). On the other hand the effects of periodic images and dimensionality are explicitly treated only in some of the first-principles effortsGunst2016a (); Kaasbjerg2012a (); Kaasbjerg2013 (). Charging is most often treated as a rigid shift of the Fermi energy in the computation of the transport properties, but not included in the computation of the matrix elements themselves. Consequently, the influence of doping on the EPC is generally neglected. In some instances, electronic screening is accounted for analyticallyKaasbjerg2012a (); Kaasbjerg2013 (), but this necessarily entails the use of models and approximations (not obvious in 2D). Also, as will be shown here, screening is not the only effect of doping on EPC matrix elements. Last, computational accuracy is often limited by the very expensive nature of the Brillouin zone (BZ) integrals. For all these reasons, a full treatment of doping and periodicity is necessary, together with an efficient and automatic implementation of all BZ sums.

All these points are addressed in this paper, that is structured as follows. In the first section, we describe the formal framework of the BTE and identify the quantities needed to solve it. In the second section, we present the workflow of first-principles calculations used to compute the physical quantities associated with phonon-limited transport. Then we discuss the results of this workflow applied to a set of known prototypical materials, i.e., three electron-doped TMDs in their 2H form (MoS, WS, WSe), as well as electron-doped arsenene, and both hole- and electron-doped phosphorene.

## Ii Boltzmann Transport Equation

In this first section, we review the framework of the Boltzmann transport equation for electrons in an effort to settle the context and introduce the various quantities used in the rest of the paper. Similar derivations can be found in the literature, e.g. in Ref. Ziman (); Ashcroft1976 (); Nag1980 (). We consider any 2D material to be in the plane with an applied electric field in the same plane, in the direction. The electric field favors states with wavevectors in the opposite direction, , thus taking the electronic distribution out of its equilibrium Fermi-Dirac ground-state . Phonon scattering acts to bring the system back towards its unperturbed equilibrium state; then, a steady-state regime is reached and a net electric current emerges. The conductivity, or the inverse of the resistivity, is then defined asZiman (); Ashcroft1976 (); Nag1980 ():

 σ=1ρ=∑α2e|E|∫k∈αdk(2π)2f(k)v(k)⋅uE (1)

where the factor 2 accounts for spin degeneracy, is the Coulomb charge, the electric field, is an index representing the different valleys in the Brillouin zone, is the steady-state occupation function, and the velocity of the electronic state. Here, rather than integrating over all the Brillouin zone, we limit ourselves to the relevant electron or hole pockets that are occupied by temperature or doping. In semi-conductors, those pockets form valleys which we will further define later, and wavevectors are taken from these valleys. Mobilities can then be obtained from the Drude model as , where is the electron density (replace with hole density in case of hole doping). The central quantity to obtain is thus the occupation distribution . Assuming this to be spatially uniform and time-independent (steady-state), the Boltzmann transport equation states that the change of the occupation distribution driven by the electric field must be compensated by scatteringZiman (); Ashcroft1976 (); Nag1980 ():

 −eEℏ⋅∂f∂k=(∂f∂t)scatt(k) (2)

The collision integral (right-hand side) can be found using Fermi’s golden rule:

 (3)

where is the scattering probability from state to state ; should include in general all relevant scattering processes. Here, we sum over all electron-phonon scattering probabilities associated to the phonons modes of the system., i.e. , where is the sum of phonon emission and absorption termsGrimvall () for phonon mode of momentum :

 Pkk+q,ν=2πℏ1N|gkk+q,ν|2{nq,νδ(εk+q−εk−ℏωq,ν)+(nq,ν+1)δ(εk+q−εk+ℏωq,ν)}. (4)

In the expression above, is the electron-phonon coupling (EPC) matrix element and is the phonon occupation, which we assume to be the equilibrium Bose-Einstein distribution. The functions stem from the energy selection rules involved in the scattering. They imply the necessity of evaluating the EPC on very fine momentum grids, making the calculations challenging. Phonon scattering as described in Eq. 4 involves three quasi-particles: an initial electronic state, a final electronic sate, and a phonon. The EPC matrix element connecting the initial and final states can be computed within DFPTBaroni2001 (); Giustino2017 () and reads:

 gk,k+q,ν=∑a,iea,iq,ν√ℏ2Maωq,ν⟨k+q|∂VKS(r)∂ua,i(q)|k⟩ (5)

where is an atomic index, a cartesian index, is the phonon eigenvector, is the mass of atom , are the initial and final electronic states, and is the derivative of the Kohn-Sham potential with respect to a periodic displacement of atom in direction .

We adopt a perturbation approach at first order in electric field and write :

 f(k)=f0(k)+f1(k) (6)

where is the Fermi-Dirac function and is the linear perturbation proportional to the electric field. By replacing the occupation distribution of Eq. 6 in Eq. 11, and keeping only first-order terms in electric fields while using the detailed balance condition , we getAshcroft1976 (); Nag1980 ():

 −e|E|ℏ∂f0∂εℏv(k)⋅uE=∑k′Pkk′1−f0(k′)1−f0(k)×{f1(k′)f0(k)(1−f0(k))f0(k′)(1−f0(k′))−f1(k)} (7)

We now need to solve the above equation to find . It is convenient to make the following general ansatz

 f1(k)=e|E|uE⋅F(k)∂f0(k)∂ε (8)

where is a vectorial quantity with units of length that can be understood as a mean free displacementLi2015 (). In the following the only approximation we make is to assume that this mean free displacement is along the band velocity, that is we write:

 F(k)=v(k)τ(k) (9)

where is unknown and has the dimension of time. We will call it scattering time, although it is not the scattering time as often understood in the context of the relaxation-time approximation. Now, putting the ansatz Eq. 8 and 9 in Eq. 7, we get:

 (10)

Knowing that , we can further simplify to get what we will refer to as the linearized Boltzmann transport equation (BTE):

 (11)

where we keep the term on the left-hand side to avoid having it on the right-hand side as a denominator, which would bring some numerical instability in the process of solving the equation numerically. Putting back the distribution Eq. 6 in the expression for the conductivity Eq. 1 we get:

 σ=1ρ=∑p2e2∫k∈αdk(2π)2(v(k)⋅uE)2τ(k)∂f0∂ε (12)

A very minimal set of approximations have been made to get to Eq. 11 and 12 (linear order in electric field, steady state, equilibrium phonon distribution, mean free displacement along band velocity). We will compute the conductivity at this level of approximation. This goes beyond the relaxation time approximation, which would correspond to neglecting the second term in brackets on the right-hand side of Eq. 11, and differs from the approaches in Ref. Li2015, and Fiorentini2016, only in the assumption of Eq. 9. Our approach to solve the BTE (Eq. 11) will be to discretize over electronic states and solve a linear system, equivalent to inverting a matrix. The general formulation of this “matrix inversion” solution is described in the following.

## Iii Workflow / Calculations

In this section we describe the computational workflow developed to calculate the transport properties of 2D materials. Electron-doped WS ( cm) is used as a case study; this is a relatively complex system due to its multi-valley nature. As illustrated in Fig. 1a), the process can be separated into the two following steps:

• EPC: We compute the linear response of the system with respect to a set of phonon momenta (defined in the following section, along with all other momenta sets). The resulting EPC matrix elements are then projected on a set of relevant initial states and interpolated on a set of final states . The initial states are where we want to evaluate the scattering time. Final states are all the states accessible from the initial states via phonon scattering.

• Transport: The matrix elements are then used to compute the scattering probabilities and solve the BTE, which yields the scattering time for each initial state in . The scattering time is then interpolated for all states in , as shown in Fig. 1 c). The integration of the scattering times gives the transport quantities, like the temperature-dependent mobilities and resistivities shown in Fig. 1 d).

In the following we first detail and justify the sampling choices for the momentum sets , and . Then we describe each of the two steps outlined above.

### iii.1 Sampling

Phonon scattering involves three quasiparticles: an initial electronic state, a final electronic state and a phonon. We use a different momentum sampling over the Brillouin zone (BZ) for each of those quasiparticles, according to their use in the workflow and the cost of the associated calculations.

Only a subset of the electronic states are relevant for transport: we see from the expression of conductivity (Eq. 12) that one needs to find the dependency of the perturbed distribution only for electronic states for which is significantly different from zero: this represents a set of electronic states with energies in a range of a few around the Fermi level. Furthermore, we see in the collision integral that the scattering time at depends on the scattering time of possible final states at . This means that the energy range where we need to evaluate must be extended by the maximum phonon energy above and below. Such constraints define pockets or valleys in the BZ, and in WS we keep all electronic states up to an energy eV above the bottom of the conduction band. Panel c) of Fig. 1 shows these valleys in the Brillouin zone for WS. We distinguish two kinds of valleys: two K-valleys around the high-symmetry points K and K’, and 6 valleys around the Q points, situated approximately halfway between and K.

We define a fine grid on the pockets, see Fig. 2, on which quantities related to the band structure, i.e. energies and velocities are computed. A fine grid is needed to evaluate the velocities (as gradients of the energies) properly and to converge the integrals in Eq. 11 and Eq. 12. To compute the eigenenergies on a fine grid, the ground-state electronic density is first computed on a relatively coarse grid (); then, non self-consistent calculations are performed on a finer grid. Since the second step is relatively inexpensive, one can afford very fine sampling. We always define the grid such that the distance between two consecutive k-points is Å (or as close as possible to this value using a Monkhorst-Pack mesh). For example, in WS, this corresponds to 226226 k-points. The scattering times are interpolated on this grid and integrated to get the conductivity, resistivity, or mobility, as shown in Fig. 1c-d). In this work, more focused on accurate linear response, band structures are computed within DFT. However, the same approach could use bands computed at a higher level of theory, such as hybrid functionals, many-body perturbation theory (e.g. GW associated to a Wannier interpolation technique Marzari2012 ()), or applying a correction to the fundamental band-gap obtained from an evaluation of the derivative discontinuity as in GLLBsc functionalsgllb1 (); gllb2 (), or even bands fitted on experiments Cohen1966 (); Wang1995 ().

As written in Eq. 11, the BTE gives the value of at an initial state as a function of at all possible final states linked to initial states by a phonon-scattering event. As will be explained in Sec. III.3, the sampling of initial and final states needs not be the same, and we simply need a map from the final states to their closest symmetry-equivalent initial state. The set of final states ultimately defines the set of scattering events we account for, represented by the scattering probability . These quantities have low symmetry since they combine the initial and final states with the phonon that links them. Thus, the final states of the BTE need to cover the full extent of the valleys. Furthermore, one needs to account for the sharp variations of due to the energy selection rules in Eq. 4. We thus use the fine sampling of the pockets for the final states in the BTE. This simply means that the EPC eventually needs to be interpolated on this grid. The set of momenta is where we want to compute to scattering times by solving the BTE, and it is convenient to note that the scattering times have the symmetry of the electronic eigenenergies . This is not obvious because the computation of the resistivity involves the electric field and phonons which break the symmetries of the band structure. The integrand in Eq. 12, for example, has the symmetry of the band structure plus the electric field. This quantity and the aforementioned have thus lower symmetry than the electronic states. On the other hand, the form of the ansatz for the perturbed distribution (Eqs. 8 and 9) is such that the effect of the electric field is separated out and essentially represents the dependency of the perturbed distribution on the electronic states only. Further analysis of the BTE (Eq. 11) shows that depends on a sum over all possible final states. While each term of the sum might not have the symmetries of the band structure (, e.g., does not), the sum is invariant under the symmetry transformations associated with the BZ. So, since has the symmetry of the band structure, the initial states are sampled from an irreducible representation of the valleys. This sampling should be chosen fine enough to capture the variations of the scattering time; however, the size of defines the size of the linear system to solve for the BTE and it is a factor in the total number of EPC matrix elements to compute. Those calculations have a non-negligible computational cost, and using the irreducible states in the fine grid of the pockets would be unnecessarily expensive. As shown in Fig. 2, we use a coarser grid: for all materials in this work we found that a grid with a k-points spacing approximatively times larger than the fine grid leads to a converged solution to the Boltzmann transport equation at a reasonable computational cost. For example, in WS2, this leads to 126 irreducible initial states, while using the fine pockets grid would have led to more than 800 initial states.

The electronic wave functions in Eq. 5 can be easily recomputed from the ground-state charge density, while the rest is associated with a given phonon perturbation, and is the bottleneck in terms of computational cost (in particular the calculation of ). There are methodsGiustino2007 (); Calandra2010 (); Ponce2016 (); Sjakste2015 () to interpolate the perturbation to the Kohn-Sham potential, but the spirit of the current approach is to compute this quantity directly in DFPT. We will, however, take full advantage of symmetries and use a reasonable sampling of phonon momenta. Each scattering event from an initial to a final states defines a phonon momentum , where spans the pockets , and spans the irreducible states in the pockets. We use the same support as for and , but define a coarser grid, as shown in Fig. 2. We then obtain the set of phonons momenta , where we actually compute the linear response to the phonon perturbation, by the following process: first, we find all possible vectors connecting all final states and initial states in the coarse grid. For these, there are many duplicates, since a single momentum can link several pairs of initial and final states, and so we remove those duplicates. The points are then reduced by symmetry. Indeed, from one -point, the linear response code (as implemented in Quantum ESPRESSO) allows, as a post-process, the computation of for all points in the set defined as all where is a symmetry operation of the crystal. The symmetry reduction is done with a tolerance of about the grid step, since the momenta in do not always fall exactly on the grid. The number of phonons to compute thus depends on the initial grid chosen for the electronic states in a non-straightforward way. We go through the above process several times until we reach a number of phonons that is largely sufficient to capture the variations of EPC while staying reasonable in terms of computational cost. For example, in WS, this results in 200 q-points where to compute the linear response . For materials with less valleys (such as phosphorene studied in Section IV), less than 100 phonons are needed.

In Fig 3 we show both the set of irreducible phonon momenta and the relevant phonon momenta in that lead to final states in the BZ for at least one of the initial states.

### iii.2 Phonons and EPC calculations

Phonons are computed using the recent implementation of DFPT for gated 2D heterostructures Sohier2017 () in Quantum ESPRESSOGiannozzi2017 (); Giannozzi2009 () (QE). This development includes two important modifications: i) a cutoff of the Coulomb interactions in the non-periodic directionRozzi2006 (); IsmailBeigi2006 () and ii) the inclusion of gates to simulate charging of the material in a field-effect setup (FET). The cutoff is necessary to properly account for dimensionality effects Sohier2016 (); Sohier2017nl (): as QE relies on 3D periodic-boundary conditions, there would always be artificial periodic images of the 2D system, and the Coulomb cutoff suppresses spurious interactions between them. The FET setup allows the simulation of the charging of the material in an electrostatic environment that aims to be more realistic than the standard approach of introducing a compensating background of charge uniformly distributed over the full simulation cell (vacuum included). So, we compute the linear response to the phonon perturbations and the dynamical matrices for using these two unique features.

These quantities are obtained for all in by applying symmetry transformations. Then, we compute the EPC matrix elements in Eq. 5 for the states that are in the irreducible initial states . This process is automatized via the AiiDA materials’ informatics infrastructurePizzi2016 () to manage the calculations and store the data. Each symmetry transformation of and its application to the initial states represents a relatively fast run of the Phonon code of QE. Thanks to AiiDA, such operations can be performed in a highly parallelized fashion. For each initial state, a list of the vectors and the corresponding values of for each mode is stored; the are then interpolated linearly to obtain for on the finer grid of the pockets . Fig. 1b) shows the result of this procedure for electron-doped WS. The data storage and provenance provided by AiiDA ensure that all the information collected from the computationally demanding phonon calculations is safely kept and easily re-used for multiple applications (e.g. one can study phonon-mediated superconductivity starting from the same electron-phonon matrix elements).

### iii.3 Transport

Suppose sampling yields irreducible initial states, while sampling yields final states. Writing the BTE for the irreducible states yields equations and unknown variables . To solve the corresponding system, we first need to fold it back on the irreducible states, using the fact that where is the symmetry equivalent of :

 (1−f0(~k))v(~k)⋅uE=∑k′P~kk′(1−f0(k′))v(~k)⋅uE×τ(~k)−∑~k′⎛⎝∑k′≡~k′P~kk′(1−f0(k′))×v(k′)⋅uE⎞⎠τ(~k′). (13)

Considering all the leads to equations for unknown variables. The size of the system to solve thus depends on the sampling of the initial states.

To arrive at the final form of the BTE, we address the fact that its solution is ill-defined when approaches zero. Indeed, considering , any satisfies the above equation. In practice, this brings numerical noise. This situation can happen: i) if , which is relatively rare in practice (e.g. only when is very close to the extremum of a valley); ii) if , where is the direction of the electric field. Situation i) can be treated approximately, without much consequences on the transport results (indeed, states with zero velocity do not contribute to the transport), and we compute the scattering time using the closed form of the BTE shown in App. B. In situation ii), we do need a consistent evaluation of . Indeed, even though does not contribute to the conductivity because its velocity is perpendicular to the field, this is not necessarily true for all that are equivalent to according to the symmetries of the BZ. We thus use the following technique: since has the symmetry of the band structure, it does not depend on the direction of the electric field. In other words, any choice for the direction of the electric field gives the same ’s. Thus, we are free to use any direction for the electric field to solve the BTE, and we can choose different directions for each of the equations corresponding to the points. For each , we choose the direction , where is taken from the extremum of its valley, since virtually never happens.

We then have a well-behaved set of equations with unknowns. Written as a matrix-vector product, it reads:

 ∀~k∈I, F~k=∑~k′S~k,~k′τ~k′, (14)

where

 F~k =(1−f0(~k))v(~k)⋅u~k (15) S~k,~k′=(∑k′′P~kk′′(1−f0(k′′))v(~k)⋅u~k)δ~k,~k′ −⎛⎝∑k′≡~k′P~kk′(1−f0(k′))×v(k′)⋅u~k⎞⎠ (16)

Solving the BTE gives us the values of on the grid. It is then rotated according to the symmetries of the band structure, and interpolated linearly on the fine grid of the pockets . Note that once we have the angular and energy dependent , which does not depend on the electric field in itself, we can compute the conduction integral for an arbitrary direction of the electric field, thus probing any diagonal element of the conductivity tensor. The functions in Eq. 4 are replaced by Gaussians, and a convergence test is performed on the associated smearing.

## Iv Results

We have applied our approach to a set of 6 different cases: 5 electron-doped materials (WS, MoS WSe, arsenene, and phosphorene), as well as hole-doped phosphorene. These are common 2D systems often praised for their potential for transport applications. Also, this set leads to an interesting diversity of band structures, in terms of valleys, their symmetries and their energetic accessibility. Let us stress that we are working here in the framework of DFT, rather than higher levels of theory, and without spin-orbit coupling (SOC). Both approximations can of course affect the bands of the materials studied here. For SOC, we expect the variations to be relatively small for electron-doped arsenene and phosphoreneQiao2014 (). The case of TMDs is a more delicate one: the relative positions of the K and Q valleys seems to be sensitive to many aspects of the calculations, like the inclusion of SOCBrumme2015 (), choice of pseudo potentials, lattice parameters, level of theoryShi2013a () (DFT versus GW) or dopingBrumme2015 (). Thus, while we will be able to compare between the three TMDs and discuss qualitative trends, we do not claim to be quantitative with respect to experiments. That being said, the methodological approach presented in this work is very flexible. One can combine the electron-phonon coupling matrix elements found in DFPT with the band structure found by any mean; e.g., one can use a band structure computed at the highest level of theory, and/or correct it to better fit experiments, before inserting it in the workflow. Regarding SOC, it can be easily included in the band structure calculations and EPC calculations; further development is planned to include it in the solution to the BTE.

As previously mentioned, first-principles calculations are performed with the Quantum ESPRESSO distribution Giannozzi2017 (); Giannozzi2009 () including a 2D Coulomb cutoff and the possibility to charge the material with gates Sohier2017 (), and using the SSSP Accuracy libraryPrandini2018a (); Hamann2013 (); Vanderbilt1990 (); DalCorso2014 (). Structures are taken from the database described in Ref. Mounet2018, . To build this database, structural relaxations were performed in the neutral material, using the SSSP library as well and a k-point sampling corresponding to a spacing of 0.2 Å in each direction. We use a symmetric double-gate setup to charge the materials with a density of cm electrons or holes for all systems except for electron-doped phosphorene, where we choose a lower density of cm to avoid raising the Fermi level too high in the conduction band, where additional valleys appear. Each gate carries half the opposite charge of the material, such that the electric field has equal norm but opposite direction on each side of the material. We could also have used a single-gate setup, with an electric field only on one side, but no large difference is expected since the impact of the electric field setup should come mainly from its effect on spin-orbit coupling, which is not included here. Barrier potentials are added to avoid leakage of electrons towards the gates. These also lead to a hardening of the ZA phonons in the long-wavelength limit, with a non-zero frequency at around cm for the materials considered. Ground state and linear-response calculations on the charged materials are performed with a k-point grid and Ry Methfessel-Paxton smearing to sample the Fermi surface. Note that we are working at relatively high doping, such that the Fermi surface is large enough to be sampled correctly. Non self-consistent calculations are performed to obtain the band structure on the fine grid .

The interpolated EPC matrix elements for electron-doped WS2 are shown in Fig. 1; for electron-doped arsenene in Fig. 4; and for hole-doped phosphorene in Fig. 5. Equivalent plots for electron-doped MoS and WSe are in App. A; they are very similar to WS. The angular dependencies of the EPC are non-trivial: the EPC can undergo some rather sharp variations, most often for intravalley transitions via acoustic modes. For intervalley transitions, we observe overall smoother variations if we ignore the discontinuities coming from phonon crossings. However, intervalley scattering is activated or not depending on the valleys and the mode in rather non-trivial ways. These plots also serve to give a visual confirmation that we use sufficient sampling to capture the variations of the EPC. One important aspect to keep in mind when interpreting these plots and the transport properties of the system is that energy and momentum conservation conditions drastically reduce the final states effectively relevant for a given initial state.

The BTE is solved in all these systems, giving the mobilities shown in Fig. 6. These cover three orders of magnitude, and the hole side of phosphorene shows mobilities ten times larger than the rest, almost reaching cm/Vs at room temperature. For comparison, electron-doped graphene has mobilities of the order of to cm/Vs depending on the density. Table 1 summarizes our findings, focusing on room-temperature results, and compares them with some values available in the literature. Below, we identify the general trends and discuss three important factors in the prediction of transport properties in 2D materials. Finally, we argue that those three factors might account for a large part of the discrepancies observed in Table 1.

### iv.1 Intervalley scattering

The ranking of the above materials in terms of mobility reflects the importance of intervalley scattering. This is clearly demonstrated by considering the three TMDs: these have essentially the same type of electronic and phonon band structures, as well as similar EPC matrix elements. Yet, we obtain an order of magnitude variation in the mobilities. This stems from the position of the Q valley: indeed, as the Q valley approaches the Fermi level, it offers an additional scattering channel for the electrons. As shown in Fig. 7, the scattering time gets shorter (more scattering) at energies close to the bottom of the Q valleys.

The largest contributions to the resistivity of TMDs (1) comes from the scattering with the LA and TA modes: LA at (intravalley scattering) and M (KQ scattering) and TA at K (KK’ scattering). This is not obvious from the plots of the electron-phonon matrix elements, because the dispersions of the three acoustic modes cross each other between and K, and between and Q. However, looking at the phonon displacements, it is quite clear that in-plane acoustic modes are associated with the regions of strong electron-phonon coupling in the first three sub-plots of Fig. 1b, and in Figs. 10, 12 and 11 of the Appendix.

The multi-valley nature of a material does not necessarily deteriorate the mobility in itself. The presence of multiple valleys increases both the accessible phase space for scattered states and the density of states, and the corresponding effects on the mobility cancel each other. Rather, it is the existence of a strong EPC between the valleys that increases scattering and lowers the mobility. For the small subset studied here, all multi-valley materials showcase strong intervalley EPC. It may be argued, in general, that intervalley EPC is often strong compared to intravalley EPC. This could first be explained by the fact that intervalley EPC is not bound to vanish at long wavelengths, contrary to the coupling to acoustic phonons. Second, it involves larger phonon momenta and the EPC tends to be less screened, since the dielectric function goes to in the short-wavelength limit.

### iv.2 Symmetries of the valleys

Effective masses are among the most influential features to consider when studying mobilities, and anisotropic effective masses can be very beneficial. Indeed, small effective masses have the benefit of bringing large carrier velocities while large effective masses have the benefit of bringing high carrier densities. Thus, one might combine those benefits having a small effective mass in the transport direction and a large one in the perpendicular direction. This contributes to phosphorene’s good transport performance. In fact, phosphorene is the only material in the present study showing significant transport anisotropy. The transport properties of TMDs are isotropic because the bottom of the K and Q valleys are roughly isotropic. In the case of arsenene, each single valley is quite anisotropic, but when summing up the six equivalent valleys the angular dependency of the transport averages out and vanishesPizzi2016nc (). Fig. 8 shows the mobility of P (h), As (e) and WS (e) as a function of the direction of the electric field. The mobility of phosphorene is highly directional, while all others are isotropic.

### iv.3 Doping effects on electron-phonon interactions

The effects of doping on electron-phonon scattering are many. The first one is to move the Fermi surface within the electronic landscape. This leads to variations of the density of states and determines whether certain valleys are accessible via phonon scattering or not. As we saw in TMDs, the activation of intervalley scattering can have drastic consequences on transport. Those effects come from the energy selection rules of Eq. 4. Doping also has consequences on the strength of EPC matrix elements themselves (Eq. 5). A well-known and important effect is the additional screening coming from the electrons added in the conduction band or holes in the valence band. In our computational framework, this is inherently accounted for since we compute the linear-response of the charged system. Any EPC related to a periodic variation of the effective scalar potential in which the electrons move will be screened. This includes: i) a variation of the charge state via a variation of the area of the unit cell, as induced by longitudinal acoustic phonons; ii) any EPC related to dipole or Born effective charges, like Fröhlich or piezo-electric EPC, in which phonons interact with electrons via the generation of electric fields; iii) other less straightforward mechanisms, like the gate-induced coupling to flexural phonons in graphene Gunst2017 (); Sohier2017 (). In general, considering only the above types of EPC is largely insufficient in doped materials, as electronic screening strongly reduces their contribution, and bare interactions or weakly screened intervalley couplings dominate.

Doping also affects EPC beyond screening: as the occupations of the valleys vary, certain orbitals/bands in the material get populated or depleted, which can directly change the amplitude of the coupling. This happens in TMDs, as shown in Fig. 9 in which we computed intravalley scattering for several doping levels. In particular, we compute the average of the long-wavelength coupling along a fixed iso-energetic line at , by taking a few initial states on the iso-energetic section of each the K and Q valleys, and six phonon momenta with a fixed small norm. Linear-response calculations are then performed for each doping to capture non-trivial dependencies of the EPC matrix elements. We average on the phonon momenta and on the initial states to get for each valley . We sum the contributions from acoustic and optical phonon modes separately. Fig. 9 shows that as doping increases in WS, the intravalley coupling in K increases while it decreases in Q. The couplings to both acoustic and optical phonons show variations of . Similar trends where observed in the literatureChakraborty2012 (), but this kind of effect is poorly understood and difficult to predict, highlighting the importance of explicitly including doping in the linear-response calculations.

Last, we note that the way doping is induced can be important. An example is the gate-induced coupling to flexural phonons in graphene Gunst2017 (); Sohier2017 (), related to a broken mirror symmetry with respect to the graphene plane. However, we do not expect such effects to be significant for the subset of materials studied here.

### iv.4 Comparison with literature

Due to the diversity of the techniques employed, both at the BTE and EPC level, discrepancies between different first-principles results can have many different explanations. However, comparison with the literature points to the three factors discussed above being quite relevant. Indeed, intervalley scattering, doping and anisotropy, which we identified to be essential in determining the transport properties of 2D materials, also turn out to be the least well treated. Many models include intravalley coupling only (e.g. ”deformation potential” models), thus neglecting intervalley transitions. Furthermore, given the complexity of EPC in multivalley materials, most analytical models of EPC with fitted parameters from first-principles are incomplete. It is no coincidence that works employing such approaches show the greatest differences with the current results, often largely overestimating mobility, like the results on TMDs reported in Ref. Zhang2014, , on arsenene in Ref. Pizzi2016nc, , or phosphorene in Ref. Qiao2014, . When intervalley scattering is accounted for, it can still be a source of discrepancy depending on the relative positions of the valleys. This is the case for TMDs, in which the position of the Q valley with respect to the bottom of the K valley is difficult to determine, as mentioned before, from first-principles. For example, the energy separation between the two extrema in MoS ranges between 70Li2013 () and 260 meVLi2015 (); Kaasbjerg2012a () with the latter results closer to our calculations. This difference brings a large intervalley contribution for MoS in Ref. Li2013, that we don’t observe, and that we find instead in the W-based TMDs. The relative energy separation of the K and Q valley affects also our ranking in terms of electron mobility among the TMDs, with MoS performing better than WS, contrary to the findings in Ref. Jin2014, . In Ref. Jin2014, the energy separation between K and Q is 80 meV in MoS and 67 meV in WS resulting in a similar intervalley scattering between the two materials and a better mobility for WS on the basis of its lighter effective mass. In our case instead the K-Q energy separation shows a larger variation, from 250 meV in MoS to only 150 meV in WS resulting in an increased intervally scattering that undermines the advantage of a lighter effective mass.

In all the cases references mentioned in Table 1, first-principles calculations are done in neutral materials. Doping is sometimes entirely ignored or included only a posteriori as a shift of the Fermi level in the computation of transport properties. In some instances, analytical models of screening are used, but those are not well established in 2D materials. In any case, we have seen that the effect of doping goes beyond shifting the Fermi level and electronic screening. Thus, performing EPC calculations in neutral material might not be appropriate in general. Properly treating doping appears to be particularly important in phosphorene, where we find a larger mobility compared to most of the results in literatureJin2016 (); Rudenko2016 (); Liao2015 (); Gaddemane2018 (). In particular we find that, due to the doping induced screening, the LA mode is less effective in scattering electrons than what previously reported and only contributes for approximatively 50% of the total scattering processes. Note also that comparing transport properties is difficult when the corresponding doping is not given. Many simplified models in the literature yield mobilities that do not depend on carrier density, as shown in Table 1.

Among the materials studied here, those with isotropic valleys seem to be insensitive to the use of the closed form described in App. B to solve the BTE. Namely, we found negligible errors () on room-temperature mobilities for MoS, WS, WSe. In the case of arsenene, where the valleys are individually anisotropic but the transport properties are isotropic, we find a moderate error of . For anisotropic phosphorene, we find an error of . Although these errors can be acceptable in most cases, a full solution to the BTE beyond the relaxation-time approximation is valuable for quantitative comparison with experiments. In any case, given that the errors associated to different approximate solutions to the BTE are not obvious to determine a priori, and that our numerical solution to the full linearized BTE is comparable in terms of computational cost, we consider it a useful method in general.

## V Conclusions

We have developed an automated procedure to determine the transport properties of 2D materials from first-principles. We aim for the highest accuracy achievable within the framework of density-functional perturbation theory, with as few assumptions and simplifications as possible. The method includes several strengths and improvements with respect to existing approaches. Electron-phonon coupling matrix elements are directly computed from density functional-perturbation theory in the correct dimensionality framework and with the correct electrostatics of field-effect doping. The linearized Boltzmann transport equation is solved numerically in full, without the use of the relaxation-time approximation or any other closed-form expressions for the scattering time. The implementation of this entire transport workflow within the AiiDA infrastructure provides great flexibility to improve or adapt the method to different applications, as well as the data storage and provenance necessary to build and disseminate databases. Here, we studied in detail a small test set of 6 systems (electron-doped MoS, WS, WSe, arsenene and phosphorene as well as hole-doped phosphorene) presenting different characteristics. Our results point out the crucial role of intervalley scattering, band anisotropy and doping to the transport properties of 2D materials, in turn underscoring the importance of an accurate treatment of these aspects in first-principles simulations. Hole-doped phosphorene is found to yield the highest mobility, thanks to its mono-valley and anisotropic nature. Electron-doped arsenene shows a lower mobility than could be expected, due to the importance of intervalley scattering. The transport properties of electron-doped TMDs are found to be very sensitive to the relative positions of the K and Q valleys: these quantities are still subject of study at the highest levels of theory and experimentally. While this work is based on DFT band structures, more accurate predictions would be reached by including spin-orbit interactions, GW corrections and by fitting at least the most important features of the valleys to experimental data.

## Acknowledgements:

This work has been in part supported by NCCR MARVEL (N.M. and D.C.). Calculations were performed on the Marconi - KNL supercomputer in Cineca under PRACE project PRA15_3963. D.C. acknowledges support from the âEPFL Fellowsâ fellowship programme co-funded by Marie Sklodowska-Curie, Horizon 2020 grant agreement no. 665667. M.G. acknowledges support from the Swiss National Science Foundation through the Ambizione career programme.

## Appendix A additional EPC plots:

Here we show additional plots for the EPC matrix elements of MoS, WS, WSe and phosphorene.

## Appendix B Closed Form of BTE

Re-writing Eq. 11 as

 1 =∑k′Pkk′1−f0(k′)1−f0(k)×{τ(k)−τ(k′)v(k′)⋅uEv(k)⋅uE}, (17)

we see that a closed form can be obtained by assuming :

 1τ(k) =∑k′Pkk′1−f0(k′)1−f0(k)×{1−v(k′)⋅uEv(k)⋅uE}, (18)

This closed form expression can be further simplified by replacing with to obtain the so-called momentum relaxation-time approximation (mRTA) for the scattering time:

 1τ\footnotesizemRTA(k) =∑k′Pkk′1−f0(k′)1−f0(k)×{1−v(k′)⋅v(k)v(k)2} (19)

If the second term in brackets on the right-hand side is neglected, we finally obtain the scattering time within the energy relaxation-time approximation (eRTA):

 1τ\footnotesizeeRTA(k) =∑k′Pkk′1−f0(k′)1−f0(k) (20)

### Footnotes

1. This number is estimated by solving the BTE for each mode successively, setting the coupling to the other modes to zero. This is simply an educated estimation. Indeed, the process is not strictly valid quantitatively since the solution to the BTE including all phonon modes is not the sum of the contributions from each mode.

### References

1. D. K. Efetov and P. Kim, Controlling Electron-Phonon Interactions in Graphene at Ultrahigh Carrier Densities, Physical Review Letters 105, 256805 (2010).
2. B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Single-layer MoS2 transistors, Nature Nanotechnology 6, 147 (2011).
3. S. Lebègue, T. Björkman, M. Klintenberg, R. M. Nieminen, and O. Eriksson, Two-Dimensional Materials from Data Filtering and Ab Initio Calculations, Physical Review X 3, 031002 (2013).
4. D. Hanlon et al.Liquid exfoliation of solvent-stabilized few-layer black phosphorus for applications beyond electronics, Nature Communications 6 (2015).
5. P. Yasaei et al.High-Quality Black Phosphorus Atomic Layers by Liquid-Phase Exfoliation, Advanced Materials 27, 1887 (2015).
6. C. Gibaja et al.Few-Layer Antimonene by Liquid-Phase Exfoliation, Angewandte Chemie International Edition 55, 14345 (2016).
7. P. Ares et al.Mechanical Isolation of Highly Stable Antimonene under Ambient Conditions, Advanced Materials 28, 6332 (2016).
8. D. A. Bandurin et al.High electron mobility, quantum Hall effect and anomalous optical response in atomically thin InSe, Nature Nanotechnology 12, 223 (2016).
9. M. Ashton, J. Paul, S. B. Sinnott, and R. G. Hennig, Topology-Scaling Identification of Layered Solids and Stable Exfoliated 2D Materials, Physical Review Letters 118, 106101 (2017).
10. G. Cheon, K.-A. N. Duerloo, A. D. Sendek, C. Porter, Y. Chen, and E. J. Reed, Data Mining for New Two- and One-Dimensional Weakly Bonded Solids and Lattice-Commensurate Heterostructures, Nano Letters 17, 1915 (2017).
11. N. Mounet et al.Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds, Nature Nanotechnology 13, 246 (2018).
12. Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nature Nanotechnology 7, 699 (2012).
13. M. Chhowalla, D. Jena, and H. Zhang, Two-dimensional semiconductors for transistors, Nature Reviews Materials 1, 16052 (2016).
14. J. Ziman, Electrons and phonons: the theory of transport phenomena in solids (Clarendon, Oxford, 1960).
15. P. B. Allen, New Method for Solving Boltzmann’s equation for electrons in metals, Physical Review B 17, 3725 (1978).
16. G. Grimvall, The electron-phonon interaction in metals (North-Holland, Amsterdam, 1981).
17. F. Giustino, Electron-phonon interactions from first principles, Reviews of Modern Physics 89, 015003 (2017).
18. S. Y. Savrasov, D. Y. Savrasov, and O. K. Andersen, Linear-response calculations of electron-phonon interactions, Physical Review Letters 72, 372 (1994).
19. A. Y. Liu and A. A. Quong, Linear-response calculation of electron-phonon coupling parameters, Physical Review B 53, R7575 (1996).
20. F. Mauri, O. Zakharov, S. de Gironcoli, S. G. Louie, and M. L. Cohen, Phonon Softening and Superconductivity in Tellurium under Pressure, Physical Review Letters 77, 1151 (1996).
21. R. Bauer, A. Schmid, P. Pavone, and D. Strauch, Electron-phonon coupling in the metallic elements Al, Au, Na, and Nb: A first-principles study, Physical Review B 57, 11276 (1998).
22. S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Reviews of Modern Physics 73, 515 (2001).
23. M. M. Dacorogna, M. L. Cohen, and P. K. Lam, Self-Consistent Calculation of the q Dependence of the Electron-Phonon Coupling in Aluminum, Physical Review Letters 55, 837 (1985).
24. T. Gunst, T. Markussen, K. Stokbro, and M. Brandbyge, First-principles method for electron-phonon coupling and electron mobility: Applications to two-dimensional materials, Physical Review B 93, 035414 (2016).
25. T. Sohier, M. Calandra, and F. Mauri, Two-dimensional Fröhlich interaction in transition-metal dichalcogenide monolayers: Theoretical modeling and first-principles calculations, Physical Review B 94, 085415 (2016).
26. B. Kozinsky and N. Marzari, Static Dielectric Properties of Carbon Nanotubes from First Principles, Physical Review Letters 96, 166801 (2006).
27. T. Sohier, M. Calandra, and F. Mauri, Density-functional calculation of static screening in two-dimensional materials: The long-wavelength dielectric function of graphene, Physical Review B 91, 165428 (2015).
28. T. Sohier, M. Gibertini, M. Calandra, F. Mauri, and N. Marzari, Breakdown of Optical Phonons’ Splitting in Two-Dimensional Materials, Nano Letters 17, 3758 (2017).
29. T. Sohier, M. Calandra, and F. Mauri, Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene, Physical Review B 96, 075448 (2017).
30. Y. S. Lee, M. Buongiorno Nardelli, and N. Marzari, Band structure and quantum conductance of nanostructures from maximally localized Wannier functions: The case of functionalized carbon nanotubes, Physical Review Letters 95, 076804 (2005).
31. J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, Spectral and Fermi surface properties from Wannier interpolation, Physical Review B 75, 195121 (2007).
32. F. Giustino, M. L. Cohen, and S. G. Louie, Electron-phonon interaction using Wannier functions, Physical Review B 76, 165108 (2007).
33. M. Calandra, G. Profeta, and F. Mauri, Adiabatic and nonadiabatic phonon dispersion in a Wannier function approach, Physical Review B 82, 165111 (2010).
34. N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized Wannier functions: Theory and applications, Reviews of Modern Physics 84, 1419 (2012).
35. S. Poncé, E. Margine, C. Verdi, and F. Giustino, EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions, Computer Physics Communications 209, 116 (2016).
36. J. Sjakste, N. Vast, M. Calandra, and F. Mauri, Wannier interpolation of the electron-phonon matrix elements in polar semiconductors: Polar-optical coupling in GaAs, Physical Review B 92, 054307 (2015).
37. C. Verdi and F. Giustino, Fröhlich Electron-Phonon Vertex from First Principles, Physical Review Letters 115, 176401 (2015).
38. S. Y. Savrasov and D. Y. Savrasov, Electron-phonon interactions and related physical properties of metals from linear-response theory, Physical Review B 54, 16487 (1996).
39. O. D. Restrepo, K. Varga, and S. T. Pantelides, First-principles calculations of electron mobilities in silicon: Phonon and Coulomb scattering, Applied Physics Letters 94, 212103 (2009).
40. R. S. Shishir and D. K. Ferry, Intrinsic mobility in graphene., Journal of physics. Condensed matter : an Institute of Physics journal 21, 232204 (2009).
41. K. M. Borysenko, J. T. Mullen, E. A. Barry, S. Paul, Y. G. Semenov, J. M. Zavada, M. Buongiorno Nardelli, and K. W. Kim, First-principles analysis of electron-phonon interactions in graphene, Physical Review B 81, 121412 (2010).
42. Z. Jin, X. Li, J. T. Mullen, and K. W. Kim, Intrinsic transport properties of electrons and holes in monolayer transition-metal dichalcogenides, Physical Review B 90, 045422 (2014).
43. O. D. Restrepo, K. E. Krymowski, J. Goldberger, and W. Windl, A first principles method to simulate electron mobilities in 2D materials, New Journal of Physics 16 (2014).
44. C.-H. Park, N. Bonini, T. Sohier, G. Samsonidze, B. Kozinsky, M. Calandra, F. Mauri, and N. Marzari, Electron-Phonon Interactions and the Intrinsic Electrical Resistivity of Graphene., Nano letters 14, 1113 (2014).
45. Z. Jin, J. T. Mullen, and K. W. Kim, Highly anisotropic electronic transport properties of monolayer and bilayer phosphorene from first principles, Applied Physics Letters 109, 053108 (2016).
46. A. N. Rudenko, S. Brener, and M. I. Katsnelson, Intrinsic Charge Carrier Mobility in Single-Layer Black Phosphorus, Physical Review Letters 116, 246401 (2016).
47. Y. Trushkov and V. Perebeinos, Phonon-limited carrier mobility in monolayer black phosphorus, Physical Review B 95, 075436 (2017).
48. G. Gaddemane, W. G. Vandenberghe, M. L. Van de Put, S. Chen, S. Tiwari, E. Chen, and M. V. Fischetti, Theoretical studies of electronic transport in mono- and bi-layer phosphorene: A critical overview,  (2018).
49. K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen, Phonon-limited mobility in n-type single-layer MoS 2 from first principles, Physical Review B 85, 115317 (2012).
50. K. Kaasbjerg, K. S. Thygesen, and A.-p. P. Jauho, Acoustic phonon limited mobility in two-dimensional semiconductors: Deformation potential and piezoelectric scattering in monolayer MoS 2 from first principles, Physical Review B 87, 235312 (2013).
51. W. Li, Electrical transport limited by electron-phonon coupling from Boltzmann transport equation: An ab initio study of Si, Al, and MoS 2, Physical Review B 92, 075405 (2015).
52. M. Fiorentini and N. Bonini, Thermoelectric coefficients of -doped silicon from first principles via the solution of the Boltzmann transport equation, Physical Review B 94, 085204 (2016).
53. G. Fugallo, M. Lazzeri, L. Paulatto, and F. Mauri, Ab initio variational approach for evaluating lattice thermal conductivity, Physical Review B 88, 045430 (2013).
54. G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, and F. Mauri, Thermal Conductivity of Graphene and Graphite: Collective Excitations and Mean Free Paths, Nano Letters 14, 6109 (2014).
55. A. Cepellotti, G. Fugallo, L. Paulatto, M. Lazzeri, F. Mauri, and N. Marzari, Phonon hydrodynamics in two-dimensional materials, Nature Communications 6, 6400 (2015).
56. T. Sohier, M. Calandra, C.-H. Park, N. Bonini, N. Marzari, and F. Mauri, Phonon-limited resistivity of graphene by first-principles calculations: Electron-phonon interactions, strain-induced gauge field, and Boltzmann equation, Physical Review B 90, 125414 (2014).
57. T. Ando, Screening Effect and Impurity Scattering in Monolayer Graphene, Journal of the Physics Society of Japan 75, 074716 (2006).
58. E. H. Hwang, S. Adam, and S. Das Sarma, Carrier Transport in Two-Dimensional Graphene Layers, Physical Review Letters 98, 186806 (2007).
59. J. L. Manes, Symmetry-based approach to electron-phonon interactions in graphene, Physical Review B 76, 045430 (2007).
60. G. G. Samsonidze, E. B. Barros, R. Saito, J. Jiang, G. Dresselhaus, and M. S. Dresselhaus, Electron-phonon coupling mechanism in two-dimensional graphite and single-wall carbon nanotubes, Physical Review B 75, 155420 (2007).
61. T. Stauber, N. Peres, and F. Guinea, Electronic transport in graphene: A semiclassical approach including midgap states, Physical Review B 76, 205423 (2007).
62. S. Fratini and F. Guinea, Substrate-limited electron dynamics in graphene, Physical Review B 77, 195415 (2008).
63. E. H. Hwang and S. Das Sarma, Acoustic phonon scattering limited carrier mobility in two-dimensional extrinsic graphene, Physical Review B 77, 115449 (2008).
64. E. H. Hwang, S. D. Sarma, and S. Das Sarma, Screening-induced temperature-dependent transport in two-dimensional graphene, Physical Review B 79, 165404 (2009).
65. E. Mariani and F. von Oppen, Temperature-dependent resistivity of suspended graphene, Physical Review B 82, 195403 (2010).
66. S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Electronic transport in two-dimensional graphene, Reviews of Modern Physics 83, 407 (2011).
67. K. Kaasbjerg, K. K. Thygesen, and K. K. Jacobsen, Unraveling the acoustic electron-phonon interaction in graphene, Physical Review B 85, 165440 (2012).
68. S. Das Sarma and E. H. Hwang, Density-dependent electrical conductivity in suspended graphene: Approaching the Dirac point in transport, Physical Review B 87, 035415 (2013).
69. T. Gunst, K. Kaasbjerg, and M. Brandbyge, Flexural-Phonon Scattering Induced by Electrostatic Gating in Graphene, Physical Review Letters 118, 046601 (2017).
70. J. Qiao, X. Kong, Z.-X. Hu, F. Yang, and W. Ji, High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus, Nature Communications 5, 4475 (2014).
71. B. Liao, J. Zhou, B. Qiu, M. S. Dresselhaus, and G. Chen, Ab initio study of electron-phonon interaction in phosphorene, Physical Review B 91, 235419 (2015).
72. Y. Wang and Y. Ding, Electronic Structure and Carrier Mobilities of Arsenene and Antimonene Nanoribbons: A First-Principle Study, Nanoscale Research Letters 10, 254 (2015).
73. G. Pizzi, M. Gibertini, E. Dib, N. Marzari, G. Iannaccone, and G. Fiori, Performance of arsenene and antimonene double-gate MOSFETs from first principles, Nature Communications 7, 1258 (2016).
74. Y. Wang et al.Many-body Effect, Carrier Mobility, and Device Performance of Hexagonal Arsenene and Antimonene, Chemistry of Materials 29, 2191 (2017).
75. N. W. Ashcroft and N. D. Mermin, Solid State Physics (Brooks Cole, Belmont (USA), 1976) ISBN 0030839939.
76. B. Nag, Electron Transport in Compound Semiconductors (Springer Berlin Heidelberg, Berlin, Heidelberg, 1980) ISBN 978-3-642-81418-1.
77. O. Gritsenko, R. van Leeuwen, E. van Lenthe, and E. Baerends, Self-consistent approximation to the Kohn-Sham exchange potential, Physical Review A 51, 1944 (1995).
78. M. Kuisma, J. Ojanen, J. Enkovaara, and T. Rantala, Kohn-Sham potential with discontinuity for band gap materials, Physical Review B 82, 115106 (2010).
79. M. L. Cohen and T. K. Bergstresser, Band Structures and Pseudopotential Form Factors for Fourteen Semiconductors of the Diamond and Zinc-blende Structures, Physical Review 141, 789 (1966).
80. L.-W. Wang and A. Zunger, Local-density-derived semiempirical pseudopotentials, Physical Review B 51, 17398 (1995).
81. P. Giannozzi et al.Advanced capabilities for materials modelling with Quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).
82. P. Giannozzi et al.QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials., Journal of Physics: Condensed Matter 21, 395502 (2009).
83. C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Exact Coulomb cutoff technique for supercell calculations, Physical Review B 73, 205119 (2006).
84. S. Ismail-Beigi, Truncation of periodic image interactions for confined systems, Physical Review B 73, 233103 (2006).
85. G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, AiiDA: automated interactive infrastructure and database for computational science, Computational Materials Science 111, 218 (2016).
86. T. Brumme, M. Calandra, and F. Mauri, First-principles theory of field-effect doping in transition-metal dichalcogenides: Structural properties, electronic structure, Hall coefficient, and electrical conductivity, Physical Review B 91, 155436 (2015).
87. H. Shi, H. Pan, Y.-W. Zhang, and B. I. Yakobson, Quasiparticle band structures and optical properties of strained monolayer MoS 2 and WS 2, Physical Review B 87, 155304 (2013).
88. G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations, http://www.materialscloud.org/sssp/ (2018).
89. D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Physical Review B 88, 085117 (2013).
90. D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Physical Review B 41, 7892 (1990).
91. A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Computational Materials Science 95, 337 (2014).
92. X. Li, J. T. Mullen, Z. Jin, K. M. Borysenko, M. Buongiorno Nardelli, and K. W. Kim, Intrinsic electrical transport properties of monolayer silicene and MoS2from first principles, Physical Review B 87, 115418 (2013).
93. W. Zhang, Z. Huang, W. Zhang, and Y. Li, Two-dimensional semiconductors with possible high room temperature mobility , Nano research 7, 1731 (2014).
94. Z. Yu et al.Realization of Room-Temperature Phonon-Limited Carrier Transport in Monolayer MoS2by Dielectric and Carrier Screening, Advanced Materials 28, 547 (2016).
95. Z. Yu et al.Towards intrinsic charge transport in monolayer molybdenum disulfide by defect and interface engineering , Nature Communications 5 (2014).
96. B. Radisavljevic and A. Kis, Mobility engineering and a metal–insulator transition in monolayer MoS2, Nature Materials 12, 815 (2013).
97. S. Jo, N. Ubrig, H. Berger, A. B. Kuzmenko, and A. F. Morpurgo, Mono- and Bilayer WS2 Light-Emitting Transistors, Nano Letters 14, 2019 (2014).
98. D. Ovchinnikov, A. Allain, Y.-S. Huang, D. Dumcenco, and A. Kis, Electrical Transport Properties of Single-Layer WS2, ACS Nano 8, 8174 (2014).
99. A. S. Aji, P. Solís-Fernández, H. G. Ji, K. Fukuda, and H. Ago, High Mobility WS2 Transistors Realized by MultilayerGraphene Electrodes and Application to High Responsivity Flexible Photodetectors, Advanced Functional Materials 27 (2017).
100. Y. Cui et al.High-Performance Monolayer WS2 Field-Effect Transistors on High-k Dielectrics, Advanced Materials 27, 5230 (2015).
101. M. Z. Iqbal, M Waqasand Iqbal, M. F. Khan, M. A. Shehzad, Y. Seo, J. H. Park, C. Hwang, and J. Eom, High-mobility and air-stable single-layer WS2 field-effect transistors sandwiched between chemical vapor deposition-grown hexagonal BN films, Scientific Reports 5 (2015).
102. J.-K. Huang et al.Large-Area Synthesis of Highly Crystalline WSe2 Monolayers and Device Applications, ACS Nano 8, 923 (2014).
103. A. Allain and A. Kis, Electron and Hole Mobilities in Single-Layer WSe2, ACS Nano 8, 7180 (2014).
104. This number is estimated by solving the BTE for each mode successively, setting the coupling to the other modes to zero. This is simply an educated estimation. Indeed, the process is not strictly valid quantitatively since the solution to the BTE including all phonon modes is not the sum of the contributions from each mode.
105. B. Chakraborty, A. Bera, D. V. S. Muthu, S. Bhowmick, U. V. Waghmare, and A. K. Sood, Symmetry-dependent phonon renormalization in monolayer MoS 2 transistor, Physical Review B 85, 161403 (2012).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters