Multi-scale approach for strain-engineering of phosphorene

TL;DR

This paper introduces a multi-scale theoretical framework combining mechanical and electronic models to analyze strain effects in phosphorene, enabling large-scale device simulations such as exciton funnels.

Contribution

It presents a novel integrated model that links macroscopic strain to microscopic electronic properties in phosphorene, suitable for large-scale device analysis.

Findings

Successfully modeled inhomogeneous deformation in phosphorene

Predicted electronic property changes under strain

Demonstrated potential for exciton funnel applications

Abstract

A multi-scale approach for the theoretical description of deformed phosphorene is presented. This approach combines a valence-force model to relate macroscopic strain to microscopic displacements of atoms and a tight-binding model with distance-dependent hopping parameters to obtain electronic properties. The resulting self-consistent electromechanical model is suitable for large-scale modeling of phosphorene devices. We demonstrate this for the case of an inhomogeneously deformed phosphorene drum, which may be used as an exciton funnel.

Tables

Table 1. Table 1: Hopping parameters and strain-displacement parameters used in this study

$t_1$ (eV)	$t_2$ (eV)	$t_3$ (eV)	$t_4$ (eV)	$t_6$ (eV)	$t_8$ (eV)	$\beta$	${\kappa }_1$	${\kappa }_2$	${\kappa }_3$	${\kappa }_4$	${\kappa }_5$
-1.25	4.38	-0.106	-0.34	-0.47	0.09	2	0.71	0.27	1.26	-0.39	-0.16

Table 2. Table 2: Decomposition of the nine nearest-neighbor vectors within the unit cell and connecting to the neighboring unit cells according to Eq. ( 20 ).

$(i,j)$	$n_{ij}$	$m_{ij}$	$C_{ij}^{\left(12\right)}$	$C_{ij}^{\left(23\right)}$	$C_{ij}^{\left(34\right)}$
$(1,2)$	$0$	$0$	$1$	$0$	$0$
$(2,3)$	$0$	$0$	$0$	$1$	$0$
$(3,4)$	$0$	$0$	$0$	$0$	$1$
$(1,5)$	$-1$	$0$	$1$	$1$	$1$
$(2,6)$	$0$	$-1$	$-1$	$0$	$0$
$(3,7)$	$0$	$-1$	$0$	$0$	$1$
$(4,8)$	$-1$	$0$	$1$	$1$	$1$
$(4,9)$	$1$	$0$	$-1$	$-1$	$-1$
$(1,10)$	$0$	$1$	$1$	$0$	$0$

Equations

$${\mathcal{H}}_{\mathrm{el}}=\sum _{i\ne j}{t}_{ij}{c}_i^{†}{c}_j,$$

$$E_1=$$

$$E_2=$$

$$E_3=$$

$$E_4=$$

$$E_{\mathrm{g}}\equiv E_3-E_2=4t_1+2t_2+4t_4+2t_6+4t_8,$$

$$E_4+E_1=0,E_4-E_1\equiv \Delta E=4(t_2+t_6)-E_{\mathrm{g}}.$$

$$E_2(k_x,0)\approx$$

$$E_3(k_x,0)\approx$$

$$E_2(0,k_y)\approx$$

$$E_3(0,k_y)\approx$$

$$t_1+t_4=$$

$$t_2+t_6=$$

$$t_3=$$

$$t_8=$$

$$t_{ij}^{\prime }=t_{ij}\exp [-\beta (|{\mathbf{R}}_{ij}^{\prime }|/|{\mathbf{R}}_{ij}|-1)],$$

$$\mathfrak{u}=\left(\begin{array}{ccc}{u}_{xx}& u_{xy}& 0\\ u_{xy}& u_{yy}& 0\\ 0& 0& 0\end{array}\right).$$

$${\mathbf{b}}_{12}^{\prime }=$$

$${\mathbf{b}}_{23}^{\prime }=$$

$${\mathbf{b}}_{34}^{\prime }=$$

$$\mathbf{v}=\left(\begin{array}{c}{\kappa }_1{u}_{xx}+{\kappa }_2{u}_{yy}\\ {\kappa }_3{u}_{xy}\\ {\kappa }_4{u}_{xx}+{\kappa }_5{u}_{yy}\end{array}\right),$$

$${\mathbf{R}}_{ij}=n{\mathbf{a}}_1+m{\mathbf{a}}_2+\sum _k{C}_{ij}^k{\mathbf{b}}_k,$$

