Phonon-limited carrier mobility in monolayer Black Phosphorous

TL;DR

This paper investigates the electron-phonon interactions and resulting carrier mobility in monolayer black phosphorus, revealing anisotropic transport properties and field-dependent drift velocities relevant for electronic applications.

Contribution

It provides a detailed calculation of phonon-limited mobility in monolayer black phosphorus using tight-binding and Boltzmann methods, highlighting anisotropic behavior and high-field effects.

Findings

Mobility ranges from 300-1000 cm^2/Vs in armchair direction at room temperature.

Mobility ranges from 50-120 cm^2/Vs in zig-zag direction at room temperature.

Drift velocity reaches 10^7 cm/s at high fields in the armchair direction.

Abstract

We calculate an electron-phonon scattering and intrinsic transport properties of black phosphorus monolayer using tight-binding and Boltzmann treatments as a function of temperature, carrier density, and electric field. The low-field mobility shows weak dependence on density and, at room temperature, falls in the range of 300 - 1000 cm^2/Vs in the armchair direction and 50 - 120 cm^2/Vs in the zig-zag direction with the anisotropy due to an effective mass difference. At high fields, drift velocity is linear with electric field up to 1 - 2 V/micron reaching values of 10^7 cm/s in the armchair direction, unless self-heating effects are included.

Equations

$$\begin{array}{rl}{E}_g& =2t_2-4|t_1|=1.54\mathrm{eV}\\ m_{ex}& =m_{hx}={\hslash }^2{E}_g/(|t_1|t_2{a}^2)=0.15m_e\\ m_{ey}& =m_{hy}=2{\hslash }^2/(|t_1|b^2)=1.2m_e\end{array}$$

$$\begin{array}{rl}H& =\sum _{<ij>}\frac{k_1}{2}\delta r_{i,j}^2+\sum _{<ijk>}\frac{k_2}{2}\delta {\theta }_{i,j,k}^2\\ & \delta {\theta }_{ijk}=\arccos \frac{\overrightarrow{r_{ij}}\overrightarrow{r_{jk}}}{r_{ij}{r}_{jk}}-\arccos \frac{{\overrightarrow{r_{ij}}}^0{\overrightarrow{r_{jk}}}^0}{r_{ij}^0{r}_{jk}^0}\\ & \delta r_{ij}=|(\overrightarrow{r_j}-{\overrightarrow{r_j}}^0)-(\overrightarrow{r_i}-{\overrightarrow{r_i}}^0)|\end{array}$$

$$H_{el-ph}=\sum_{q\mu kij}M_{q\mu kij}c_{k+qj}^{{\dagger}}c_{ki}\Big{(}a_{q\mu}+a_{-q\mu}^{\dagger}\Big{)}$$

$$\frac{1}{{\tau }_{k\alpha }}$$

$${\mu }_{\alpha }$$

$$e{F}_{\alpha }\frac{\partial f_k}{\partial \hslash k_{\alpha }}=-\sum _{k^{\prime }}{S}_{k{k}^{\prime }}{f}_k(1-f_{k^{\prime }})-S_{k^{\prime }k}{f}_{k^{\prime }}(1-f_k)$$

$${\mu }^{-1}={\mu }_{ac}^{-1}\frac{k_{B}T}{k_B{T}_{rt}}+{\mu }_{op}^{-1}\frac{\exp \left(\hslash {\omega }_{op}/k_B{T}_{rt}\right)-1}{\exp \left(\hslash {\omega }_{op}/k_{B}T\right)-1}$$

$$V=\frac{\mu F}{1+\mu F/V_{sat}}$$

$$\frac{d{E}_g}{d{r}_2}$$

$$\frac{d{E}_g}{d{r}_1}$$

$$eF_{x}\frac{\partial f_{k}}{\partial p_{x}}=I_{c}=-\sum_{\vec{k^{\prime}}}\Big{[}S_{kk^{\prime}}f_{k}(1-f_{k^{\prime}})-S_{k^{\prime}k}f_{k^{\prime}}(1-f_{k})\Big{]}$$

