Coulomb interactions and screening effects in few-layer black phosphorus: a tight-binding consideration beyond the long-wavelength limit

TL;DR

This paper investigates dielectric screening and Coulomb interactions in few-layer black phosphorus using a tight-binding approach that extends beyond the long-wavelength limit, revealing anisotropic screening and plasmon behavior.

Contribution

It introduces a systematic method combining tight-binding and Coulomb interactions to analyze screening effects beyond the long-wavelength approximation in black phosphorus.

Findings

Dielectric function shows strong anisotropy in static limit.

On-site Coulomb interaction is moderately reduced by p_z states.

Full plasmon spectrum exhibits features beyond long-wavelength limit.

Abstract

Coulomb interaction and its screening play an important role in many physical phenomena of materials ranging from optical properties to many-body effects including superconductivity. Here, we report on a systematic study of dielectric screening in few-layer black phosphorus (BP), a two-dimensional material with promising electronic and optical characteristics. We use a combination of a tight-binding model and rigorously determined bare Coulomb interactions, which allows us to consider relevant microscopic effects beyond the long-wavelength limit. We calculate the dielectric function of few-layer BP in the random phase approximation and show that it exhibits strongly anisotropic behavior even in the static limit. We also estimate the strength of effective local and non-local Coulomb interactions and determine their doping dependence. We find that the $p_z$ states responsible for low-energy excitations in BP provide a moderate contribution to the screening, weakening the on-site Coulomb interaction by less that a factor of two. Finally, we calculate the full plasmon spectrum of few-layer BP and discuss the effects beyond long-wavelengths.

Tables

Table 1. Table 1: Computed on-site and nearest-neighbor bare and screened Coulomb interactions in BP within the same layer.

$W_{ij}$(eV)	Bare	Screened			$R_{ij}$(Å)
		1 layer	2 layers	3 layers
$W_{00}$	10.85	6.74	5.96	5.64	0.0
$W_{01}$	5.88	3.27	2.73	2.63	2.22
$W_{02}$	6.20	5.06	4.30	4.07	2.24
$W_{03}$	4.03	2.37	1.88	1.84	3.31
$W_{04}$	3.82	2.88	2.30	2.21	3.47

Equations

$$V_{ij}=e^2\int d\mathbf{r}d{\mathbf{r}}^{\prime }|w_i(\mathbf{r}){|}^2|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^{-1}|w_j({\mathbf{r}}^{\prime }){|}^2,$$

$$V_{ij}(\text{q})=\sum _{\mathbf{R}}{V}_{ij}(\text{R})e^{-i\text{q}\text{R}}+\frac{2\pi e^2}{|\mathbf{q}|}-\sum _{\mathbf{R}}\frac{e^2}{|{\mathbf{r}}_{ij}+\mathbf{R}|}{e}^{-i\text{q}\text{R}},$$

$${ϵ}_{ij}(\text{q},\omega )={\delta }_{ij}-\sum _k{V}_{ik}(\mathbf{q}){\Pi }_{kj}(\text{q},\omega ),$$

$$\begin{array}{rl}{\Pi }_{ij}& (\text{q},\omega )=\frac{g_s}{{\Omega }_0}\sum _{\mathbf{k},mn}\frac{f_m(\text{k})-f_n(\text{k}+\text{q})}{E_m(\text{k})-E_n(\text{k}+\text{q})+\hslash \omega +i\eta }\\ & \times C_{i,m}(\text{k})C_{i,n}^{\ast }(\text{k}+\text{q})C_{j,m}^{\ast }(\text{k})C_{j,n}(\text{k}+\text{q}),\end{array}$$

$${\psi }_{m\mathbf{k}}(\mathbf{r})=\sum _i{C}_{i,m}(\mathbf{k}){\varphi }_{i\mathbf{k}}(\mathbf{r}),$$

$${\varphi }_{i\mathbf{k}}(\mathbf{r})=\sum _{\mathbf{R}}{e}^{i\mathbf{k}\cdot \mathbf{R}}{w}_{i\mathbf{R}}(\mathbf{r}).$$

$$H=\sum _{i\ne j}{t}_{ij}^{||}{c}_i^{†}{c}_j+\sum _{i\ne j}{t}_{ij}^{\perp }{c}_i^{†}{c}_j,$$