$$|{\mathbf{R}}_{ij}^{\prime }|=|{\mathbf{R}}_{ij}|+\frac{{\mathbf{R}}_{ij}\cdot \mathfrak{u}\cdot {\mathbf{R}}_{ij}}{|{\mathbf{R}}_{ij}|}+\frac{\left[C_{(23)}^{ij}\mathbf{v}-(C_{(12)}^{ij}+C_{(34)}^{ij}){\mathbf{v}}_{||}\right]\cdot {\mathbf{R}}_{ij}}{|{\mathbf{R}}_{ij}|}.$$

$$E_{\mathrm{g}}^{\prime }=4t_1^{\prime }+2t_2^{\prime }+4t_4^{\prime }+2t_6^{\prime }+4t_8^{\prime }.$$

$$\Delta E_{\mathrm{g}}\approx \beta u_{xx}(0.11t_1+0.10t_2-4.63t_4-0.65t_6+0.02t_8)+\beta u_{yy}(-1.50t_1+0.13t_2-1.64t_4+0.33t_6-3.54t_8).$$

$$u_{xx}(\mathbf{r})$$

$$u_{yy}(\mathbf{r})$$

$$u_{xy}(\mathbf{r})$$

$${\mathbf{a}}_i^{\prime }={\mathbf{a}}_i+\mathfrak{u}\cdot {\mathbf{a}}_i.$$

$${\mathbf{R}}_i=m_i{\mathbf{a}}_1+n_i{\mathbf{a}}_2+\sum _{k=1}^{N_{\mathrm{bb}}}{C}_i^k{\mathbf{b}}_k,$$

$${\mathbf{R}}_{ij}=m_{ij}{\mathbf{a}}_1+n_{ij}{\mathbf{a}}_2+\sum _{k=1}^{N_{\mathrm{bb}}}{C}_{ij}^k{\mathbf{b}}_k,$$

$${\mathbf{R}}_{ij}^{\prime }={\mathbf{R}}_{ij}+\mathfrak{u}\cdot {\mathbf{R}}_{ij}+\sum _{k=1}^{N_{\mathrm{bb}}}{C}_{ij}^k{\delta \mathbf{b}}_k,$$

$${\mathbf{b}}_k^{\prime }={\mathbf{b}}_k+\mathfrak{u}\cdot {\mathbf{b}}_k+{\delta \mathbf{b}}_k.$$

$$E_{\mathrm{VFM}}=\frac{1}{2}\sum _{ij,i^{\prime }{j}^{\prime }}\left({\mathbf{R}}_{ij}^{\prime }-{\mathbf{R}}_{ij}\right)\cdot {\mathfrak{K}}_{ij,i^{\prime }{j}^{\prime }}\cdot \left({\mathbf{R}}_{i^{\prime }{j}^{\prime }}^{\prime }-{\mathbf{R}}_{i^{\prime }{j}^{\prime }}\right),$$

$$\frac{1}{Y_y}$$

$$({\partial }_x^2\chi )({\partial }_y^{2}w)-2({\partial }_x{\partial }_y\chi )({\partial }_x{\partial }_{y}w)=\left({\kappa }_x{\partial }_x^{4}w+{\kappa }_y{\partial }_y^{4}w+2\sqrt{{\kappa }_x{\kappa }_y}{\partial }_x^2{\partial }_y^{2}w\right)-P_z.$$

Full text

I Introduction

Atomically thin semiconductors, like molybdenum disulfide and phosphorene, have recently attracted a lot of interest due to their potential for optoelectronic applicationsmiau+14 ; chja14 ; liwa+15 . In particular, phosphorene - a single layer of black phosphorous - has been in the focus of interest due to its anisotropic elastic and electronic propertiesxiwa+14 ; feya14a ; qiko+14 ; wepe14 ; waku+15 ; wajo+15 . Moreover, the ability to engineer the physical properties of this material by applying strain offers a range of new possibilitiesRoldan2015 . One example in this regard is the exciton funnel effectfeqi+12 ; sapa+16 , which allows for steering excitons in a desired direction by means of an inhomogeneous band-gap. Exploiting this principle may offer a route towards the development of more efficient solar cells.

The theoretical description of the elastic and electronic behavior of realistic phosphorene-based devices is, in general, computationally far too demanding for ab initio calculations. Valence-force models (VFM) and electronic tight-binding models (TB) provide viable semi-empirical alternatives. A consistent description of the influence of deformations on optoelectronic properties is naturally given by using the VFM to find energetically optimal positions of the atoms in combination with a TB model with distance-dependent hopping parameters. Within this approach electronic transport or optical response can be calculated retaining an accurate description of the system at the atomistic scale.