$$f_k=f_{0k}+f_{1k}$$

$$\begin{split}I_{c}=&-\sum_{\vec{k^{\prime}}}f_{1k}\Big{[}S_{kk^{\prime}}(1-f_{0k^{\prime}})+S_{k^{\prime}k}f_{0k^{\prime}}\Big{]}\\ &+\sum_{\vec{k^{\prime}}}f_{1k^{\prime}}\Big{[}S_{kk^{\prime}}f_{0k}+S_{k^{\prime}k}(1-f_{0k})\Big{]}\\ \end{split}$$

$$\begin{split}I_{c}&=-\sum_{\vec{k^{\prime}}}S_{kk^{\prime}}\Big{(}f_{1k}\frac{1-f_{0k^{\prime}}}{1-f_{0k}}-f_{1k^{\prime}}\frac{f_{0k}}{f_{0k^{\prime}}}\Big{)}\\ \end{split}$$

$$f_{1k}=eF_{x}\frac{\partial E_{k}}{\partial p_{x}}\Big{(}-\frac{\partial f_{0k}}{\partial E_{k}}\Big{)}\tau(\vec{k})$$

$$\begin{split}I_{c}&=-\frac{eF_{x}}{kT}\sum_{\vec{k^{\prime}}}S_{kk^{\prime}}f_{0k}(1-f_{0k^{\prime}})\Big{[}v_{kx}\tau(\vec{k})-v_{k^{\prime}x}\tau(\vec{k^{\prime}})\Big{]}\end{split}$$

$$e{F}_x\frac{\partial E}{\partial p_x}\frac{\partial f_{0k}}{\partial E_k}=-\frac{e{F}_x}{kT}{v}_{kx}{f}_{0k}(1-f_{0k})$$

$$\begin{split}1=\sum_{\vec{k^{\prime}}}S_{kk^{\prime}}\frac{1-f_{0k^{\prime}}}{1-f_{0k}}\Big{[}\tau(\vec{k})-\frac{v_{k^{\prime}x}}{v_{kx}}\tau(\vec{k}^{\prime})\Big{]}\end{split}$$

$$\begin{split}&\frac{1}{\tau_{x}(\vec{k})}=\sum_{\vec{k^{\prime}}}S_{kk^{\prime}}\frac{1-f_{0k^{\prime}}}{1-f_{0k}}\Big{(}1-\frac{v_{k^{\prime}x}}{v_{kx}}\Big{)}\\ \mu_{x}&=\frac{j}{enF_{x}}=\frac{\int d\vec{k}v_{kx}f_{1k}}{F_{x}\int d\vec{k}f_{0}}=\frac{e\int d\vec{k}\left(-\frac{\partial f_{0k}}{\partial E_{k}}\right)v_{kx}^{2}\tau_{kx}}{\int d\vec{k}f_{0k}}\end{split}$$

Full text

I Introduction

Recently, black phosphorous emerged as a new two-dimensional (2D) material, which has attracted a great deal of attention due to its excellent transport and optical properties Xia et al. (2014a); Liu et al. (2015a); Ling et al. (2015); Perello et al. (2015); Xia et al. (2014b); Engel et al. (2014); Youngblood et al. (2015); Guo et al. (2016); Buscema et al. (2014); Li et al. (2014a); Chen et al. (2015); Liu et al. (2014); Koenig et al. (2014). Black phosphorous exhibits direct bandgap ranging from 0.3 eV to 1.8 eV for bulk and monolayer films respectively, filling the space between the gapless graphene and wide-gap transition metal dichalcogenides. High mobility at room temperature Perello et al. (2015); Xia et al. (2014b), strong coupling with light, and high anisotropy all contribute to huge potential in application of black phosphorous for infrared imaging Engel et al. (2014), detection Youngblood et al. (2015); Buscema et al. (2014); Guo et al. (2016), and electronic applications Li et al. (2014a); Koenig et al. (2014); Liu et al. (2014); Chen et al. (2015). However, little is known about the interplay of the acoustic and optical phonon scattering as a a function temperature, carrier density, and electric field.