$${ϵ}_M(\text{q},\omega )=\frac{1}{N}\sum _{\text{ij}}{ϵ}_{ij}(\text{q},\omega ),$$

$$W_{ij}(\text{q})=\sum _p{ϵ}_{ip}^{-1}(\text{q},0)V_{pj}(\mathbf{q}).$$

$$W_{ij}(\mathbf{q})=N_k^{-2}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }}⟨{\varphi }_{i\mathbf{k}}(\mathbf{r}){\varphi }_{j{\mathbf{k}}^{\prime }}(\mathbf{r})|W|{\varphi }_{i\mathbf{k}+\mathbf{q}}({\mathbf{r}}^{\prime }){\varphi }_{j{\mathbf{k}}^{\prime }-\mathbf{q}}({\mathbf{r}}^{\prime })⟩,$$

$$⟨...|W|...⟩=\int d\mathbf{r}d{\mathbf{r}}^{\prime }{\varphi }_{i\mathbf{k}}^{\ast }(\mathbf{r}){\varphi }_{i\mathbf{k}+\mathbf{q}}(\mathbf{r})W{\varphi }_{j{\mathbf{k}}^{\prime }}^{\ast }({\mathbf{r}}^{\prime }){\varphi }_{j{\mathbf{k}}^{\prime }-\mathbf{q}}({\mathbf{r}}^{\prime })$$

$$\int ϵ(\omega ,\text{r},{\text{r}}^{\prime })\varphi (\omega ,{\text{r}}^{\prime })d{\text{r}}^{\prime }=0,$$

$$\det [{ϵ}_{mn}(\mathbf{q},\omega )]=0,$$

Full text

I Introduction

A few-layer black phosphorus (BP) has recently attracted significant attention as a prospective material for electronic and optical applications [1, 2, 3, 4]. Two-dimensional (2D) BP features a direct layer-dependent energy gap, strongly anisotropic transport characteristic, and high carrier mobility [5, 6, 7, 8, 9, 10, 11, 12, 13]. These properties appear to be particularly promising for plasmonics, where BP is expected to show highly confined, low-loss, and tunable plasmon polaritons [14, 15, 16, 17], as well as it offers the possibility to explore new plasmonic effects, such as hyperbolic optical response [18, 19, 20].

Electron-electron interactions and their screening are known to play a key role in determining the properties of materials. Those include optical properties ranging from the fundamental gap to collective excitations [21], charge carrier transport [22], as well as more exotic examples such as superconductivity [23, 24] and $sp$-magnetism [25, 26, 27]. Some aspects of the Coulomb screening in BP have already been addressed in the literature. Low et al. [14] studied screening in $n$-doped single-layer (1L) BP and despite highly anisotropic band dispersion found essentially isotropic screening in the static limit. In contrast, dynamical screening was found to exhibit strong anisotropy, enabling the existence of plasmons with anisotropic dispersion. Such prediction has been confirmed by electron energy loss spectroscopy (EELS), yet in the context of bulk BP crystal [28]. Similar theoretical results have been reported by Jin et al. [29] with special emphasis on the effects of disorder in 1L- and 2L-BP. Other aspects of electron screening in 2D BP related to strain engineering of plasmons and to the Coulomb drag have been respectively addressed in Refs. 30 and 31.

Up to now, screening in BP has been studied on the basis of the random phase approximation (RPA) for the dielectric function [14, 29, 30, 31]. Moreover, simple form of the Coulomb interaction in reciprocal space ($V\sim 1/q)$ was always assumed, which corresponds to the long-wavelength limit ($q\rightarrow 0$). Within this approximation, dielectric function is being considered as a scalar $\varepsilon(q)=1-V(q)\Pi(q)$, meaning that microscopic effects relevant at short distances are neglected. Although RPA polarizability $\Pi(q)$ is treated as a matrix within the tight-binding (TB) consideration [29, 30], this approximation is apparently not sufficient to capture all microscopic effects in $\varepsilon(q)$. At the same time, screening effects beyond the long-wavelength limit represent a problem of its own significance. For instance, they can be relevant for scattering processes involving large momentum transfer and are relevant in the context of many-body Hamiltonians. In 2D plasmonics, microscopic effects are especially relevant for interlayer interactions, giving rise to the existence of multiple plasmon modes [22, 32]. To date, local Coulomb interaction in BP has only been considered as a parameter [26].