On the other hand, there are many situations of practical interest where the distribution of (macroscopic) strains is known and one would like to use this information to infer the electronic properties of the systemRoldan2015 . For a given quantity of interest, like the electronic band-gap, one can fit the corresponding strain-induced modification using experimental data or ab initio calculationsvowa+15 . A more comprehensive description of the structural and electronic properties, for example in the case of non-uniform deformations, requires a relation between microscopic displacements of the atoms and the macroscopic strain. Often the Cauchy-Born rule, which states that the atomic positions within the crystal lattice follow the overall strain of the medium, is usedRoldan2015 . However, this approximation applies only to Bravais lattices with a monoatomic basiser08 ; mile+16 and thus fails for phosphorene. Instead, this work is based on the strain-displacement relations obtained by minimizing the VFM of a strained unit-cellmile+16 , which generalize the standard Cauchy-Born approximation.

In this article we develop a theoretical framework that consistently treats both the elastic and electronic degrees of freedom of phosphorene micr16a ; ruka14 ; ruyu+15 and is suitable for large-scale modeling of phosphorene devices. To this end, we put forward a minimal TB model with that accurately describes the electronic band-structure of phosphorene in the vicinity of the $\Gamma$-point. Introducing distance-dependent hopping parameters and using strain-displacement relations obtained from the VFM, we derive analytical expressions that relate the renormalization of the hopping parameters with the modifications of the band gap and effective masses caused by strain. (For moderate values of strain we find that the band-gap correction obtained from the Cauchy-Born approximation differs by a factor $3$ from our results.)

As an application, we study the mechanical and electronic properties of a phosphorene drumhead subjected to uniform pressure for both the small and large deformation regimes. For different central deflections values (the ratio between the maximal deformation height and the drum radius) we determine the non-uniform strain both analytically using the elasticity theory and numerically from the VFM. The agreement is remarkable. We calculate the local band-gap an find that it depends very strongly on the deformation profile, and can be easily enhanced by $10$ up to $15$%.

The paper is structured as follows. In Sec. II we briefly present the lattice structure of phosphorene, introduce a TB model that accurately describes the low energy band structure of phosphorene, and discuss how to account for strain effects on the electronic properties. Sec. III begins with the study of the uniform strain case, which allows it to determine the parameters of the TB model. Next, we analyze the situation of non-uniform strain arising in a pressurized phosphorene drum. Finally, in Sec. IV we present our conclusions and outlook.

II Model

The crystal structure of phosphorene consists of an orthorhombic lattice with lattice vectors $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ and four basis atoms arranged in a puckered structure, as depicted in Fig. 1(a). We denote the atomic positions by $\mathbf{R}_{i}$ with subindex $i$ running from $1$ to $N$, where $N$ is the number of atoms in the phosphorene sheet. The interatomic bond vectors are $\mathbf{b}_{ij}=\mathbf{R}_{j}-\mathbf{R}_{i}$ and the angle between $\mathbf{b}_{ij}$ and $\mathbf{b}_{ik}$ is denoted by $\theta_{jik}$. The equilibrium structure is characterized by the interatomic spacing $d\approx 2.22{\rm\AA}$ and intra- and inter-pucker angles $\theta_{1}\approx 96.5^{\circ}$ and $\theta_{2}\approx 101.9^{\circ}$ kaka+86 . The primitive lattice vectors are thus $\mathbf{a}_{1}=2d\left(\cos(\theta_{1}/2)-\cos(\theta_{2})/\cos(\theta_{1}/2),0,0\right)$ and $\mathbf{a}_{2}=2d\left(0,\sin(\theta_{1}/2),0\right)$, whereas the equilibrium bond vectors are $\mathbf{b}_{12}=d\left(\cos(\theta_{1}/2),-\sin(\theta_{1}/2),0\right)$, $\mathbf{b}_{23}=d\left(-\cos(\theta_{2})/\cos(\theta_{1}/2),0,\sqrt{1-[\cos(\theta_{2})/\cos(\theta_{1}/2)]^{2}}\right)$ and $\mathbf{b}_{34}=d\left(\cos(\theta_{1}/2),\sin(\theta_{1}/2),0\right)$.

II.1 Tight-binding model

To calculate the electronic structure we consider an effective four-band TB Hamiltonian ruka14 ; ruyu+15 , containing four $p_{z}$-like orbitals per unit cell (one per atom), namely