Here we calculate the electron-phonon interactions and drift velocity in a monolayer black phosphorous, within a standard tight binding approach. Our results provide a detailed microscopic picture of phonon scattering. In particular, we find that two optical phonons with energies of around 400 cm*-1* dominate scattering at room temperatures with acoustic phonon contribution being around 30%. The carrier density dependence of the low-field mobility is weak, such that mobilities at rooms temperature vary by a factor of three as carrier density varies by three orders of magnitude from 1011 cm*-2* to 1014 cm*-2*. This is a consequence of a constant density of states in 2D.

II Theoretical Model

For the transport calculations, we use Fermi Golden rule to obtain the electron-phonon scattering rates. The model ingredients include: single-particle band structure fitted to GW calculations from Ref. Rudenko and Katsnelson (2014); valence force phonon model fitted to DFT calculations from Ref. Zhu and Tománek (2014); and electron-phonon coupling constants fitted to our DFT calculations shown in Appendix A.

The single-particle band structure of black phosphorous monolayer can be obtained within the nearest-neighbor interactions in-plane $t_{1}=-1.15$ eV and out-of-plane $t_{2}=3.07$ eV hopping integrals, respectively, as shown in Fig. 1a. The bandgap $E_{g}$ and effective masses in our model are:

[TABLE]

where $m_{e}$ is the electron mass and the lattice parameters $a=4.627$ Å  and $b=3.308$ Å  are along the armchair ($x$) and the zig-zag ($y$) directions, respectively. The bandgap is consistent with the experiment Li et al. (2017) and GW calculations Rudenko and Katsnelson (2014); Tran et al. (2014); Li et al. (2017) and effective masses with the density functional theory calculations Li et al. (2014b); Qiao et al. (2014); Peng et al. (2014). The overall bandstructure of our parametrization of the tight binding model agrees fairly well with that from GW calculations in Ref. Rudenko and Katsnelson (2014) near the top of the valence bands and bottom of the conduction bands, as shown in Fig. 2.

For the phonon model, we use a two parameter Keating-type model Keating (1966):

[TABLE]

where $\delta r_{ij}$ is a distance variation between the neighbouring atoms $i$ and $j$ and $\delta\theta_{ijk}$ is a variation of the angle made up by the three neighbouring atoms $i$, $j$ and $k$. Variables marked with superscript zero correspond to the equilibrium lattice structure Morita (1986); Liu et al. (2015a); Peng et al. (2014). The best fit to the phonon dispersion from the DFT calculations Zhu and Tománek (2014), shown in Fig. 1c, uses the bond stretching stiffness $k_{1}=7.2$ eV/Å2 and the bending stiffness $k_{2}=4.4$ eV.

The electron-phonon interaction is modelled by the Su-Schrieffer-Heeger (SSH) hamiltonian Su et al. (1979) with distance dependent hopping integrals $t_{1}(r_{ij})=t_{1}-g_{1}\delta r_{ij}$ and $t_{2}(r_{ij})=t_{2}-g_{2}\delta r_{ij}$. Using our DFT calculations of the bandgap modulation as a function of an applied in-plane and out-of-plane strain to a black phosphorus monolayer as shown in Appendix A, we find e-ph coupling constants $g_{1}=1.53$ eV/Å  and $g_{2}=2.46$ eV/Å. The Fourier transformed SSH Hamiltonian is

[TABLE]