In this contribution, we study screening in a few-layer BP considering effects beyond the long-wavelength limit. We calculate static microscopic dielectric function of black phosphorus for different number of layers ($n=$1–3), as well as determine both short-range and long-range screening of the Coulomb interaction between the relevant BP orbitals. To this end, we use a tight-binding formalism in conjunction with the bare Coulomb interaction determined from realistic charge density distribution in real space. Based on the calculated microscopic dielectric function resolved over the whole Brillouin zone, we restore the plasmon spectra for each system considered and discuss their features arising beyond the limit $q\rightarrow 0$.

The paper is organized as follows. In Sec. II, we consider bare Coulomb interactions in few-layer BP calculated from first-principles. Sec. III is devoted to the static dielectric function calculated at the level of the TB model. Results on the static screened Coulomb interactions are presented in Sec. IV. In Sec. V, we analyze the spectrum of plasmon excitations in BP. In Sec. VI, the paper is concluded.

II Bare Coulomb interaction

The classical form of the 2D Coulomb interaction, $V(q)=2\pi e^{2}/q$, is inapplicable for wave vectors $q\gtrsim 1/a_{0}$, where $a_{0}$ is the unit cell size. To adequately describe the interactions beyond the long-wavelength limit, both radial and angular dependence of the electronic density distribution at short distances should be taken into account. Given that the relevant low-energy states in black phosphorus are predominantly composed of $p_{z}$ orbitals (schematically shown in Fig. 1), the bare Coulomb interaction between single-occupied orbitals residing at sites $i$ and $j$ can be calculated as

[TABLE]

where $w_{i}(\bf r)$ is the Wannier function obtained by projecting the Bloch functions on the orbitals of $p_{z}$ symmetry localized on phosphorus atoms. Here, the Wannier functions are obtained within the formalism of projected Wannier functions [33] on the basis of $GW_{0}$ calculations performed in Ref. 34. To this end, we use vasp package [35, 36, 37] in conjunction with wannier90 code [38]. Integrals appearing in Eq. (1) are evaluated numerically.

Fig. 2 shows the bare Coulomb interaction calculated between the $p_{z}$-like orbitals of phosphorus atoms in bulk BP for different distances. Keeping in mind strong angular dependence of $p_{z}$-like orbitals depicted schematically in Fig. 1, we distinguish between the in-plane and out-of-plane interactions. One can see that at short distances, calculated in-plane interaction (blue line) diverge from $2e^{2}/r$ (black line). To be able to use this data in our calculation of the dielectric function, we transform $V_{ij}$ to the reciprocal space,

[TABLE]

where $i$ and $j$ label now the atoms within the unit cell, and the sum runs over a finite cluster of unit cells separated by distance $|{\bf R}|$. In Eq. (2), the first term describes the short-range part of the interaction, whereas the second and third term ensure numerically accurate behavior of $V_{ij}({\bf q})$ in the long-wavelength limit (${\bf q}\rightarrow 0$).

III Dielectric function

Within the random-phase approximation (RPA), frequency-dependent dielectric matrix $\epsilon_{ij}(\textbf{q},\omega)$ defined in terms of the momentum transfer vector ${\bf q}$ can be written as [39]:

[TABLE]

where $V_{ik}({\bf q})$ is the matrix of Coulomb interactions defined above in Sec. II A, and $\Pi_{kj}(\textbf{q},\omega)$ is the single-partice polarizability matrix. We note that in the present consideration, local field effects are neglected. In the basis of localized orbitals or Wannier functions, the polarizability matrix reads [39]:

[TABLE]

where $C_{i,m}(\textbf{k})$ is the contribution of the $i$-th Wannier function $w_{i{\bf R}}({\bf r})$ to the Hamiltonian eigenstate $\psi_{m{\bf k}}({\bf r})$ with energy $E_{m}({\bf k})$:

[TABLE]

where

[TABLE]

In Eq. (4), $g_{s}=2$ is the spin degeneracy factor, $\Omega_{0}$ is the unit cell volume, $f_{m}({\bf k})=(\mathrm{exp}[(E_{m}({\bf k})-\mu)/k_{B}T]+1)^{-1}$ is Fermi-Dirac occupation factor, $\mu$ is the chemical potential determined by the carrier concentration $n$, and $\eta$ is a broadening term. In our calculations, we used $T=300$ K, $\eta=5$ meV, and $n=10^{13}$ cm*-2* both for electron and hole doping unless stated otherwise. Brillouin zone integration has been performed on a grid of $\sim$106 k-points.