[TABLE]

where indices $i$ and $j$ run over all sites of the phosphorene sheet and $t_{ij}=t_{ji}$ are the hopping parameters. Our aim is to develop a simple TB model that accurately describes the strain-induced reconstruction of the electronic bands close to the $\Gamma$-point.

Let us start by putting forward a minimal model that captures the energy dispersion close to the $\Gamma$-point in the absence of strain, as schematically shown in Fig. LABEL:fig:disp-sketch(b). To describe the electron-hole asymmetry and the anisotropy of the electronic dispersion relation around the $\Gamma$-pointline+14 ; qiko+14 ; feya14 ; pewe+14 , we consider both nearest and next-nearest neighbor couplings within a cutoff radius of $4\,{\rm\AA}$ to obtain the five-parameter TB model put forward in Ref. ruka14, . The hopping parameters are denoted by $t_{1},\cdots,t_{5}$. Fig. 1(a) shows the correspondence between the hopping parameters and interatomic matrix elements connecting orbitals centered at the $i$ and $j$ lattice atoms. We note that $t_{5}$ does not affect the dispersion of the valence and conduction bands and, hence, set $t_{5}=0$ to simplify the model. The agreement between the resulting four-parameter model and DFT band-structure calculations at the $\Gamma$-point with respect to the corresponding effective electron and hole massesline+14 is not very good. We find that the accuracy of the TB model is significantly improved by introducing a single additional hopping parameter $t_{8}$, see Fig. 1(a). Finally, we introduce another hopping parameter ($t_{6}$ in Fig. 1) to allow for reproducing the strain-dependence of the electronic band-gap reported in recent ab initio calculations pewe+14 .

At the $\Gamma$-point our six-parameter TB Hamiltonian has the following eigenenergies

[TABLE]

(We have subtracted a common energy shift $E_{0}=2t_{3}$ from all eigenvalues.) The band gap is given by the difference of the energy of the conduction band ($E_{3}$) and the valence band ($E_{2}$)

[TABLE]

whereas for the other two bands we obtain

[TABLE]

Hence, using the band-gap and $\Delta E$ as an input one can only determine two hopping parameters.

Additionally, close to the $\Gamma$-point we find that

[TABLE]

along the armchair direction and

[TABLE]

in the zigzag direction.

Our results are conveniently cast in terms of the energy-band curvatures (or inverse effective masses), namely, $c_{i,\alpha}=\partial^{2}E_{i}/\partial k_{\alpha}^{2}$ with $i=2,3$ and $\alpha=x,y$. Using $E_{\rm g}$, $\Delta E$, $c_{2,y}$ and $c_{3,y}$ as input parameters, we obtain the following simple analytical relations

[TABLE]

These relations show that one needs more input information to uniquely determine the model parameters. We choose to use the curvature in armchair direction and the strain-dependence of the band-gap energy, a quantity that we address in the following section.

II.2 Strained phosphorene

The application of strain to a phosphorene sheet induces a shift in the interatomic distances, which also changes the material band-structure. In the TB model, the modifications of the electronic properties can be accounted for by hopping parameters $t_{ij}$ which depend on the interatomic distances. Here, we assume that after a deformation the hopping parameters become

[TABLE]

where $\mathbf{R}^{\prime}_{ij}$ ($\mathbf{R}_{ij}$) is the modified (original) vector connecting atoms at sites $i$ and $j$ and $\beta$ quantifies the decay.

Next, we need to establish a relation between the applied strain, which is a macroscopic quantity, and the resulting microscopic shift in atomic positions. As already mentioned, since phosphorene has four atoms in the unit cell, the Cauchy-Born rule does not applyer08 ; mile+16 . Instead, the position of the atoms within the unit cells of a strained sample is obtained by minimizing the microscopic elastic energy with the constraint $\mathbf{a}^{\prime}_{i}=\mathbf{a}_{i}+\mathfrak{u}\cdot\mathbf{a}_{i}$ where $\mathbf{a}_{i}$ is the unstrained primitive lattice vector, $\mathbf{a}^{\prime}_{i}$ is the strained one, and $\mathfrak{u}$ is the strain tensor. The latter has three independent elements and is given by

[TABLE]

These elements lead to the following set of strain-displacement relations for phosphorenemile+16

[TABLE]

where $\mathbf{b}_{ij}$ are the bond vectors connecting the atoms $i$ and $j$ within the unit cell. The vector $\mathbf{v}$ reads