where $M_{q\mu kij}$ is the e-ph coupling 111Here we neglect terms corresponding to the virtual electron-hole pair excitations across the bandgap., $c_{k+qj}^{{\dagger}}$($c_{ki}$) denotes creation (annihilation) of an electron in the conduction band with index $i,j=1,2$, $a_{-q\mu}^{{\dagger}}$ is a phonon creation operator with wavevector $-q$ and phonon band index $\mu=1\dots 12$. Previous studies of the intrinsic mobilities used the deformation potentials to describe acoustic phonon Qiao et al. (2014); Rudenko et al. (2016) and flexural phonon modes Rudenko et al. (2016) scattering. Inclusion of the optical phonon scattering Liao et al. (2015) predicted lower mobilities.

The low-field electron mobility can be calculated according to Appendix B:

[TABLE]

where $\alpha=x,y$ labels direction of an electric field, $S_{kk^{\prime}}\propto|M|^{2}$ is the electron-phonon scattering rate, $f_{0k}(E_{F}-E_{k},T)$ is the equilibrium Fermi Dirac distribution at temperature $T$ and Fermi energy $E_{F}$, $v_{k\alpha}$ is the band velocity at wavevector index $k$, which absorbs both wavevector direction and electron band index. The hole mobility is identical, due to the electron-hole symmetry in our model.

III Results

In Fig. 1d we show relative contributions of different phonons to the scattering rate in Eq. (4), weighted average by a distribution $\partial f_{0k}/\partial E_{k}$ for $n=10^{13}$ cm*-2* and $T=300$ K. The weighted average of the scattering times in $x$ and $y$ directions are very similar, as shown in Fig. 1d. The areas under both curves in Fig. 1d correspond to the scattering times $\bar{\tau}_{x}=63$ fs and $\bar{\tau}_{y}=62$ fs, which would suggest Drude mobilities to be $\mu_{x}=e\bar{\tau}_{x}/m_{x}=740$ cm2/Vs and $\mu_{y}=e\bar{\tau}_{y}/m_{y}=90$ cm2/Vs, respectively. This is consistent with the full calculations of the low-field mobilities according to Eq. (5) $\mu_{x}=625$ cm2/Vs and $\mu_{y}=82$ cm2/Vs.

For a high field $F$, drift velocity can be obtained from the numerical solution of the Boltzmann transport equation 222The Brillouin zone mesh of 1000 $\times$ 500 and energy cut off 0.5 eV above the bottom of the conduction band were employed to transform Boltzmann transport equation to a matrix form and to solve it iteratively. for the distribution function $f_{k}$:

[TABLE]

We find that mobility from the Boltzmann transport equation solution Eq. (6) differs by 5%-10% from the mobility calculated from Eq. (5).

Fig. 3 shows temperature dependence of the low-field mobility for different concentrations. As expected, phonon-limited mobility decreases with temperature but exhibits several peculiarities. In the armchair direction, the room temperature mobility values are not as large as in graphene, but higher than that in silicon and other currently employed materials confirming high potential of black phosphorus in electronics and optoelectronics. In the zig-zag direction, mobility values are about eight times lower due to the differences in the effective masses.

The overall temperature dependence appears as a power law in Fig. 3a and 3b. However, a more detailed analysis based on the scattering rate from Fig. 1d suggests contributions from both acoustic and optical phonons, which can be accounted for by an empirical expression as motivated by the Matthiessen rule Ashcroft and Mermin (1976):

[TABLE]

where $T_{rt}=300$ K and three fitting parameters are acoustic phonon limited mobility $\mu_{ac}$ at $T_{rt}$, optical phonon limited mobility $\mu_{op}$ at $T_{rt}$, and optical phonon energy $\hbar\omega_{op}$. Using $\mu_{ac}=2570$ cm2/Vs, $\mu_{op}=1140$ cm2/Vs for $x$ direction and $\mu_{ac}=260$ cm2/Vs, $\mu_{op}=160$ cm2/Vs for $y$ direction, and phonon energy of $\hbar\omega_{op}=450$ cm*-1* (or 56 meV), Eq. (7) reproduces very well calculated mobility temperature dependence, as shown in Fig. 3c.