Numerical determination of the polarizablilty matrix $\Pi_{ij}(\textbf{q},\omega)$ requires the knowledge of eigenvalues $E_{m}({\bf k})$ and eigenstates $\psi_{m{\bf k}}({\bf r})$ of a single-particle Hamiltonian. Those must be defined on a dense k-point grid, ensuring accurate evaluation of the Brollouin zone integrals. Such calculations cannot be easily done directly from first principles. The method of model Hamiltonians offers an efficient alternative for studying electronic properties in BP. In our work, we use a tight-binding model proposed for multilayer BP in Ref. 34, which is given by an effective Hamiltonian, describing one electron per lattice site,

[TABLE]

where $i$ and $j$ run over the lattice sites, $t^{||}_{ij}$ ($t^{\perp}_{ij}$) is intralayer (interlayer) hopping integral between the $i$ and $j$ sites, and $c_{i}^{{\dagger}}$ ($c_{j}$) is the creation (annihilation) operator of electrons at site $i$ ($j$). The hopping integrals $t_{ij}$ are parametrized in Ref. 34 on the basis of first-principles $GW_{0}$ calculations, so that the model accurately describes the quasiparticle bands in the vicinity of the gap formed predominantly by the states of $p_{z}$ symmetry. The model is valid in a wide (up to visible light) spectral region and it is applicable to BP with arbitrary number of layers. The model band structure and possible excitations in the vicinity of a band gap are shown in Fig. 3 for 1L-, 2L-, and 3L-BP.

The macroscopic dielectric function $\epsilon_{M}({\bf q},\omega)$ can be calculated from the full dielectric matrix [Eq. (3)] by averaging over all the orbital states [39],

[TABLE]

where $N$ is the number of states per unit cell.

Fig. 4 shows static macroscopic dielectric function $\epsilon_{M}({\bf q})\equiv\epsilon_{M}({\bf q},0)$ calculated for undoped, as well as $n$- and $p$-doped monolayer BP. Overall, one can see that $\epsilon_{M}$ is highly anisotropic, which primarily determined by the anisotropy of the polarizability matrix $\Pi_{ij}({\bf q},0)$. At the edges of the BZ, the screening anisotropy reaches its maximum, yielding $\epsilon_{M}(Y)/\epsilon_{M}(X)\approx 4$ in all cases. The screening anisotropy is therefore more pronounced for wave vectors comparable with inverse unit cell size, $|{\bf q}|\sim 1/a_{0}$. On the contrary, in the vicinity of the zone center ($\Gamma$ point) the dielectric function is essentially isotropic with $\epsilon_{M}({\bf q}\rightarrow 0)\rightarrow\infty$ for doped BP and $\epsilon_{M}({\bf q}\rightarrow 0)=1$ in the absence of doping, which is consistent with generic dielectric properties of 2D materials [40]. Indeed, for the doped case the result for small $q$ corresponds to the Thomas-Fermi approximation yielding $\varepsilon({\bf q})=1+\kappa/|{\bf q}|$, where $\kappa=2\pi e^{2}N(\varepsilon_{F})$ is the screening wave vector and $N(\varepsilon_{F})$ is the density of states. Without doping, the polarization function can be represented at small $q$ in the general form as $\Pi({\bf q})=\sum_{\alpha\beta}C_{\alpha\beta}q_{\alpha}q_{\beta}$, where $C_{\alpha\beta}$ is some finite tensor. Noting that $V\sim 1/q$, one can see from Eq. (3) that $\varepsilon_{M}({\bf q}\rightarrow 0)\rightarrow 1$. We note that the behavior of the macroscopic dielectric function over the whole BZ does not change qualitatively with increasing the number of layers in BP.

IV Screened Coulomb interaction

With the knowledge of the bare Coulomb interaction $V_{ij}(\bf q)$ and the dielectric matrix $\epsilon_{ij}({\bf q},\omega)$, the screened Coulomb interaction matrix $W_{ij}(\textbf{q})$ can be routinely calculated in the static limit ($\omega=0$) as

[TABLE]

$W_{ij}({\bf q})$ should be understood as

[TABLE]

where $N_{k}$ is the number of k-points in the Brillouin zone, and

[TABLE]