[TABLE]

and $\mathbf{v}_{||}$ is the projection of $\mathbf{v}$ onto the plane of the monolayer. The parameters $\kappa_{i}$ are obtained by minimizing the elastic energy per unit cell with the constraints imposed by Eqs. (10). Note that the component $v_{z}$ accounts for a transversal Poisson effect, i.e., a change of thickness due to strain in the $x-y$ plane. In phosphorene, this effect has recently been addressed by using an ad hoc “3D Cauchy-Born” relation, where the strain-tensor has a non-zero component ${u}_{zz}$ Roldan2015 . In effect, this procedure correctly accounts for the component $v_{z}$, but neglects the components $v_{x}$ and $v_{y}$.

To estimate the modification of the hopping terms due to an applied strain, Eq. (8), it is convenient to write $\mathbf{R}_{ij}$ as a linear combination of the lattice vectors and the bond vectors in the unit cell, namely

[TABLE]

where the index $k$ runs over the pairs of sites $(12),\;(23)$, and $(34)$. $C_{ij}^{k}$ is a matrix of integers, which is explicitly given in Appendix A. After a little algebra, we write $|\mathbf{R}^{\prime}_{ij}|$ as

[TABLE]

We comment at this point on the implications of Eq. (13) on the modification of the dispersion due to an applied shear. Since the component $v_{y}$ is proportional to the shear, the last term in Eq. (13) implies that shear, in contrast to what one expects from the Cauchy-Born rule, does not preserve the lengths of bond vectors with a non-vanishing $y$-component. On the other hand, the inversion symmetry of phosphorene implies that shear cannot contribute to the electronic structure to linear order. Consistent with this, we find that the contributions to the electronic structure from an applied shear cancel to linear order, leaving the band structure unchanged for small shearsali+15 .

The band structure of strained phosphorene is obtained by using the modified hopping terms $t_{ij}^{\prime}$ from Eq. (8) in the TB Hamiltonian. In particular, the band-gap of strained phosphorene is given by

[TABLE]

To lowest order in strain, the modified effective masses are obtained by inserting $t^{\prime}_{ij}$ into Eqs. (5) and (6), which describe the electronic dispersion close to the $\Gamma$-point.

III Results

III.1 Parameter estimation

We have described the general procedure for obtaining a TB model to calculate strain-induced changes on the electronic structure of phosphorene. Let us now determine the model parameters introduced in the previous section.

We obtain the strain-displacement parameters $\kappa_{i}$ by minimizing the VFM put forward in Ref. micr16a, . The results are given in Table 1. Next, with the help of Eq. (14) we write the strain-induced modification of the band-gap, $\Delta E_{\rm g}=E_{\rm g}^{\prime}-E_{\rm g}$, as

[TABLE]

The hopping parameters are estimated by fitting the main features of the low energy band structure ruyu+15 , namely, $E_{\rm g}=1.84\;{\rm eV}$ and the effective masses $m_{\rm c}(X)\approx 0.2m_{\rm e}$, $m_{\rm v}(X)\approx 0.2m_{\rm e}$, $m_{\rm c}(Y)\approx 1.2m_{\rm e}$, and $m_{\rm v}(Y)\approx 3.9m_{\rm e}$. Further, we set $E_{4}=-E_{1}=6.9\;{\rm eV}$. As explained in Sec. II.1, these inputs are not sufficient to determine all six hopping parameters and $\beta$. To this end, we compare the strain-induced modification of the band-gap predicted by Eq. (15) with ab initio results pewe+14 . The latter show that the band-gap increases linearly by approximately $0.1\,{\rm eV}$ with strains up to $\approx 4\%$. The increase is larger for strain in armchair direction. By using $\beta=2$ and the hopping parameters given in Table 1 our TB model nicely reproduces this behavior, as shown in Fig. 2(b). The model is also in good quantitative agreement with the electronic dispersion of Ref. ruyu+15, at zero strain, see Fig. 2(a).

We verify the accuracy of our analytical results by comparing the band gaps predicted by Eq. (15) with those obtained from numerical calculations. For that, we consider systems periodic boundary conditions with supercells with $15\times 15$ unit cells. For a given strain, the system is allowed to relax using the VFM developed in Ref. micr16a, . The band gaps obtained by diagonalization of the TB Hamiltonian for the resulting atomic configurations agree very well with Eq. (15), as shown by Fig. 2(b). The 2D and 3D Cauchy-Born relations (dotted and dashed lines) lead to significant deviations from the numerical results, with the latter giving rise to the largest errors.