The mobility variation with carrier density at room temperature is weak, as shown in Fig. 4. We attribute this behaviour to the constant density of states in 2D and dominant contribution of the optical phonon scattering at room temperature. At carrier density of about $10^{13}$ cm*-2*, the Fermi level coincides with the optical phonon energy. This causes suppression of the optical phonon emission due to the lack of the final density of states. As a result, mobility increases by roughly a factor of two in the limit of low carrier density. Note, that the carrier density dependence exhibits different behaviour at low temperatures below $70$ K, where acoustic phonon contribution dominates, such that mobility slightly increases as the carrier density increases.

For electronic applications velocity saturation is important. At high electric fields, drift velocity for both directions is shown in Fig. 5a and 5b. We find that for a constant phonon temperature, i.e. isothermal calculations, there is a weak sign of the velocity saturation due to the optical phonon scattering at the experimentally accessible field strengths. This is opposite to the case of graphene, where high energy optical phonon scattering is much stronger then the acoustic phonon scattering which leads to the velocity saturation Meric et al. (2008); Perebeinos and Avouris (2010). Dashed curves in Fig. 5 are fits to the velocity saturation expression:

[TABLE]

where $\mu$ is the low-field mobility and $V_{sat}$ is the saturation velocity. A reference field $F_{ref}$, at which velocity shows saturation, can be defined as $F_{ref}=V_{sat}/\mu$. For the armchair direction, we find $F_{ref}=1.7$ V/$\mu$m and 3.3 V/$\mu$m for $n=10^{12}$ cm*-2* and $n=10^{13}$ cm*-2*, correspondingly, and $F_{ref}=6.3$ V/$\mu$m and 14 V/$\mu$m for the zig-zag direction, respectively. Such high value sof $F_{ref}$ imply an absence of the velocity saturation, except for $n=10^{12}$ cm*-2* case in the armchair direction, when $V_{sat}=1.9\times 10^{7}$ cm/s from the fit to Eq. (8).

Nevertheless, the Joule losses lead to the temperature rise as $T=T_{amb}+jF/r$, where $T_{amb}=300$ K is an ambient temperature, $j=enV$ is the electric current, and $r$ is the thermal conductance. The upper bound for $r$, for a monolayer device on an insulating substrate, corresponds to the zero contact thermal resistance between the black phosphorus and the substrate. Such that $r$ is determined by the thermal conductivity and thickness of an insulating substrate. For example, for an oxide thickness with $h=300$ nm of SiO2 substrate with thermal conductivity of $\kappa=1.4$ W/(mK) Ref. Shackelford and Alexander (2000), a maximum value of $r=\kappa/h\approx 0.47$ kW/(K cm2) is expected. Calculations with the self-heating effects included lead to the velocity saturation at much smaller field of $F_{ref}=0.6$ V/$\mu$m and saturated velocity of $V_{sat}=5.6\times 10^{6}$ cm/s, as shown in Fig. 5a. We also show in Fig. 5a an isothermal calculations at the burning temperature of black phosphorous few tens of nm thick sample Engel et al. (2014). However, in the high bias experiments on few tens of nm samples in Ref. Engel et al. (2014) and Ahmed et al. current was almost linear up to the break down voltage. This may imply an important role of the electron-hole pair excitation across the small bandgap of $E_{g}\approx 0.3$ eV in multilayers at high temperatures and large electric fields. Such that increased $n$ would overwhelm a reduction of $V$ due to the saturation.

IV Conclusions