is the matrix element of the screened static interaction $W=W({\bf r},{\bf r}^{\prime};0)$, describing the interaction of electrons with orbital indices $i,j$ and momenta ${\bf k}$,${\bf k}^{\prime}$, which involves momentum transfer ${\bf q}$.

Fig. 5 shows calculated screened Coulomb interaction $W_{00}$ in undoped 1L-BP. The screening is more efficient at short distances ($W_{00}/V_{00}\approx$ 0.6), while at distances $r\gtrsim 10$ Å, the interaction is virtually unscreened. Screening anisotropy is less pronounced compared to the reciprocal space (Fig. 4), making the interactions along the zigzag ($x$) direction slightly smaller than those along the armchair ($y$) direction. On-site ($W_{00}$) and nearest-neighbor $W_{01}$ screened interaction in 1L-BP is estimated to be 6.7 and 3.3 eV, respectively, which is somewhat larger than the fully screened interaction predicted for graphene [41]. We note, however, that here we are focused on the screening effects originating exclusively from the $p_{z}$ states of phosphorus, whereas other states of $p$ symmetry as well as high energy states are neglected. In Fig. 6, the dependence of $W_{00}$ and $W_{01}$ on the electron doping $n_{e}$ in 1L-BP is shown. One can see that doping enhances the screening significantly. At experimentally achievable gate doping of $n_{e}=10^{14}$ cm*-2*, $W_{00}$ and $W_{01}$ reach 3.5 and 0.7 eV, respectively. Qualitatively the same results are obtained for the case of hole doping (not shown here).

Table 1 shows calculated screened Coulomb interactions for BP with different (1–3) number of layers. As expected, the interaction strength decreases with the number of layers. This can be attributed to a larger screening associated with a reduced band gap in multilayer BP. Indeed, for smaller band gaps $\Delta$ interband transitions provide larger contributions to the polarization function as $\Pi_{ij}({\bf q},0)\sim 1/\Delta$ [Eq. (4)].

V Plasmons

In the absence of external field, the following equation serve as a criterion for the existance of self-sustained oscillations of the charge density (plasmons),

[TABLE]

where $\epsilon(\omega,\textbf{r},\textbf{r}^{\prime})$ is the dielectric function in a continuous representation and $\phi(\omega,\textbf{r}^{\prime})$ is the potential, induced by charge oscillations. In the discrete representation, Eq. (12) can be rewritten in the form of a secular equation,

[TABLE]

which defines the plasmon spectrum of a system. In the presence of the plasmon damping, it is more convenient to obtain the plasmon spectrum by finding the maxima of the loss funcion $L({\bf q},\omega)=-\mathrm{Im}[1/\prod_{i}\epsilon_{i}(\omega)]$, with $\epsilon_{i}(\omega)$ being the $i$-th eigenvalue of $\epsilon_{mn}({\bf q},\omega)$.

Fig. 7 show the loss function $L({\bf q},\omega)$ of 1-, 2- and 3L-BP calculated for wave vectors up to 0.1 Å*-1* resolved over the $x$ and $y$ directions, and for excitation energies $\omega>\Delta$, where $\Delta$ is the energy gap. For each number of layers considered, one can see the prominent $\omega\sim{\sqrt{q}}$ dependence, which is the classical plasmon dependence originating from the intraband transitions (type 1 in Fig. 3) in 2D electron gas. The corresponding dispersion is anisotropic, being suppressed in the y-direction as a consequence of the anisotropy of the BP electronic structure. Similar results have been previously obtained in the long-wavelength limit for 1L- [14, 30, 29] and 2L-BP [29, 31] using low-energy continuum and TB Hamiltonians.

Deviations from the $\omega\sim{\sqrt{q}}$ dependence occur at $q>0.05$ Å*-1* and become more pronounced as the number of layer increases. Apart from the main plasmon mode, there are additional damped excitations at $\omega<{\sqrt{q}}$, whose intensity $L(q,\omega)\sim q$ in both crystallographic directions, as can be seen from Fig. 8 for all the systems considered. This kind excitations emerge predominantly at $q>0.05$ Å*-1* and thus are typical to short wavelengths. Technically, the origin of those modes is related to the presence of nondiagonal elements in the polarizability matrix $\Pi_{ij}({\bf q},\omega)$ [Eq. (4)], and thus can be interpreted as out-of-phase oscillations of electronic density within or between the layers. Similar behavior has been predicted in metallic bilayers of transition metal dichalcogenides [42]. The description of such features in BP requires to go beyond the continuum low-energy models and is possible at the TB level [29].