Interestingly, we find that the magnitude of $t_{4}$ (or $t_{6}$ depending on which parameter is undetermined) is crucial to correctly describe the behavior of the band gap versus strain. Maintaining the constraints for the band-gap and the effective masses, we use Eq. (15) to calculate the ratio of $\Delta E_{\rm g}$ for stretching along the armchair and zigzag directions. Figure 2(c) shows that if $t_{4}\gtrsim-0.24\,{\rm eV}$, the ratio becomes smaller than one, which disagrees with the findings of Ref. pewe+14, . In the present case, $t_{4}=-0.34\,{\rm eV}$ leads to a ratio of $\approx 1.64$. Further investigations and, in particular, experimental results are needed to better assess $\Delta E_{\rm g}$.

III.2 Nonuniform strain

Now we turn to an example of practical interest: a suspended phosphorene drum and the influence of nonuniform strain. If the drumhead is subjected to uniform pressure, it will be statically deformed, which leads to a nonuniform strain distribution. This system is conveniently described by continuum elasticitymicr16a . The phosphorene monolayer is effectively modeled as a thin anisotropic plate. Its macroscopic elastic properties are characterized by two bending rigidities and four stiffness constants. We find the shape of the deformed drum and the strain distribution by solving the equations corresponding to the out-of-plane displacement-field $w(\mathbf{r})$ and for the Airy stress function $\chi(\mathbf{r})$ micr16a . In following we focus our discussion on approximate analytic solutions (details are given in Appendix B) and compare them with numerical results.

For sufficiently small pressure the bending contribution to the elastic energy of the drum dominates and its shape is well approximated by $w(x,y)=w_{0}[1-(x/R)^{2}-(y/R)^{2}]^{2}$, where $w_{0}$ is the deflection at the center and $R$ is the radius of the drum. In contrast, for large pressures stretching plays a major role and the shape is approximately given by $w(x,y)=w_{0}[1-(x/R)^{2}-(y/R)^{2}]$. As shown in Appendix B, the strain distribution in both regimes becomes

[TABLE]

where $P_{1}$, $P_{2}$ and $P_{3}$ are polynomials containing even powers of $x/R$ and $y/R$ (see Appendix B). The tensile strains are maximal and the shear vanishes at the center of the drum. From Eq. (15) we find that the corresponding change in band-gap is maximal and proportional to $(w_{0}/R)^{2}$, with a factor of proportionality of approximately $6.4\,{\rm eV}$ for both deformation regimes.

In order to verify the analytical results, we use the VFM of Ref. micr16a, to numerically calculate the shape of a pressurized drum. To validate our approach in the continuum limit while maintaining computational convenience we simulate a drum of radius $103.16$ Å. The out-of-plane deformation field $w(\mathbf{r})$ is obtained by considering the midpoints of each primitive unit cell. In Figs. 3(a) and 3(b) the deformation field is plotted along $x=0$ and $y=0$, respectively, for three central deflections $w_{0}/R$, together with the analytical expressions for the deformation in the bending and stretching regimes. We find that the drum shape is indeed very close to being radially symmetric and shows the expected crossover from the bending to the stretching dominated regime as $w_{0}/R$ increases.

Further, we consider the vectors connecting the midpoints of adjacent primitive unit cells to estimate the local strain distribution. We test the accuracy of our strain-displacement relations, by calculating the expected lengths of the bond vectors using Eq. (10), and comparing them to the numerically obtained ones. The maximal error of the strain-displacement relations occurs at the center of the drum, where the strain is maximal. Fig. 3(c) shows the relative errors as a function of central deflection together with the relative error given by the 3D Cauchy-Born rule, ignoring $\mathbf{v}_{||}$. At the largest deflection, the accuracy of the Cauchy-Born approximation is about $3\%$ at the largest deflection, whereas our method leads to bond length relative errors of less than $0.5\%$. These differences have a significant impact on the hopping matrix elements $t_{ij}$ and, hence, on the tight-binding electronic structure calculations.

For each unit cell we also calculate the (local) electronic band-structure using the TB model and a $3\times 3$ super-cell centered at the unit cell of interest and using periodic boundary conditions. This allows us to obtain the local band-gap $E_{\rm g}(\mathbf{r})$. In Fig. 4(a) we show maps of the band-gap for three different central deflections, indicating the transition from a bending-dominated regime at low deflections to a stretching dominated regime at high deflections. Figure 4(b) shows a comparison between the numerically calculated band-gap at the center of the drum with our analytical result. The agreement is excellent, suggesting that the model can be used to give qualitative predictions of local electronic properties in complex geometries.