In conclusion, using a tight-binding model with parameters benchmarked by the DFT calculations, we have predicted phonon limited transport of black phosphorous as a function of temperature, carrier density, and electric field. The low-field mobility is highly anisotropic with values along the armchair direction about eight times higher than that of the zig-zag direction. This difference stems from the high ratio of the effective masses and almost isotropic scattering rate. At room temperatures, optical phonons dominate the scattering leading to decreasing mobilities with increasing carrier densities. While transport data on monolayers are scare for a direct comparison, the available low bias measurements on black phosphorous multilayer samples are consistent with our predictions. We hope our work would stimulate further experiments on monolayer samples, which requires inert gas environments to prevent oxidation Koenig et al. (2014); Favron et al. (2015); Edmonds et al. (2015); Liu et al. (2015b).

Appendix A

DFT derived parameters for the electron-phonon scattering

Using WIEN2k density functional theory (DFT) code Blaha et al. (2001) within the PBE functional Perdew et al. (1996), for a monolayer black phosphorous we obtained the internal parameter $v=0.061$ and the equilibrium lattice constants of $a=4.63$ Å  and $b=3.31$ Å  along the armchair ($x$) and the zig-zag ($y$) directions, respectively. The layer separation in the $z$ direction was set at 1.5 inter layer distance, i.e. $c=17.4$ Å. (For the bulk black phosphorus we obtain values $v=0.110$ and $c=11.12$ Å.) The bandgap $E_{g}$ of 0.9 eV for a monolayer black phosphorus is consistent with the previous theoretical reports Liu et al. (2014); Qiao et al. (2014); Peng et al. (2014); Li et al. (2014b); Tran et al. (2014).

As we apply an out-of-plane strain in (001) direction, the out-of-plane bondlength $r_{2}$ increases, while the in-plane $r_{1}$ bondlength is constant. As a result, the bandgap reduces, as shown in Fig. (6). For the in-plane strain along the (110) direction, the in-plane bondangle is constant and the bandgap increases, as shown in Fig. (7). Using a relationship between the bandgap and the hoping integrals from Eq. (1) in the main text, i.e. $E_{g}=2t_{2}-4|t_{1}|$, we can obtain:

[TABLE]

where $dr_{2}/dr_{1}\approx 0.13$ is due to the geometrical factor, i.e. the out-of-plane P-P bond is not perpendicular to the plane. The best fits of the bandgap change with an applied strain to Eq. (9) and (10) result in the values of $g_{1}=1.53$ eV/Å  and $g_{2}=2.46$ eV/Å.

Appendix B

Relaxation time approximation for an anisotropic material

We are looking for a solution of the Boltzmann transport equation (BTE) for a carrier distribution function $f_{k}$ in an external electric field, without loss of generality, along the $x$ direction:

[TABLE]

where $I_{c}$ is a collision integral, $S_{kk^{\prime}}$ is an electron-phonon scattering rate, and $p_{x}=\hbar k_{x}$ is carrier momentum.

The distribution function can be written in the form

[TABLE]

where $f_{0k}$ is the equilibrium Fermi Dirac distribution. We assume that $f_{1k}$ is small and linear in $F_{x}$. In the collision integral, the zeroth order terms in $F_{x}$ sum to zero and the second order terms in field strength are discarded, such that:

[TABLE]

Using the detailed balance relationship $S_{kk^{\prime}}f_{0k}(1-f_{0k^{\prime}})=S_{k^{\prime}k}f_{0k^{\prime}}(1-f_{0k})$, Eq. (B) transforms to

[TABLE]

Using an ansatz for the distribution function

[TABLE]

the collision integral takes the form:

[TABLE]

where $v_{kx}=\partial E_{k}/\partial p_{x}$ is the band velocity. The left hand side of the BTE, which is linear in $F_{x}$, is given by:

[TABLE]

such that the BTE equation reduces to:

[TABLE]

If the scattering is nearly elastic and scattering time has weak angle dependence, i.e. $\tau(\vec{k}^{\prime})\approx\tau(\vec{k})$, we can arrive to the following expression for the relaxation rate and mobility:

[TABLE]

where $j$ is an electric current density and $n$ is the carrier density. Eq. (18) coincides with Eq. (4)-(5).