In the long-wavelength limit ($q\rightarrow 0$), a charged layer induces a Coulomb potential $v(q)$ decaying with distance $z$ as $e^{-qz}$ [43]. In the context of 2D materials, this results in the emergence of an acoustic plasmon mode with dispersion $\omega\sim q$ in bilayer materials [44]. Such mode has been shown to exist in 2L-BP [29, 31], though it turns out to be strongly damped similar to bilayer graphene [44, 45]. In our calculations, we do not make any assumptions regarding the behavior of the intralayer Coulomb potentials. Having rigorously calculated the Coulomb interaction matrix, $V_{ij}({\bf q})$ [Eq. (2)], we observe strongly damped acoustic plasmon mode at small $q$ for all the cases under consideration including 1L-BP. It is worth noting that the energies reported here for plasmon excitations associated with intraband transitions are somewhat underestimated. This is due to the underestimation of screening effects neglecting the transitions between the states not included into the model Hamiltonian.

For multilayer BP, another type of excitations comes into play, namely, weakly dispersed optical plasmon mode appearing at $\omega\approx 0.6$ eV for 2L-BP, as well as at $\omega_{1}\approx 0.4$ eV and $\omega_{2}\approx 0.8$ eV for 3L-BP (see Figs. 7 and 8). The corresponding frequencies are close to the intraband resonance resulting from the transitions between the subbands (type 2 on Fig. 3), whose splitting is governed by the interlayer interactions. Optical plasmon modes of comparable frequencies have been previously reported for bilayer graphene [45]. It is interesting to note that the damping of optical plasmon excitations in BP is strongly anisotropic. While in $y$-direction the corresponding excitations decay at $q\sim 0.01$ Å*-1*, they are preserved up to $q\sim 0.06$ Å*-1* in the $x$-direction. Although even at small $q$ optical plasmon is damped, its spectral weight is comparable with excitations in the particle-hole continuum ($\omega>\Delta$). Moreover, the spectral weight can be increased by increasing carrier doping or choosing a proper dielectric environment [45].

Optical excitations appearing in Figs. 7 and 8 at $\omega>\Delta$ are not directly related to the screening effects. They are observable without doping and originate from the dipole transitions. Strong anisotropy of this kind of excitations in 1L- and multilayer BP has been extensively analyzed previously using different theoretical methods [46, 47, 34], as well as observed experimentally in bulk BP [28].

VI Conclusion

In this work, we studied dielectric screening and related properties of few-layer BP. We determined dielectric function matrix elements in the orbital subspace using the tight-binding model and a rigorous form of the bare Coulomb interaction computed in real space. Our consideration does not impose any restriction on the wave vector length and, therefore, goes beyond the long-wavelength limit studied before.

Using the random phase approximation, we calculated the static dielectric function of few-layer BP over the whole Brillouin zone, which exhibits strongly anisotropic behavior, especially pronounced at the zone edges. The effective local and non-local Coulomb interactions screened by the $p_{z}$ orbitals included in the TB model are also estimated. Screening is shown to be more efficient at short distances, where the bare interaction is reduced by factor of two, and can be further increased by doping. In real-space, however, the anisotropy of the screened interaction is less evident.

Finally, we calculated the full plasmon spectrum for few-layer BP and classified the origin of different types of excitations. Short-wavelength effects related to the interactions between different sublattices are clearly observable in the plasmon spectrum of all the systems considered including 1L-BP. Apart from the classical $\omega\sim{\sqrt{q}}$ plasmon excitations, we observe additional quasi-linear “acoustic” plasmon mode originating from out-of-phase charge oscillations, as well as a weakly dispersed anisotropic “optical” resonance associated with interband transitions in multilayer BP.

The results presented here provide insights into the dielectric screening in BP at the microscopic level and can serve as a starting point for the analysis of many-body phenomena related to electron-electron coupling, such intrinsic charge carrier transport and superconductivity.

Acknowledgements.

The work was supported by the grant program of the Russian Science Foundation 15-12-20021. Funding from the European Union’s Horizon 2020 Programme under Grant No. 696656 Graphene Core1 is also gratefully acknowledged.