Interplay between in-plane and flexural phonons in electronic transport of two-dimensional semiconductors

TL;DR

This paper develops a theory for electron-phonon scattering in 2D semiconductors that includes anharmonic coupling between phonon modes, revealing flexural phonons are less impactful on mobility than previously thought.

Contribution

It introduces a systematic approach to account for anharmonic phonon coupling in electron transport, applied to group V 2D semiconductors, showing flexural modes have limited effect on mobility.

Findings

Flexural phonons' role is suppressed by in-plane coupling.

Mobility reduction is less than 30% at low doping levels.

Flexural phonons are less significant than in-plane phonons for transport.

Abstract

Out-of-plane vibrations are considered as the dominant factor limiting the intrinsic carrier mobility of suspended two-dimensional materials at low carrier concentrations. Anharmonic coupling between in-plane and flexural phonon modes is usually excluded from the consideration. Here we present a theory for the electron-phonon scattering, in which the anharmonic coupling between acoustic phonons is systematically taken into account. Our theory is applied to the typical group V two-dimensional semiconductors: hexagonal phosphorus, arsenic, and antimony. We find that the role of the flexural modes is essentially suppressed by their coupling with in-plane modes. At dopings lower than 10$^{12}$ cm$^{-2}$ the mobility reduction does not exceed 30\%, being almost independent of the concentration. Our findings suggest that compared to in-plane phonons, flexural phonons are considerably less important in the electronic transport of two-dimensional semiconductors, even at low carrier concentrations.

Tables

Table 1. Table 1: Relevant material constants estimated from first principles for group V hexagonal (A7 phase) 2D semiconductors: Electron (hole) effective mass m e ​ ( h ) superscript 𝑚 e h m^{\text{e}(\text{h})} (in units of free electron mass), for the case of hole doped arsenene averaged effective mass is given, electron (hole) deformation potential g e ​ ( h ) superscript 𝑔 e h g^{\text{e}(\text{h})} (in eV), Young modulus Y 𝑌 Y (in eV ⋅ ⋅ \cdot Å -2 ), Poisson ratio ν 𝜈 \nu , flexural rigidity κ 𝜅 \kappa (in eV), and a characteristic wave vector q ∗ superscript 𝑞 q^{*} at T = 𝑇 absent T=  300 K (in Å -1 ).

System	$m^{\text{e}}$	$m^{\text{h}}$	$g^{\text{e}}$	$g^{\text{h}}$	$Y$	$\nu$	$\kappa$	$q^{\ast }$
P	0.40	1.72	2.26	1.34	4.53	0.11	0.70	0.12
As	0.30	∗0.33∗	1.34	6.13	3.07	0.16	0.51	0.14
Sb	0.29	0.11	0.86	5.99	1.86	0.19	0.33	0.17

Equations

$$\sigma =\frac{e^2}{2S}\sum _{\mathbf{k}}{\tau }_{\mathbf{k}}{v}_{\mathbf{k}}^2\left(-\frac{\partial f({\epsilon }_{\mathbf{k}})}{\partial {\epsilon }_{\mathbf{k}}}\right).$$

$${\epsilon }_{\mathbf{k}}={\epsilon }_0+\frac{{\hslash }^2(\mathbf{k}-{\mathbf{k}}_0{)}^2}{2m},$$

$$\mu =T\ln \left[\exp \left(\frac{2\pi {\hslash }^{2}n}{g_{\text{v}}{g}_{\text{s}}mT}\right)-1\right],$$

$$\sigma =\frac{g_{\text{v}}{g}_{\text{s}}{e}^2\tau }{4\pi {\hslash }^2}\left\{\mu +2T\ln \left[2\cosh \left(\frac{\mu }{2T}\right)\right]\right\},$$

$$\frac{1}{\tau }=\frac{2\pi }{\hslash S}\frac{{\sum }_{\mathbf{k}{\mathbf{k}}^{\prime }}\left(-\frac{\partial f({\epsilon }_{\mathbf{k}})}{\partial {\epsilon }_{\mathbf{k}}}\right)\delta ({\epsilon }_{\mathbf{k}}-{\epsilon }_{{\mathbf{k}}^{\prime }})(1-\cos {\theta }_{\mathbf{k}{\mathbf{k}}^{\prime }})⟨|V_{\mathbf{k}{\mathbf{k}}^{\prime }}{|}^2⟩}{{\sum }_{\mathbf{k}}\left(-\frac{\partial f({\epsilon }_{\mathbf{k}})}{\partial {\epsilon }_{\mathbf{k}}}\right)},$$

$${\tau }^{-1}=\frac{m}{2\pi {\hslash }^{3}f(0)}{\int }_0^{\infty }d\epsilon \left(-\frac{\partial f(\epsilon )}{\partial \epsilon }\right){\int }_0^{2k}\frac{dq{q}^2⟨|V_q{|}^2⟩}{k^2\sqrt{k^2-q^2/4}},$$

$$V(\mathbf{r})=g{u}_{\alpha \alpha }(\mathbf{r}),$$

$$u_{\alpha \beta }(\mathbf{r})=\frac{1}{2}\left[{\partial }_{\alpha }{u}_{\beta }(\mathbf{r})+{\partial }_{\beta }{u}_{\alpha }(\mathbf{r})+{\partial }_{\alpha }h(\mathbf{r}){\partial }_{\beta }h(\mathbf{r})\right]$$

$$\begin{array}{rl}& u_{\alpha }\left(\mathbf{r}\right)\to u_{\alpha }(\mathbf{r})-\delta {\phi }_{\alpha }h(\mathbf{r})\\ & h(\mathbf{r})\to h(\mathbf{r})+\delta {\phi }_{\alpha }\left(r_{\alpha }+u_{\alpha }(\mathbf{r})\right).\end{array}$$

$$V_{\mathbf{q}}=ig{q}_{\alpha }{u}_{\alpha }(\mathbf{q})+\frac{g}{2}{f}_{\alpha \alpha }(\mathbf{q}),$$