IV Conclusions

In summary, we have presented a multi-scale approach to calculate electronic properties of deformed phosphorene. To this end, we developed a TB model to describe the electronic structure. We found that six hopping parameters are required to get a good quantitative description of the electronic bands, which includes the band-gap and the (anisotropic) effective masses at the $\Gamma$-point. The influence of deformations is described by considering distance-dependent hopping-parameters. The crucial point of our approach is the fact that in almost all cases of interest, the macroscopic strain-distribution can be inferred, but the microscopic positions of the atoms after deformation are not known. Thus a relation between strain and atomic displacements is required to make use of the TB model. As we have shown, the simple Cauchy-Born relation - often employed in this context - is not a good approximation as it tends to overestimate the strain. Instead, we propose to use strain-displacement relations, which are obtained by minimizing the elastic energymile+16 . This approach is easily combined with the microscopic TB model and yields good quantitative agreement. The central result is given by Eq. (15), which gives an expression for the strain-induced modification of the band-gap for a given (homogeneous) strain.

We demonstrate our method for the relevant case of a phosphorene drum. The deformation due to pressure leads to inhomogeneous strain distributions. Introducing a local band-gap, one then finds a spatially changing strain-induced contribution, which is essential for the so-called inverse funneling effectsapa+16 . We obtain a very good agreement between fully numerical results (TB and VFM) with analytical estimates resulting from continuum elasticitymicr16a combined with the derived strain-induced modification of the band-gap given by Eq. (15).

The multi-scale approach presented in this article can also be used for other 2D materials, provided suitable TB and VF models are known. It offers a quantitative and efficient procedure to study the structural and electronic properties of single-layered materials with arbitrary deformation landscapes.

Acknowledgements.

C.H.L. acknowledges financial support of the Brazilian funding agencies CNPq and FAPERJ.

Appendix A Strain-displacement relation for phosphorene

A detailed description of the procedure to obtain strain-displacement relations, Eqs. (10), from a VFM is given in Ref. mile+16, . Here we revisit the theory to explain and provide further insight on some of the expressions used in Sec. II.2.

For a homogeneous strain the primitive lattice-vectors change according to

[TABLE]

The position $\mathbf{R}_{i}$ of any arbitrary lattice site can be expressed in terms of the primitive lattice-vectors and the $N_{\rm bb}$ bond-vectors of the basis,

[TABLE]

where $m_{i}$, $n_{i}$, and the entries of the matrix $C_{i}^{k}$ are integers. The vector $\mathbf{R}_{ij}$ connecting the atoms $i$ and $j$ is then given by

[TABLE]

where for all variables $X_{ij}=X_{i}-X_{j}$ holds. Since we assume that the deformation preserves the lattice structure, $m_{ij}$, $n_{ij}$ and $C_{ij}^{k}$ remain unchanged, upon strain

[TABLE]

which implies that the basis bond-vectors transform as

[TABLE]

For a monoatomic basis, where $N_{\rm bb}=0$, Eq. (20) is reduced to the standard Cauchy-Born relation. The lattice symmetries relate the vectors $\mathbf{\delta b}_{k}$ to each other and thus reduce the number of unknown componentsmile+16 . In the case of phosphorene, one has $\mathbf{\delta b}_{23}=\mathbf{v}$ and $\mathbf{\delta b}_{12}=\mathbf{\delta b}_{34}=\mathbf{v}_{\parallel}$, which are explicitly given by Eq. (11). The nearest neighbors of the atoms within a unit cell are shown in Fig. 5. The corresponding values of $m_{ij}$, $n_{ij}$ and $C_{ij}^{k}$ are given in Table 2.

In general, the vectors $\mathbf{\delta b}_{k}$ can be determined by minimizing the elastic energy of the deformed structure, that for small deformations readsasme76

[TABLE]

where the sum runs over all pairs of atoms and $\mathfrak{K}_{ij,i^{\prime}j^{\prime}}$ denotes the matrix of force-field parameters. The latter is typically restricted to neighboring bonds. Using Eq. (20) one finds that $E_{\rm VFM}$ is a second order polynomial in $\mathfrak{u}$ and $\mathbf{\delta b}_{k}$. By requiring that $\partial E_{\rm VFM}/\partial\mathbf{\delta b}_{k}=0$ one obtains a set of linear relations between $\mathbf{\delta b}_{k}$ and $\mathfrak{u}$, which constitutes the strain-displacement relations.