$$H=\frac{1}{2}\int d\mathbf{r}\left\{\kappa [{\nabla }^{2}h(\mathbf{r}){]}^2+\lambda u_{\alpha \alpha }^2(\mathbf{r})+2G{u}_{\alpha \beta }^2(\mathbf{r})\right\},$$

$$⟨|V_{\mathbf{q}}{|}^2⟩=g^2\frac{T}{Y}(1-{\nu }^2)+g^2\frac{(1-\nu {)}^2}{4}⟨f^T(\mathbf{q})f^T(-\mathbf{q})⟩,$$

$$\begin{array}{rl}⟨f^T(\mathbf{q})& f^T(-\mathbf{q})⟩\simeq \Pi \left(\mathbf{q}\right)\\ & \equiv 2\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}{\left[\frac{(\mathbf{q}\times \mathbf{k}{)}^2}{q^2}\right]}^{2}G(\mathbf{k})G(\mathbf{q}-\mathbf{k}),\end{array}$$

$$\Pi (\mathbf{q})=\{\begin{array}{ll}\frac{3}{8\pi }\frac{T^2}{{\kappa }^2{q}^2}& \text{, }q\gg q^{\ast }\\ C{\left(\frac{T}{\kappa }\right)}^{2-\eta }{\left(\frac{\kappa }{Y}\right)}^{\eta }{q}^{-2+2\eta }& \text{, }q\ll q^{\ast }\end{array},$$

$$⟨f^T(\mathbf{q})f^T(-\mathbf{q})⟩\simeq \frac{\Pi (\mathbf{q})}{1+\frac{Y}{4T}\Pi (\mathbf{q})}.$$

$$⟨f^T(\mathbf{q})f^T(-\mathbf{q})⟩=\{\begin{array}{ll}\frac{3}{8\pi }\frac{T^2}{{\kappa }^2{q}^2}& \text{, }q\gg q^{\ast }\\ \frac{4T}{Y}-\frac{16}{C}\frac{T^2}{{\kappa }^2}{\left(\frac{TY}{{\kappa }^2}\right)}^{-2+\eta }{q}^{2-2\eta }& \text{, }q\ll q^{\ast }\end{array}$$

$$E_{\mathrm{elas}}=\frac{\kappa }{2}\int d\mathbf{r}({\nabla }^{2}h{)}^2.$$

$$\Delta E_{\mathrm{tot}}=\frac{\kappa }{4}\frac{1}{R^2},$$

$$E_{\mathrm{elas}}=\frac{Y\nu }{2(1-{\nu }^2)}{u}_{\alpha \alpha }^2+\frac{Y}{2(1+\nu )}{u}_{\alpha \beta }^2.$$

$$E_{\mathrm{elas}}=\frac{1}{2}Yu_{\mathrm{\alpha\alpha}}^{2}\bigg{|}_{\sigma_{\beta\beta}=0},$$

$$\nu=-\frac{\Delta l_{\beta\beta}}{\Delta l_{\alpha\alpha}}\bigg{|}_{\sigma_{\beta\beta}=0},$$

$$E_{\mathrm{int}}=2g^{\mathrm{e}(\mathrm{h})}u=\{\begin{array}{l}\Delta E^{\mathrm{CBM}}\text{(electrons)}\\ \Delta E^{\mathrm{VBM}}\text{(holes).}\end{array}$$

$$⟨|V_{\mathbf{q}}{|}^2⟩=Z^{-1}\int Dh\left(\mathbf{r}\right)\int D\mathbf{u}\left(\mathbf{r}\right){\left|V_{\mathbf{q}}\right|}^2{e}^{-\beta H},$$

$$Z=\int Dh\left(\mathbf{r}\right)\int D\mathbf{u}\left(\mathbf{r}\right)e^{-\beta H}.$$

$$\begin{array}{r}{f}_{\alpha \beta }\left(\mathbf{q}\right)=i\left[q_{\alpha }{\varphi }_{\beta }\left(\mathbf{q}\right)+q_{\beta }{\varphi }_{\alpha }\left(\mathbf{q}\right)\right]+P_{\alpha \beta }^T\left(\mathbf{q}\right)f^T\left(\mathbf{q}\right).\end{array}$$

$${\tilde{u}}_{\alpha }\left(\mathbf{q}\right)\equiv u_{\alpha }\left(\mathbf{q}\right)+{\varphi }_{\alpha }\left(\mathbf{q}\right)-\nu \frac{i}{2}\frac{q_{\alpha }}{q^2}{f}^T\left(\mathbf{q}\right),$$

$$\begin{split}&H=\sum_{{\bf q}}\bigg{\{}\frac{1}{2}\kappa q^{4}\left|h_{{\bf q}}\right|^{2}+\frac{Y}{8}f^{T}\left({\bf q}\right)f^{T}\left(-{\bf q}\right)\\ &+\frac{1}{2}\left[(\lambda+G)q_{\alpha}q_{\beta}+G\,q^{2}\delta_{\alpha\beta}\right]\tilde{u}_{\alpha}\left({\bf q}\right)\tilde{u}_{\beta}\left(-{\bf q}\right)\bigg{\}},\\ \end{split}$$

$$V_{\mathbf{q}}=ig{q}_{\alpha }{\tilde{u}}_{\alpha }\left(\mathbf{q}\right)+\frac{g}{2}\left(1-\nu \right)f^T\left(\mathbf{q}\right).$$

$$\begin{array}{rl}⟨{\left|V_{\mathbf{q}}\right|}^2⟩& =g^2{q}_{\alpha }{q}_{\beta }⟨{\tilde{u}}_{\alpha }\left(\mathbf{q}\right){\tilde{u}}_{\beta }\left(-\mathbf{q}\right)⟩\\ & +g^2\frac{(1-\nu {)}^2}{4}⟨f^T\left(\mathbf{q}\right)f^T\left(-\mathbf{q}\right)⟩.\end{array}$$

$$⟨{\tilde{u}}_{\alpha }\left(\mathbf{q}\right){\tilde{u}}_{\beta }\left(-\mathbf{q}\right)⟩=\frac{T}{\left(\lambda +2G\right)q^2}{P}_{\alpha \beta }^L\left(\mathbf{q}\right)+\frac{T}{G{q}^2}{P}_{\alpha \beta }^T\left(\mathbf{q}\right),$$

Full text

††preprint: APS/123-QED

I Introduction

Charge carrier transport in two-dimensional (2D) materials is different from their three-dimensional analogs Katsnelson-Book ; Avouris . Among other reasons Peres ; DasSarma , this is mainly related to the existence of the flexural phonon modes (out-of-plane vibrations), whose energy dispersion is quadratic in wave vector compared to the linear dispersion of conventional (in-plane) modes Katsnelson-Book ; Lifshitz ; AccChemRes . As the charge carrier scattering on phonons is the main factor limiting the intrinsic mobility, understanding of the interplay between flexural and in-plane phonons is of vital importance for 2D material’s electronic transport.

The contributions of in-plane and flexural phonons are usually considered independently of each other, at the level of the harmonic approximation. Moreover, the role of flexural phonons is often completely neglected Kaasbjerg2013 ; Qiao ; InSe , assuming that their contribution is suppressed by the presence of a substrate Amorim . At the same time, flexural phonons are considered to play the dominant role in the electron-phonon scattering of free-standing graphene Morozov2008 ; Mariani2008 ; Castro2010b , and may contribute significantly in the presence of a dielectric substrate Gunst ; Sohier17 . Transport properties of conventional 2D semiconductors like black phosphorus are less affected by flexural modes in the regime of high doping Rudenko2016 , whereas at low dopings the situation is unclear. Furthermore, flexural modes could significantly renormalize the electronic spectrum, leading to the formation of bound states Flexuron ; Brener .

The separation between the in-plane and flexural modes becomes undefined at wave vectors $q\lesssim q^{*}$, smaller than the characteristic vector $q^{*}=\sqrt{3TY/16\pi\kappa^{2}}$, determined by the Young modulus $Y$, flexural rigidigy $\kappa$, and temperature $T$ AccChemRes ; Nelson-orig ; Nelson . This is a manifestation of the anharmonic coupling regime, which is typical to low carrier concentrations with small Fermi wave vectors, $k_{\text{F}}\lesssim q^{*}$. In this situation, the electron-phonon interaction for flexural modes diverges and is usually cut off. Several attempts have been made to systematically account for the anharmonic effects in the transport properties of graphene Mariani2010 ; Castro2010b ; Ochoa2011 ; Gorniy . However, only little attention to this problem has been given in the context of conventional semiconductors Fischetti . Despite the availability of advanced ab initio computational techniques developed to describe electron-phonon scattering Marzari ; Giustino ; Ponce ; Sohier , their applicability is limited with respect to the above-mentioned effects.

Group V elemental semiconductors is an emerging class of 2D materials with attractive electronic properties Sofer ; Ares2018 . Hexagonal (A7) 2D phases of phosphorus (P), arsenic (As), and antimony (Sb) are among the most recent materials fabricated experimentally blueP-syn ; Sb-syn ; As-syn . Apart from being mechanically stable blueP ; AsSb_stability , they are proposed to have excellent electron transport characteristics Wang ; AsSb_mobility and a high degree of gate tunability Ghosh ; Shu ; Prishchenko due to the buckled crystal structure. Mechanically, these materials are highly flexible with the flexural rigidities in the range 0.3–0.8 eV blueP ; Lugovskoi , which is significantly smaller compared to graphene kappaC . This property allows us to expect a considerable anharmonic coupling in these materials, and its possible influence on the electron transport properties.

In this paper, we focus on a theory for phonon-mediated charge carrier scattering in 2D semiconductors taking the coupling between the acoustic phonons into account. We use the diagrammatic approach to calculate the anharmonic contribution to the electron-phonon scattering rate, ending up at two plausible approximations. We apply the developed theory to calculate electron and hole mobilities in three representative group V isotropic semiconductors: hexagonal monolayers P (also known as blue phosphorene), As (arsenene), and Sb (antimonene). To this end, we estimate relevant model parameters from first-principles density functional theory (DFT) calculations. We find that although flexural phonons tend to reduce the mobility at concentrations below 1014 cm*-2* by 15–30%, their effect reaches its maximum at around 1012 cm*-2*. At lower dopings, the mobility becomes essentially independent of the concentration, being the manifestation of the anharmonic coupling between the flexural and in-plane phonon modes.

The paper is organized as follows. In Sec. II, we present a theory of charge carrier scattering in 2D semiconductors, explicitly considering in-plane and out-of-plane phonons, as well as their anharmonic coupling. In Sec. III, we use first-principles calculations to estimate model parameters for hexagonal P, As, and Sb monolayers. The results of numerical calculations are presented in Sec. IV, where we investigate the interplay between in-plane and flexural phonons in the context of charge carrier scattering, as well as estimate upper limits for the carrier mobility in the materials under consideration. In Sec. V, we summarize our results, and conclude the paper.

II Theory

II.1 Carrier mobility and relaxation time

Carrier mobility of isotropic electron gas can be expressed as $\mu_{c}=\sigma/ne$, where $n$ is the carrier concentration and $\sigma$ is the dc conductivity, which has the form Ziman

[TABLE]

Here, $S$ is the sample area, $\varepsilon_{\bf k}$ is the band energy, $v_{\bf k}$ is the group velocity, $\tau_{\bf k}$ is the momentum-dependent scattering relaxation time, and $f(\varepsilon_{\bf k})=\{1+\mathrm{exp}[(\varepsilon_{\bf k}-\mu)/T])\}^{-1}$ is the Fermi-Dirac distribution function. The energy dispersion of charge carriers is assumed to have the following form near the band edge:

[TABLE]

where $m$ is the electron ($m^{\mathrm{e}}$) or hole ($m^{\mathrm{h}}$) effective mass and $\varepsilon_{0}$ and ${\bf k}_{0}$ are the energy and wave vector determining the position of band edges, respectively. The chemical potential $\mu$ can be determined from the carrier concentration $n$ as

[TABLE]

where $g_{\text{v}}g_{\text{s}}$ is the valley-times-spin-degeneracy factor. In turn, one can recast Eq. (1) in the form,

[TABLE]

which in the limit $T\rightarrow 0$ transforms into the standard conductivity formula $\sigma=ne^{2}\tau/m$.

Here we assume the relaxation time to be averaged over thermally broadened Fermi surface states $\tau\equiv\langle\langle\tau_{\bf k}\rangle\rangle_{\mathrm{FS}}$, which is dependent on the chemical potential $\mu$ and temperature $T$.

The corresponding expression for $\tau$ can be obtained from the variational solution of the isotropic Boltzmann transport equation Ziman , yielding

[TABLE]

where ${\bf k}={\bf k}^{\prime}+{\bf q}$ with ${\bf q}$ being the scattering wave vector, $V_{{\bf k}{\bf k}^{\prime}}$ is the effective electron-phonon scattering potential defined in Sec. II.2, and $\langle...\rangle=\mathrm{Tr}(e^{-\beta H}...)/\mathrm{Tr}(e^{-\beta H})$ with $\beta=1/T$ denotes thermal averaging over the phonon states with the Hamiltonian $H$ defined below. In Eq. (5), emission and absorption of phonons are not explicitly involved in the energy conservation law, which assumes elastic scattering, i.e., $\hbar\omega_{k_{\text{F}}}\ll\varepsilon_{F}$ with $\omega_{k_{\text{F}}}$ being the phonon frequency at the Fermi wave vector. This condition is satisfied for acoustic phonons considered in our study. After some algebra, Eq. (5) can be rewritten in a more convenient form as follows:

[TABLE]

where $k=\sqrt{2m\varepsilon}/\hbar$ plays the role of the Fermi wave vector at $\mu\gg T$.

We note that Eq. (6) is only valid for conventional 2D semiconductors with quadratic dispersion of charge carriers. Deviations from the quadratic dispersion, as well as the presence of a dielectric environment (e.g., substrates), could result in a considerable modification of the energy dependence of the scattering rate (see, e.g., Refs. Hwang, ; Kaasbjerg12, ; Min, ).

II.2 Electron-phonon scattering

Here we first define the scattering potential. We restrict ourselves to the long-wavelength limit and take both in-plane and out-of-plane (flexural) phonon modes into account. We assume that electrons interact with phonons through the deformation potential of the form

[TABLE]

where

[TABLE]

is the strain tensor and $g$ is the coupling constant, whereas $u_{\alpha}({\bf r})$ and $h({\bf r})$ are in-plane and out-of-plane displacement fields, respectively.

The form of the deformation potential in Eq. (7) is justified by the following considerations. In the absence of substrates or external tension, the electronic structure is invariant under rotations of the crystal as a whole in three-dimensional space guinea_collective_2014 . Therefore, $V\left({\bf r}\right)$ must be invariant under the following transformations:

[TABLE]

To lowest order in $\partial_{\alpha}u_{\beta}\left({\bf r}\right)$ and $\partial_{\alpha}h\left({\bf r}\right)$, this implies that $V({\bf r})$ can only depend on the components of the strain tensor, which are scalar quantities under the transformation in Eq. (9) footnote . In principle, $V({\bf r})$ can depend on all components of $u_{\alpha\beta}\left({\bf r}\right)$. Here we assume the isotropic form in Eq. (7) for simplicity. Notice that, as a consequence of rotational invariance, in-plane and out-of-plane coupling constants in Eq. (7) are not independent SanJose ; guinea_collective_2014 .

In view of the absence of horizontal mirror ($\sigma_{\mathrm{h}}$) symmetry in buckled hexagonal semiconductors ($D_{3\mathrm{d}}$ point group) being studied in this work, Eqs. (7) and (8) require additional justification. Particularly, as has been shown in Ref. Fischetti, , single-phonon processes associated with flexural modes are highly relevant for non-mirror-symmetric 2D Dirac materials like silicene and germanene. Such processes correspond to additional terms in Eq. (7) that are linear in out-of-plane displacements, i.e., $\nabla h$. However, in our case those terms are irrelevant for the following reasons. Here we are focused on lightly doped semiconductors with quadratic energy dispersion represented by Eq. (2). Moreover, the intervalley scattering can be neglected because $k_{F}\ll k_{0}$. Under these assumptions, it is obvious from in-plane symmetry that $\nabla h$ terms can only enter the interacting potential via the transformation ${\bf k}_{0}\rightarrow{\bf k}_{0}+\alpha\nabla h$. However, such terms are forbidden by rotational symmetry in three-dimensional space. This consideration does not forbid the appearance of higher-gradient terms in the deformation potential. As can be shown, $\nabla^{2}h$ terms are forbidden, too, which follows from the combination of inversion and time-reversal symmetry, particularly because ${\bf k}_{0}\rightarrow-{\bf k}_{0}$. On the other hand, $\nabla^{3}h$ and higher terms acquire additional power of phonon wave vector ${\bf q}$ in comparison to $\nabla u$ terms, and since $|{\bf q}|\lesssim 2k_{F}$, such terms have smallness in electron concentration and, therefore, can be neglected.

Let us now evaluate the scattering probability $\langle|V_{\bf q}|^{2}\rangle$. The Fourier transform of the scattering potential $V({\bf r})$ is then given by

[TABLE]

where $u_{\alpha}({\bf q})$ are the Fourier components of the in-plane displacement field $u_{\alpha}({\bf r})$ and $f_{\alpha\beta}({\bf q})=-S^{-\frac{1}{2}}\sum_{\bf k}k_{\alpha}(q_{\beta}-k_{\beta})h_{\bf k}h_{{\bf q}-{\bf k}}$ is the Fourier transform of $f_{\alpha\beta}({\bf r})=\partial_{\alpha}h({\bf r})\partial_{\beta}h({\bf r})$. In the long-wavelength limit, fluctuations of the displacement fields are described by the Hamiltonian Nelson ; Nelson-orig ; le_doussal_self-consistent_1992

[TABLE]

where $\lambda$ and $G$ are the Lamé parameters. We assume that fluctuations of the displacement fields can be described classically. Neglecting quantum effects is fully justified at not-too-small temperatures, when $\hbar\omega(k_{F})\ll T$. In the calculation of the thermal average $\langle\left|V_{{\bf q}}\right|^{2}\rangle$, we can take advantage of the fact that the Hamiltonian in Eq. (11) is quadratic with respect to the in-plane displacement fields $u_{\alpha}\left({\bf r}\right)$, which allows one to integrate them out exactly Nelson ; Nelson-orig ; le_doussal_self-consistent_1992 . After some algebraic manipulations Nelson ; Nelson-orig , we obtain the following expression for the scattering matrix (see the Appendix for details):

[TABLE]

where $\nu=\lambda/(\lambda+2G)$ is the 2D Poisson ratio and $Y=\lambda(1-\nu^{2})/\nu$.

Being focused on the processes with small momentum transfer ${\bf q}={\bf k}-{\bf k}^{\prime}$, we ignored the effect of the wave-function overlap between the initial ${\bf k}$ and final ${\bf k}^{\prime}$ states, which might lead to some overestimation of the scattering probability.

In Eq. (12), the first term coincides with the usual contribution of in-plane phonons in the harmonic approximation, whereas the second term describes both out-of-plane and cross interactions generated by anharmonic coupling; $f^{T}\left({\bf q}\right)\equiv\left(\delta_{\alpha\beta}-q_{\alpha}q_{\beta}/q^{2}\right)f_{\alpha\beta}\left({\bf q}\right)$ denotes the transverse part of $f_{\alpha\beta}\left({\bf q}\right)$. In what follows, we use two different approximations to calculate the correlation function $\langle f^{T}\left({\bf q}\right)f^{T}\left(-{\bf q}\right)\rangle$.

Wick decoupling approximation. In the first approximation, we assume the validity of the Wick theorem, which yields $\langle h_{\bf k}h_{{\bf q}-{\bf k}}h_{{\bf k}^{\prime}}h_{-{\bf q}-{\bf k}^{\prime}}\rangle\simeq\left(\delta_{{\bf k},{{\bf q}+{\bf k}^{\prime}}}+\delta_{{\bf k},-{\bf k}^{\prime}}\right)\langle h_{\bf k}h_{-\bf k}\rangle\langle h_{{\bf k}-{\bf q}}h_{{\bf q}-{\bf k}}\rangle$ Katsnelson-Book . In terms of diagrams, this assumption corresponds to the neglect of vertex contributions to the two-particle Green’s function. In this approximation, the desired correlation function reads:

[TABLE]

where $G({\bf q})\equiv\langle h_{\bf q}h_{-\bf q}\rangle$ is the single-particle correlation function of out-of-plane fields. The function $\Pi\left({\bf q}\right)$ exhibits a crossover from a harmonic behavior for $q\gg q^{*}$ to an anomalous scaling behavior for $q\ll q^{*}$, which is controlled by anharmonic interactions. In this work, we addressed the calculation of $\Pi\left({\bf q}\right)$ within the self-consistent screening approximation (SCSA) le_doussal_self-consistent_1992 ; le_doussal_anomalous_2018 . In the limit of large and small wave vectors, $\Pi\left({\bf q}\right)$ behaves asymptotically as:

[TABLE]

where $\eta$ is an universal exponent and $C$ is a dimensionless constant. The value of $\eta$ is close to $0.821$ within the SCSA, in the physical case of two-dimensional membranes embedded in three-dimensional space. Through a numerical solution of the SCSA equations le_doussal_self-consistent_1992 ; le_doussal_anomalous_2018 , we determined $\Pi\left({\bf q}\right)$ for arbitrary values of the wave vector, including the intermediate scale $q\simeq q^{*}$ and we obtained the value $C\simeq 2.60$ for the dimensionless amplitude in Eq. (14).

Screening approximation. In the second approximation, we consider the following expression for the correlation function:

[TABLE]

Equation (15) can be derived diagrammatically by summing a geometric series of polarization bubble diagrams connected by bare interaction lines, the same diagrams which determine, within the SCSA, the screened interaction between flexural phonons le_doussal_self-consistent_1992 ; le_doussal_anomalous_2018 . An expression equivalent to Eq. (15) is reported in Ref. guinea_collective_2014 .

Adopting Eq. (15), one arrives at the following asymptotic limits for the correlation function:

[TABLE]

As expected, the two approximations result in the same behavior at $q\gg q^{*}$, while the result at $q\ll q^{*}$ is different. In the Wick decoupling approximation, the electron-phonon coupling $\langle|V_{\bf q}|^{2}\rangle$ scales as $1/q^{2-2\eta}$ for small $q$ and it diverges in the limit $q\rightarrow 0$, although the divergence is slow, because the exponent $\eta$ is not very different from $1$. By contrast, in the screening approximation $\langle|V_{\bf q}|^{2}\rangle$ approaches the constant value $2g^{2}\frac{T}{Y}(1-\nu)$ for $q\rightarrow 0$. This expression ensures that, in the screening approximation, scattering on in-plane phonons always gives a larger contribution than scattering on out-of-plane phonons. In the discussion below, the two approximations will be considered as the upper and lower limits for the relaxation time.

The expressions presented above allow one to estimate the dc conductivity (or mobility) of isotropic 2D electron gas in the presence of elastic scattering on acoustic phonons. Importantly, this approach is equally valid for small carrier concentrations, where scattering involves the coupling between phonons modes.

III Estimation of material constants

III.1 Calculation details

To estimate material constants, we use first-principles DFT calculations. All calculations have been carried out using the projected augmented-wave Blochl formalism as implemented in the Vienna ab initio simulation package (vasp) Kresse1996a ; Kresse1996b ; Kresse1999 . To describe exchange-correlation effects, we employed the gradient-corrected approximation in the parametrization of Perdew-Burke-Ernzerhof pbe . An energy cutoff of 500 eV for the plane waves and the convergence threshold of 10*-9* eV were used in all cases. To avoid spurious interactions between the cells, a vacuum slab of 50 Å  was added in the direction perpendicular to the 2D sheet. Structural relaxation including the optimization of in-plane lattice constants was performed until the forces acting on atoms were less than 10*-3* eV/Å. Unlike the cases of P and As monolayers, spin-orbit interaction in monolayer Sb is important Rudenko2017 and, therefore, has been taken into account.

Hexagonal monolayers studied in this work adopt honeycomb structure with vertically displaced (buckled) sublattices. We obtained the following optimized values of the lattice constant ($a$) and buckling parameter: 3.28 Å and 1.24 Å for P, 3.61 Å and 1.40 Å for As, 4.12 Å and 1.64 Å for Sb, respectively. Primitive hexagonal cells were used in most of the calculations, for which a Monkhorst-Pack kmesh $\Gamma$-centered (32$\times$32) mesh was adopted to sample the Brillouin zone. To calculate the Poisson ratio and Young modulus, a rectangular $(\sqrt{3}/2\times 1)a$ unit cell was used, together with a $(32\times 24)$ k-point mesh. To induce flexural deformations needed for the calculation of flexural rigidities, we considered rectangular (6$\sqrt{3}\times$1)$a$ supercells, such that the deformation period amounts to $l=6\sqrt{3}a\sim 40$ Å  for each system considered. We first apply a field of 1D sinusoidal out-of-plane deformations in the form $h(x,y)=h\,\mathrm{sin}(qx)$ to each system, where $h$ the deformation amplitude and $q=2\pi/l$. We then perform full structural relaxation for a series of fixed amplitudes $h$, which allows us to obtain the elastic energies. In the calculations with rectangular supercells, a ($4\times 24$) k-point mesh was used.

III.2 Effective masses

Electron and hole effective masses are estimated from the band structures by fitting the conduction band minimum (CBM) and valence band maximum (VBM), respectively, using the expression given by Eq. (2). The calculated band structure used for the fitting, and the corresponding effective masses are given in Fig. 1. In all three systems considered, CBM resides at the $\Sigma$ point along the $\Gamma$–M high symmetry line. The corresponding bands are spin ($g_{\text{s}}=2$) and valley degenerate ($g_{\text{v}}=2$). The constant-energy (Fermi) contours at small energies form elliptical pockets Lugovskoi , meaning that the electron effective mass is different in the direction perpendicular to $\Gamma$–M, denoted in Fig. 1 as $-$PΣ–PΣ. We denote the corresponding masses as longitudinal ($m_{\text{l}}^{\text{e}}$) and transverse ($m^{\text{e}}_{\text{t}}$), as shown in Fig. 1. In practical calculations, we use geometrically averaged effective masses, $m^{\text{e}}=\sqrt{m^{\text{e}}_{\text{l}}m^{\text{e}}_{\text{t}}}$, with the resulting values listed in Table 1. The obtained values overall agree with the literature data AsSb_mobility2015 ; AsSb_mobility ; AsSb_mobility2017 .

The case of holes in monolayer P is similar to the case of electrons. Since VBM is located at the $\Sigma^{\text{h}}$ point along $\Gamma$–K symmetry line, the effective mass tensor is anisotropic, and $m^{\text{h}}$ can be calculated as the average of its longitudinal and transverse components. The situation with the hole effective masses in As and Sb is different. In these systems, VBM is located at the $\Gamma$ point, making the Fermi contour isotropic with the direction-independent hole effective masses. Unlike monolayer Sb, where VBM at the $\Gamma$ point is doubly degenerate ($g_{\text{s}}=2$), VBM of monolayer As is four times degenerate due to the absence of strong spin-orbit coupling. As one can see, there are two different bands forming the valence band in monolayer As, meaning that there are two transport channels for holes in As. To take them into account, specifically for As we define the transport effective mass, which is the sum of light and heavy hole masses, $m^{\text{h}}_{\text{tr}}=m^{\text{h}}_{\text{lt}}+m^{\text{h}}_{\text{h}}$. On the other hand, the density of states effective mass is different, $m^{\text{h}}_{\text{d}}=\sqrt{m^{\text{h}}_{\text{h}}m^{\text{h}}_{\text{lt}}}$. Given that the conventional electron-phonon contribution to the mobility can be expressed through the product $m_{\text{d}}m_{\text{tr}}$ Takagi , we employ $m^{\text{h}}=\sqrt{m^{\text{h}}_{\text{d}}m^{\text{h}}_{\text{tr}}}$ in the calculations of As. In all other systems, $m_{\text{d}}$ and $m_{\text{tr}}$ are equivalent.

III.3 Elastic properties

III.3.1 Flexural rigidities

Flexural rigidities $\kappa$ can be estimated from first principles starting from a macroscopic expression for the elastic energy density of a 2D membrane, Eq. (11). In the absence of strain, it reads

[TABLE]

The elastic energy $E_{\mathrm{elas}}$ can be directly obtained from DFT calculations as the total energy difference $\Delta E_{\mathrm{tot}}$ between the corrugated and flat states of a membrane, $E_{\mathrm{elas}}=\Delta E_{\mathrm{tot}}$. By employing the sinusoidal height field in the form $h(x,y)=h\,\mathrm{sin}(qx)$, the constant $\kappa$ can be estimated by a regression fit from

[TABLE]

where $R$ is the curvature radius of a sine wave at its extrema, given by $1/R=hq^{2}$.

In Fig. 2(b), we show the calculated dependencies $\Delta E_{\mathrm{tot}}$ as a function of $hq^{2}$, which allow us to estimate $\kappa$ for the systems under consideration (see Table 1 for summary). In all cases the expression given by Eq. (17) perfectly fits DFT data.

III.3.2 In-plane elastic constants

In the absence of flexural deformations, the elastic energy density $E_{\mathrm{elas}}$ of a hexagonal strained 2D material can be expressed from Eq. (11) as

[TABLE]

The Young modulus $Y$ can be found straightforwardly from DFT calculations by applying a uniaxial deformation in the desired direction and allowing full relaxation in other directions:

[TABLE]

where $\sigma_{\beta\beta}$ is the component of the stress tensor in the relaxed direction. Similarly, the Poisson ratio can be obtained from the same set of calculations:

[TABLE]

where $\Delta l_{\alpha\alpha}$ and $\Delta l_{\beta\beta}$ are the change of lattice constants in the strained and relaxed directions, respectively.

In Fig. 2(a), we show the calculated $\Delta E_{\mathrm{tot}}$ as a function of uniaxial strain, which perfectly fits the macroscopic energy expression given above. This allows us to readily estimate the Young modulus $Y$ for the materials under consideration. The calculated Poisson ratios are listed in Table 1.

III.4 Deformation potentials

III.4.1 Interaction with in plane phonons

To estimate the coupling of charge carriers with phonons from DFT calculations, we adopt the concept of the deformation potentials BirPikus . Assuming the long-wavelength limit, we quantify linear response of the band edges to external deformation. To this end, we apply both tensile and compressive biaxial strain and obtain relative band shifts $\Delta E^{\mathrm{VBM}}$ and $\Delta E^{\mathrm{CBM}}$ as shown in Figs. 2(b) and 2(c) as a function of strain $u=u_{\alpha\alpha}=u_{\beta\beta}$. The band shifts are calculated relative to the vacuum level $E_{\text{vac}}$, determined for every deformed state, that is $\Delta E=E-E_{\text{vac}}$, where $E$ is the energy of the corresponding band edge.

The interaction energy $E_{\mathrm{int}}$ is linear in strain and takes the form

[TABLE]

Since the sign of $E_{\mathrm{int}}$ depends on the strain type (tensile or compressive), the sign of $g^{\text{e}(\text{h})}$ is the matter of convention. The coupling constants are given in Fig. 2 and also summarized in Table 1. While the effective masses, as well as hole deformation potentials ($g^{\text{h}}$) for monolayers As and Sb agree well with the literature data, larger electron deformation potentials ($g^{\text{e}}$) have been reported for these materials in Refs. AsSb_mobility2015 ; AsSb_mobility ; AsSb_mobility2017 . The discrepancy can be attributed to a different approach utilized in those works. Particularly, uniaxial strain was used to quantify the band shifts, while we utilized biaxial strain in our work. The latter approach implies that only diagonal elements of the deformation potential tensor are taken into account in our study [see Eq. (7)]. This approximation is expected to quantitatively overestimate the mobilities calculated below.

IV Results and discussion

We first analyze the role of flexural phonons and their anharmonic coupling in the scattering of charge carriers in 2D semiconductors P, As, and Sb. To this end, we decompose the total scattering rate into the two contributions $\tau^{-1}=\tau_{\parallel}^{-1}+\tau_{\perp}^{-1}$, each of which describes the scattering on pure in-plane modes [first term on the right-hand side of Eq. (12)] and on flexural modes subject to anharmonic coupling with in-plane modes [second term in the same equation]. In Fig. 3, we show the dependence $\tau^{-1}_{\perp}/\tau^{-1}$ of the carrier concentration both for electron and hole dopings calculated at $T=300$ K. As long as the effective mass approximation holds, from Eq. (6) we have $\tau_{\parallel}^{-1}=mg^{2}T(1-\nu^{2})/\hbar^{3}Y$, meaning that the doping dependence arises from the $\tau_{\perp}^{-1}$ term only. In the limit of large concentrations $\tau_{\perp}^{-1}/\tau^{-1}\rightarrow 0$ because $\tau_{\perp}^{-1}\sim 1/\varepsilon_{F}$, indicating that flexural phonons are negligible in this regime, which can be clearly seen from Fig. 3 for all the cases considered. At small concentration the behavior is quantitatively different for the two approximations describing the anharmonic coupling. Nevertheless, the trend is qualitatively similar: The contribution of flexural phonons becomes independent or weakly dependent on the doping. This behavior is consistent with ${\bf q}$-dependence of the scattering matrix $\langle|V_{\bf q}|^{2}\rangle$, which at $q\ll q^{*}$ behaves as \sim$$q^{-0.36} within the Wick decoupling approximation and approaches a constant in the screening approximation. On average, the contribution of flexual phonons at room temperature does not exceed 30% of the total scattering rate for all materials considered.

Let us now examine the temperature dependence of the scattering rate. As has been shown above $\tau_{\parallel}^{-1}\sim T$ independently of the carrier concentration. On the other hand, $\tau_{\perp}^{-1}\sim T^{2}$ in the harmonic regime ($q\gg q^{*}$), while it softens significantly in the strong anharmonic regime ($q\ll q^{*}$), yielding $\tau_{\perp}^{-1}\sim T^{1.18}$ or $\sim T$ depending on the approximation employed. The temperature dependence calculated for materials under consideration is shown in Fig. 4. The trend is similar for all the cases considered. In particular, at low temperatures $\tau^{-1}_{\perp}/\tau^{-1}_{\parallel}$ is small as the electron gas is degenerate ($T\ll\varepsilon_{F}$), resulting in a negligible flexural phonon contribution. At higher temperatures, the behavior essentially depends on the carrier concentration. At small concentration around $10^{11}$ cm*-2*, $\tau^{-1}_{\perp}/\tau^{-1}_{\parallel}$ is almost independent of the temperature, indicating the strongly anharmonic regime. At concentrations of the order of $10^{13}$ cm*-2*, the corresponding ratio $\sim T$, as is expected from the limit $q\gg q^{*}$.

As a final step, we make a quantitative estimate of the carrier mobility, presented in Fig. 5 for a wide range of electron and hole doping in P, As, and Sb. At large carrier concentrations, scattering by in-plane phonons dominates, and the mobility approaches its upper limit for a given temperature. As the concentration drops down, the mobility decreases by 15–30% for different materials in the carrier density range from 1012 to 1014 cm*-2*. This effect is attributed to the increased role of flexural phonons, whose scattering rate is inversely proportional to the carrier concentration. At lower dopings, the scattering behavior is modified by the strong anharmonic coupling between in-plane and flexural phonons. In particular, below 1012 cm*-2* the mobility behavior is nearly independent of the material and doping type, resulting in a constant value.

Quantitatively, the mobilities are strongly dependent on the carrier type. The electron mobilities are significantly higher for all the systems considered, which can be attributed to relatively small effective masses and deformation potentials, listed in Table 1. The lowest values not exceeding 1000 cm2V*-1s-1* are obtained for the hole doping of P and As, i.e., materials with relatively massive carriers. It is interesting to note that the electron-doped hexagonal (blue) phosphorus exhibits mobilities an order of magnitude larger than its orthorhombic (black) counterpart Rudenko2016 . It should be emphasized that the results presented in Fig. 5 represent an upper limit for the intrinsic mobility. A number of important factors, such as, for instance, optical phonons, dielectric screening, many-body renormalization, and other effects, discussed, for example in Refs. Ponce, and Fischetti2, , are not taken into account in the present study. This can explain some discrepancy with the results obtained from first principles AsSb_mobility ; Wang ; Sohier . Moreover, under realistic conditions, charge carriers are subject to extrinsic scattering by impurities, defects, and substrates, which could further reduce the mobility by orders of magnitude.

V Conclusion

Flexural phonons are the inherent phenomena of 2D materials, which necessarily contributes to the electron-phonon scattering. Unlike conventional in-plane modes, whose contribution to electronic transport is weakly dependent on the carrier concentration, scattering on flexural phonons becomes more effective in the regime of low carrier densities. At the same time, this regime is characterized by the strong anharmonic coupling between in-plane and flexural modes, resulting in the modification of the scattering behavior.

Here we developed a theory for the electron-phonon scattering in 2D semiconductors taking the anharmonic effects into account. We applied this theory to examine the scattering behavior in typical group V hexagonal monolayers. We showed that although the scattering rate is affected by the flexural phonons, the net effect does not exceed 30%. Compared to the harmonic approximation, where the scattering by flexural phonons is divergent in the limit of low concentrations, the anharmonic coupling suppresses this effect significantly. For all materials considered, at carrier concentrations below 1012 cm*-2* flexural phonon contribution turns out to be nearly independent of doping in a wide range of temperatures. Figure 5 summarizes our findings.

From the practical point of view, the suppression of flexural phonons (e.g., by a substrate or by strain) is considered to be one of the ways to enhance the mobility in 2D semiconductors. Our results show that this approach is highly limited and cannot lead to significant mobility gain. On the theory side, flexural phonons can be safely excluded from the consideration in the electronic transport calculations.

Acknowledgements.

A.V.L. and M.I.K. acknowledge support from the research program “Two-dimensional semiconductor crystals” Project No. 14TWOD01, which is financed by the Netherlands Organisation for Scientific Research (NWO). S.Y. acknowledges support from the National Key R&D Program of China (Grant No. 2018FYA0305800). Numerical calculations presented in this paper have been partly performed on a supercomputing system in the Supercomputing Center of Wuhan University. Computational facilities of Radboud University (TCM/IMM) funded by the FLAG-ERA JTC2017 Project GRANSPORT are also gratefully acknowledged. The work is partially supported by the Russian Science Foundation, Grant No. 17-72-20041.

Appendix A Derivation of anharmonic scattering probability

Here we provide details on the derivation of Eq. (12). Treating the phonon fields classically, the thermal average of $\left|V_{{\bf q}}\right|^{2}$ is given by the following functional integral over in-plane and out-of-plane displacement fields:

[TABLE]

where the Boltzmann weight is determined by the Hamiltonian in Eq. (11). The partition function $Z$ is given by:

[TABLE]

Before performing the Gaussian integration over in-plane fields, it is useful to decompose $f_{\alpha\beta}\left({\bf q}\right)$, the Fourier transform of $f_{\alpha\beta}\left({\bf r}\right)=\partial_{\alpha}h\left({\bf r}\right)\partial_{\beta}h\left({\bf r}\right)$, into longitudinal and transverse parts Nelson ; Nelson-orig :

[TABLE]

Here $P^{T}_{\alpha\beta}\left({\bf q}\right)=\delta_{\alpha\beta}-q_{\alpha}q_{\beta}/q^{2}$ is the transverse projector and $f^{T}\left({\bf q}\right)=P^{T}_{\gamma\delta}\left({\bf q}\right)f_{\gamma\delta}\left({\bf q}\right)$. It is then convenient to change variables in the functional integral Nelson ; Nelson-orig and to integrate over shifted in-plane fields, defined in Fourier space as:

[TABLE]

where $\nu=\lambda/(\lambda+2G)$ is the 2D Poisson ratio. In terms of the variables $\tilde{u}_{\alpha}\left({\bf q}\right)$, the Hamiltonian reads

[TABLE]

where $Y=\lambda\left(1-\nu^{2}\right)/\nu$. The Fourier transform of the deformation potential can be written as

[TABLE]

The out-of-plane fields and the in-plane modes $\tilde{u}_{\alpha}\left({\bf q}\right)$ are now decoupled in the Hamiltonian. The thermal average of $\left|V_{{\bf q}}\right|^{2}$ is therefore:

[TABLE]

Since the Hamiltonian is quadratic in the fields $\tilde{u}_{\alpha}$, the correlation function $\langle\tilde{u}_{\alpha}({\bf q})\tilde{u}_{\beta}(-{\bf q})\rangle$ can be calculated explicitly. From the Hamiltonian (26) we find:

[TABLE]

where $P^{L}_{\alpha\beta}\left({\bf q}\right)=q_{\alpha}q_{\beta}/q^{2}$ is the longitudinal projector. Combining Eqs. (28) and (29) yields:

[TABLE]

which is Eq. (12) in the main text.