Appendix B Strain distribution for a drum

Denoting the displacement field by $\mathbf{u}=(u(\mathbf{r}),v(\mathbf{r}),w(\mathbf{r}))$ and the Airy stress function by $\chi(\mathbf{r})$, the shape of the deformed drum is determined by the following equationsmicr16a

[TABLE]

and

[TABLE]

Here, $Y_{x,y}$ and $\kappa_{x,y}$ are the Young’s modulus and the bending rigidity in $x$ and $y$ direction, $G_{xy}$ is the shear modulus and $\nu_{xy}$ and $\nu_{yx}$ are the Poisson ratios.

In the following, we solve these equations for a drum of radius $R$ with vanishing prestrain under a spatially constant external pressure for the limit cases of small and large deformations. Taking into account the boundary conditions for the in-plane displacement fields, namely, $u(\tilde{x}=0,\tilde{y}=0)=0$, $v(\tilde{x}=0,\tilde{y}=0)=0$, $u(\tilde{x}^{2}+\tilde{y}^{2}=1)=0$, $v(\tilde{x}^{2}+\tilde{y}^{2}=1)=0$, we find the Airy stress function and, hence, the strain fields.

For sufficiently small deformations one can ignore the left hand side in (24) which is cubic in the deformation. The deformation is then given by $w(x,y)=w_{0}(1-\tilde{x}^{2}-\tilde{y}^{2})^{2}$ where $\tilde{x}=x/R$ and $\tilde{y}=y/R$. The maximal deflection $w_{0}$ is related to the applied pressure by

[TABLE]

where $\kappa_{\rm eff}\equiv(3\kappa_{x}+3\kappa_{y}+2\sqrt{\kappa_{x}\kappa_{y}})/8\approx 3.9$ eV for phosphorenemicr16a .

We consider a generic 8th order polynomial as Ansatz for the Airy stress function and solve Eq. (23) together with the boundary conditions for the drum in-plane displacement fields by coefficient matching. The resulting strain distributions are sixth order polynomials in $\tilde{x}$ and $\tilde{y}$ with coefficients that depend on intricate combinations of the elastic constants. Taking the explicit values for the elastic constants given in Ref. micr16a, , we find

[TABLE]

To gain insight into the interpretation of these expressions, we solve the Airy stress function for an isotropic material with hexagonal symmetry and small Poisson ratio ($\nu\ll 1$), such as graphene. In that case, we find a simpler expression for the strain fields that is independent of the elastic constants, namely

[TABLE]

Comparing the leading order coefficients of the polynomials in Eqs. (B) and (B) we find that despite the anisotropy of phosphorene, the strain distribution close to the center of the drum, $\tilde{x}\ll 1$ and $\tilde{y}\ll 1$, deviates only by about $30\%$ from what one would expect for an isotropic material.

In the opposite limit of large deformations, we ignore the bending related terms in (24) and consider the Ansatz $w(x,y)=w_{0}(1-\tilde{x}^{2}-\tilde{y}^{2})$, $\tilde{x}=x/R$ and $\tilde{y}=y/R$ for the drum deformation. Inserting this Ansatz into (23) we find that the right-hand-side of the equation becomes spatially constant, namely, $(\partial_{x}\partial_{y}w)^{2}-(\partial^{2}_{x}w)(\partial^{2}_{y}w)=-(4w_{0}^{2}/R^{4})$. By inspection of (23), this implies that the Airy function is given by a fourth order polynomial whose coefficients are obtained by matching the boundary conditions at the rim of the drum. We find the strain distribution

[TABLE]

where $\gamma_{x}$ and $\gamma_{y}$ are functions of the elastic parameters. For phosphorene one has $\nu_{xy}\ll 1$, for which one finds $\gamma_{y}\approx 0$ and $\gamma_{x}\approx(2/3)\nu_{yx}\left(C_{66}-2C_{12}\right)/\left(6C_{12}+C_{66}\nu_{yx}\right)$. Using the values for the elastic constants provided in Ref. micr16a, we find that $C_{66}\approx 2C_{12}$, and therefore $\gamma_{x}\approx 0$. Interestingly, despite phosphorene being anisotropic, this result coincides with the case of isotropic materials with small Poisson ratio ($\nu\ll 1$), where $\gamma_{x}=\gamma_{y}=0$.