Strain-induced polar discontinuities in two-dimensional materials from combined first-principles and Schr\"odinger-Poisson simulations

TL;DR

This paper demonstrates how strain-induced polar discontinuities in 2D materials like SnSe create free charge wires, using combined first-principles and Schrödinger-Poisson simulations to explore their properties and tunability.

Contribution

It introduces a novel combined simulation approach to study strain-induced charge wires in 2D materials, extending beyond first-principles limitations.

Findings

Strain interfaces in SnSe generate free electron and hole wires.

The Schrödinger-Poisson model accurately reproduces ab-initio results.

Charge wire properties can be tuned by strain parameters.

Abstract

The local application of mechanical stress in piezoelectric materials gives rise to boundaries across which the electric polarization changes. Polarization charges appear along such polar discontinuities and the ensuing electric fields drive a charge reconstruction with the accumulation of free carriers at the boundaries. This is particularly relevant for two-dimensional materials that can sustain very large strains and display record piezoelectric responses. Here we show by first-principles simulations the emergence of one-dimensional wires of free electrons and holes along strain interfaces, taking SnSe as a paradigmatic material. We complement this by developing a Schr\"odinger-Poisson approach specifically designed for two-dimensional materials that it is able to reproduce the ab-initio results and also to extend them to regimes of parameters and system sizes that would be unaffordable in first principles calculations. This model allows us to assess the degree of tunability for the free charge in the wires coming from strain values and profiles, and to obtain the critical size at which the interfaces start to be metallic.

Figures

Tables

Table 1. Supplementary Table 1: Effective masses in units of electron mass obtained by diagonalizing the Hessian of the DFT band structure of bulk SnSe close to a minimum/maximum. Data are reported for ϵ = 0.1 italic-ϵ 0.1 \epsilon=0.1 strain and for the unstrained material. The x ​ x 𝑥 𝑥 xx subscript corresponds to the confinement mass used in the Schrödinger equation and y ​ y 𝑦 𝑦 yy to the mass used in density of states (DOS) calculations.

	Valence Band			Conduction Band
	Unstrained	10%		Unstrained	10%
$m_{xx}^{\Gamma }$	1.756	0.804		2.741	0.985
$m_{yy}^{\Gamma }$	2.733	2.748		2.994	1.593
$m_{xx}^{\Gamma -X}$	0.125	0.206		0.111	0.183
$m_{yy}^{\Gamma -X}$	0.159	1.060		0.190	0.213
$m_{xx}^{\Gamma -Y}$	0.110	0.271		0.132	0.206
$m_{yy}^{\Gamma -Y}$	0.160	0.689		0.130	0.320

Equations

$${\nabla }^2\varphi (x,z)$$

$${\sigma }_{\mathrm{tot}}$$

$$\varphi (x,z)=\int \text{d}{x}^{\prime }{\sigma }_{\mathrm{tot}}(x^{\prime }){\varphi }_{\mathrm{w}}(x^{\prime }-x,z),$$

$${\lambda }_{\mathrm{P}}=-\hat{\mathbf{n}}\cdot ({\mathbf{P}}_0-{\mathbf{P}}_{ϵ}),$$

$${\varphi }_{\mathrm{w}}^{\mathrm{PBC}}(x,z)$$

$$=-\lambda \ln \left\{(x^2+z^2)\prod _{n=1}^{+\infty }\left[(x-nL{)}^2+z^2\right]\left[(x+nL{)}^2+z^2\right]\right\}$$

$$=-\lambda \ln \left\{(x^2+z^2)\prod _{n=1}^{+\infty }(x-nL+iz)(x-nL-iz)(x+nL+iz)(x+nL-iz)\right\}$$

$$=-\lambda \ln \left\{(x^2+z^2)\prod _{n=1}^{+\infty }\left[n^2{L}^2-(x+iz{)}^2\right]\left[n^2{L}^2-(x-iz{)}^2\right]\right\}$$

$$=-\lambda \ln \left\{(x+iz)(x-iz)\prod _{n=1}^{+\infty }{n}^4{L}^4\left[1-\frac{(x+iz{)}^2}{n^2{L}^2}\right]\left[1-\frac{(x-iz{)}^2}{n^2{L}^2}\right]\right\}$$

$$=-\lambda \ln \left\{(x+iz)\prod _{n=1}^{+\infty }\left[1-\frac{(x+iz{)}^2}{n^2{L}^2}\right](x-iz)\prod _{n=1}^{+\infty }\left[1-\frac{(x-iz{)}^2}{n^2{L}^2}\right]\right\}+K$$

$$=-\lambda \ln \left[4\sin \left(\pi \frac{x+iz}{L}\right)\sin \left(\pi \frac{x-iz}{L}\right)\right]+K+c=-\lambda \ln \left[4{\left|\sin \left(\pi \frac{x+iz}{L}\right)\right|}^2\right]+K+c$$

$$=-2\lambda \ln \left|2\sin \left(\pi \frac{x+iz}{L}\right)\right|+K+c$$

$$=-\lambda \ln \left[2\cosh (2\pi z/L)-2\cos (2\pi x/L)\right]+K+c$$

$$\boxed{{\varphi }_{\mathrm{w}}^{\mathrm{PBC}}(x,z)=-2\lambda \ln \left|2\sin \left(\pi \frac{x+iz}{L}\right)\right|=-\lambda \ln \left[2\cosh (2\pi z/L)-2\cos (2\pi x/L)\right]}$$

Full text

I Introduction

A polar discontinuity is a discrete change in electric polarization across an interface or surface. Classical electrostatics tells us that bound polarization charges appear at such polar discontinuities and the system responds with compensating charges in order to screen the corresponding electric fields. In many polar surfaces and interfaces this happens through atomic reconstructions and changes in stoichiometry Harrison1978 ; Noguera2000 ; Goniakowski2008 . Alternatively, a charge reconstruction can occur, with the accumulation of free electrons or holes at the polar discontinuity to screen the polarization charges. This mechanism is considered to be at the origin of the two-dimensional (2D) electron gas appearing at the polar interface between insulating perovskites like LaAlO3 and SrTiO3 Oht2004 ; Bris2014 .

Recently, the emergence of metallic states at polar discontinuities has been extended to 2D materials MartinezGordillo2015 . In this case, one-dimensional (1D) wires of free carriers have been theoretically predicted to appear at polar interfaces between different 2D materials Bristowe2013 ; Gibertini2014 , at boundaries between functionalized and pristine monolayers Gibertini2014 , along the edges of polar nanoribbons  Gul2013 ; Gibertini2014 ; Gibertini2015 , and at inversion domain boundaries Gibertini2015 . Possible applications include electronics, spintronics and solar-energy harvesting. Experimental signatures of such metallic states have been reported in nanostructures based on transition metal dichalcogenides, such as triangular islands of MoS2 Bol2001 and twin boundaries in MoSe2 Liu2014 ; Barja2016 . In the latter case, the rational value of the amount of charge per unit length in the wires (2/3) makes the system susceptible to Peierls instabilities Peierls , giving rise to charge-density waves with the opening of a small gap at the Fermi energy Barja2016 . Although being of utmost fundamental interest, charge density waves might hinder potential applications of these systems.

In this respect, piezoelectricity might offer new strategies to induce polar discontinuities that are more resilient to charge-density waves. In fact, the application of mechanical stress on piezoelectric materials induces a change in their electric polarization. When stress is applied only locally, a boundary appears between strained and unstrained regions, i.e., between regions with different electric polarizations. This means that polar discontinuities occur at strain interfaces in piezoelectric materials. This is particularly relevant for 2D materials that have been shown to withstand very large strains Lee385 ; Akinwande2017 and are endowed with sizable piezoelectric responses Duerloo2012 ; fei2015giant ; Gomes2015 ; Blonsky2015 . In particular, group-IV monochalcogenides like SnSe and SnTe exhibit colossal piezoelectric coefficients fei2015giant ; Gomes2015 and a ferroelectric polarization that can be switched with an external electric field Fei2016 ; Hanakata2016 ; Chang2016 . Here we show by accurate first-principles simulations that interfaces between strained and unstrained SnSe host 1D wires of free carriers as a consequence of the polar discontinuity associated with the piezoelectric response of the material. The amount of free electrons and holes in these wires can be controlled by the value of strain and it is in general not commensurate with the unit length, thus making these metallic states more robust against Peierls instabilities.

In the context of 3D semiconductor heterostructures the self-consistent solution of coupled Schrödinger and Poisson equations has turned out to be extremely useful to predict the charge distribution in doped multilayer systems and even in the case of polar discontinuities in undoped 3D oxides Janotti2012 . In order to gain additional insight into the properties of strain interfaces in SnSe and complement first-principles results, here we also develop a Schrödinger–Poisson (SP) solver specifically designed for 2D materials that properly takes into account the unconventional non-local screening occurring in 2D Keldysh1979 ; cudazzo2011dielectric ; Qiu2016 . The solution of the 2D SP equations has been validated against first-principles results for the charge density and spatial distribution of free carriers at strain interfaces in SnSe. The much cheaper SP approach allows us to extend ab initio calculations and investigate a wide range of strain values and systems sizes that would be otherwise unaffordable. This allows us to assess the tunability of the amount of free charge in the wires by strain and to obtain the critical size at which the interfaces start to be metallic. Finally we show that the assumption of sharp strain interfaces is not necessary and that the effect is robust against a gradual change in strain across an interface.

II Methods

II.1 First-principles simulations

First-principles simulations are performed within density-functional theory (DFT) using the Quantum ESPRESSO Gia2009 suite of codes. The Perdew–Burke–Ernzerhof PBE generalized-gradient approximation of the exchange-correlation functional is adopted. Electron-ion interactions are accounted for using ultrasoft pseudopotentials from the GBRV library GBRV that have been tested SSSP to obtain the best accuracy with respect to all-electron calculations for the elements considered in this work (Sn and Se). An energy cutoff of 40 (320) Ry is employed to expand wave-functions (densities) in all cases, with the exception of variable-cell calculations for which the cutoffs were doubled. The Brillouin zone is sampled using a $14\times 14\times 1$ Monkhorst–Pack grid Monk1976 for bulk 2D materials (both strained and unstrained) while for strain interfaces we used a $1\times 20\times 1$ grid with a 0.005 Ry Marzari-Vanderbilt smearing marzari1999thermal . Bulk structures were carefully relaxed with a threshold of 10*-3* eV/Å for forces and 0.5 kbar for stresses using a Broyden–Fletcher–Goldfarb–Shanno algorithm. 2D systems are simulated within 3D periodic boundary conditions by using the effective-screening-medium approach Otani2006 with 20 Å of vacuum to eliminate the spurious effect of artificial periodic replicas.

II.2 2D Schrödinger–Poisson Solver

In the macroscopic limit we can assume 2D materials to have vanishing thickness and, for definiteness, to lie in the $xy$ plane. We are interested in the case when the system is homogeneous and infinitely extended along the $y$ direction. The electrostatic potential $\phi$ is then a function of $x$ and $z$ only. The Poisson equation links $\phi(x,z)$ to the total charge density $\rho_{\rm tot}(x,z)$, which is the sum of a free $\rho_{\rm f}$ and polarization $\rho_{\rm pol}$ contribution. In this framework, both densities are confined to the plane, i.e. $\rho_{\rm f/pol}(x,z)=\sigma_{\rm f/pol}(x)\delta(z)$, with $\sigma_{\rm f/pol}(x)$ being the planar density in the 2D sheet. The polarization charge density depends on the polarizability $\alpha$ of the material and, generalizing the expression by Cudazzo et al. cudazzo2011dielectric to a spatially-varying $\alpha$, we obtain $\sigma_{\rm pol}(x)=\partial_{x}\left[\alpha(x)\partial_{x}\phi(x,z=0)\right]$. The relevant Poisson equation is hence:

[TABLE]

where

[TABLE]

Given the total charge density, the Poisson equation (1) can be solved using the superposition principle:

[TABLE]

where $\phi_{\rm w}(x,z)$ is the potential produced by an infinite wire along $y$ with unit linear charge density. In open boundary conditions we have $\phi^{\rm OBC}_{\rm w}(x,z)=-\log(x^{2}+z^{2})$, while in periodic boundary conditions $\phi^{\rm PBC}_{\rm w}(x,z)=-\log[2\cosh(2\pi z/L)-2\cos(2\pi x/L)]$, $L$ being the periodicity along $x$ (see Supplemental Material suppl ). Since the right-hand side of Eq. (3) depends on the in-plane electrostatic potential, this equation can be solved self-consistently with Eq. (2); we choose the free charge density as a starting point.

We use this solution of the Poisson equation to implement a multi-band SP solver suited for 2D materials. Since we assume homogeneity along $y$, the Schrödinger equation is one-dimensional and easy to solve numerically using finite-difference techniques brown1990numerical ; tan1990self . We work at zero temperature, in the parabolic-band approximation and use Pulay mixing pulay1980convergence as convergence scheme. In addition to the peculiar nature of screening in 2D described above, another key difference with respect to 3D SP solvers such as in Refs. serra1991one ; karner2007multi ; Busby2010 is that the density of states (DOS) of a one-dimensional electron gas is used to fill the bands. Finally, we introduce some cold smearing marzari1999thermal of the zero temperature Fermi–Dirac distribution for improved convergence and for fair comparison with first-principles calculations. The zero-temperature assumption is required to compare our solver with the DFT results. This assumption can be easily lifted in our solver. The correct materials parameters (like effective masses and position of the band edges), however, should be used, as these change with temperature due to effects like electron-phonon interaction. The code implementing this solver is released open-source on GitHub schrpoisson .

The parameters needed for the SP solver are obtained via DFT calculations on bulk 2D materials. The effective mass tensors are extracted by computing the Hessian matrix of energy bands close to their extrema. In the case studied in the following, effective masses are diagonal in the Cartesian coordinate system, with the $xx$ component representing the confinement mass and the $yy$ component being the DOS mass. Energy bands are aligned by referring them to the vacuum level, defined as the electrostatic potential far away from the 2D layer. Finally, the polarizability $\alpha$ can be computed using density-functional perturbation theory (DFPT) from the scaling of the in-plane dielectric constant in a repeated-layer model with respect to the interlayer separation cudazzo2011dielectric ; Sohier2016 .

Results and discussion

We consider SnSe as a paradigmatic monolayer to illustrate the emergence of polar discontinuities and 1D wires of free electrons/holes across strain interfaces in 2D materials. Owing to its large piezoelectric response fei2015giant ; Gomes2015 , SnSe offers the best platform to engineer large polar discontinuities at moderate strain and it is thus the optimal candidate material to validate our predictions in experiments.

In Fig. 1(a) we show top and lateral views of the crystal structure typical of group-IV monochalcogenides. In the case of SnSe, the rectangular unit cell has relaxed lattice parameters $a=4.41$ Å and $b=4.29$ Å along the $x$ and $y$ directions, respectively. A finite intrinsic electric polarization ${\bm{P}}_{0}$ along $x$ is present as a consequence of the horizontal displacement between Sn and Se atoms. The electromechanical coupling in SnSe allows to tune the polarization value by straining the material (piezoelectric effect) fei2015giant ; Gomes2015 , while its direction can be switched by applying an in-plane electric field (ferroelectric effect) Fei2016 ; Hanakata2016 , which has been recently validated experimentally in a related material, SnTe Chang2016 . The band structure of monolayer SnSe is reported in Fig. 1(b) (black solid lines), showing a finite gap with the valence band maximum along the $\Gamma-$X direction and two almost degenerate conduction band minima along the $\Gamma-$X and $\Gamma-$Y directions, respectively. The zero of energy is set at the vacuum level, that is at the electrostatic potential far away from the monolayer. In the same panel we also plot the band structure when a uniaxial strain $\epsilon=0.1$ is applied along the $x$ direction (red dashed lines), that is parallel to the intrinsic polarization. The energy location of conduction band minima is only marginally affected by strain, while a larger effect is manifest on effective masses (see also Supplemental Material) and on the valence band maximum, giving rise to an increase of band gap. As a consequence of its piezoelectric response, such strained SnSe will also display a different electric polarization. This means that at interfaces between strained and unstrained SnSe, as the ones shown in Fig. 1(c), a polar discontinuity will occur. The linear density of polarization charges accumulating at such strain interfaces increases with the amount of strain (see Supplemental Material) and is given by

[TABLE]

$\hat{\bm{n}}$ being the unit vector normal to the interface and pointing from strained to unstrained regions. The polarization charge density $\lambda_{\rm P}$ crucially depends on the interface orientation and it is maximal in modulus when the interface is along a direction perpendicular to the change in polarization (as in Fig. 1c and in all simulations hereafter).

In the following we consider only uniaxial strain $\epsilon$ along the polarization direction $x$, so that the strained unit cell size along $x$ is $a_{\epsilon}=(1+\epsilon)\,a$, while both strained and unstrained cells have the same length $b$ along $y$. In order to simulate strain interfaces from first-principles within periodic boundary conditions, we take alternating stripes of strained and unstrained SnSe obtained by periodically repeating the structure in Fig. 1(c) along the $x$ direction. We assume to have the same number $N$ of strained and unstrained unit cells, so that the supercell size along $x$ is $L=L_{\epsilon}+L_{0}$, with $L_{\epsilon}=Na_{\epsilon}=(1+\epsilon)Na=(1+\epsilon)L_{0}$. To complement first-principles calculations we take advantage of our SP solver and simulate the system by modelling it as a sequential repetition of strained and unstrained regions of lateral size $L_{\epsilon}$ and $L_{0}$, respectively, separated by sharp interfaces across which effective masses, band edges, and polarizabilities change abruptly. Three valence and three conduction bands, with their respective degeneracies, are taken into account and are associated with the local maxima/minima at the $\Gamma$ point and along the $\Gamma-$X and $\Gamma-$Y directions emphasized in Fig. 1(b). Band extrema are referred to the vacuum level, resulting in larger offsets in the valence bands than in conduction bands (see also Supplemental Figure 3). Polarization charges, computed from DFT using a Berry phase approach King1992 ; Resta2007 , are included as $\delta$-doping at the interfaces in the SP solver Janotti2012 . Unless otherwise stated, we use a discretization step $\Delta x=0.2$ Å to solve the Schrödinger and Poisson equations.

In Fig. 2 we plot the local density-of-states (LDOS) computed with DFT integrated over the vertical $z$ coordinate and averaged over $y$, as a function of energy and position along $x$ for two systems with $L_{0}=26.5~{}$Å and $44.1~{}$Å, corresponding respectively to $N=6$ and $N=10$. Locally a band gap is present between valence and conduction states, but the band edges become position dependent. In order to understand this behaviour we also plot the in-plane macroscopic electrostatic potential (black solid line) that can be obtained from Eq. (3) at $z=0$ by taking as $\sigma_{\rm tot}(x)$ the macroscopic average Bal1988 of the total charge density integrated over $z$ and averaged over $y$. The presence of polarization charges is clearly evident from the non-vanishing electric field appearing throughout the system. Away from the interfaces, the electric field gives rise to a simple shift of the energy states, so that the band edges follow the electrostatic potential. Close to the interfaces, more subtle quantum effects take place as a consequence of confinement associated with the approximately logarithmic potential generated by the polarization charges. To clarify this, we performed SP simulations for the same systems. In Fig. 2 red solid lines represent the SP band edge profiles and show a remarkable agreement with the DFT potential and LDOS away from the interfaces. We also report the energy (red dash-dotted line) and profile (red dashed line) of the lowest-energy state arising from conduction bands in SP simulations. This state is localized at the left interface and quantum confinement endows it with an energy that lies significantly above the SP conduction band profile. This is in very good agreement with the LDOS from DFT and explains the flat band edge profile in the LDOS close to the interfaces and its deviation from the electrostatic potential. We have thus elucidated the importance of quantum effects in the emergence of free carriers at polar discontinuities in 2D materials, contrary to their 3D counterparts where semiclassical arguments can be still applied to explain the relationship between electrostatic potentials and band profiles Bris2014 .

It is also interesting to follow the evolution of this lowest-energy conduction state as a function of the system size. For $L_{0}=26.5~{}$Å ($N=6$) this state is empty and the system is semiconducting, although the gap is considerably reduced as a consequence of the electric fields that tend to shift upwards the valence states on one interface and downwards the conduction states on the other. For $L_{0}=44.1~{}$Å ($N=10$) the system is metallic as this state lies below the Fermi energy and gets partially occupied. This signals the occurrence of an insulator-to-metal transition with increasing $N$. The partial occupation of this conduction state and the corresponding partial depletion of a valence state lead to the formation of pockets of free electrons and free holes. These pockets are localized at opposite interfaces and form 1D metallic wires of free electrons and holes that extend along the interfaces and partially screen the polarization charges. This is even more clear in Fig. 3 where we show the 3D spatial profile of the free carrier density obtained by integrating the LDOS separately for electrons (red) and holes (blue) when $L_{0}=44.1~{}$Å and $\epsilon=0.1$. Averaging these 3D free-carrier densities over $y$ and integrating over $z$ gives the profiles reported in the main panel (thin lines). From SP simulations we obtain density profiles (thick lines) that are in very good agreement with the macroscopic average of the DFT results (dashed lines), validating once more the predictiveness of the 2D SP solver.

The linear density of free electrons/holes present in the 1D wires depends on the system size, as shown in Fig. 4. DFT results (symbols) are reported for three different values of strain ($\epsilon=0.03$, 0.05, and 0.1) and are in remarkable agreement with the same quantity obtained from SP simulations (solid lines), highlighting the predictive quality of the SP solver. An insulator-to-metal transition occurs for any value of strain as the system size is increased. The critical width of the unstrained region at which free carriers start to populate the wires decreases with increasing strain, passing from $\approx 57$ Å for $\epsilon=0.03$ to $\approx 26$ Å for $\epsilon=0.1$. This can be understood from the fact that with increasing strain the change in electric polarization across the interfaces is enhanced, together with the modulus of the polarization charge density in Eq. (4) (see also Supplemental Material). As a consequence, the electric fields are larger and the overall potential drop between the interfaces increases with strain. This favours the transition to a metallic state, that is obtained at smaller critical widths when the localized conduction state on one interface and the hole state on the opposite one meet at the Fermi energy. By increasing the width beyond the critical value, the amount of free carriers accumulated in the wires increases, until asymptotically the free charge density perfectly balances the polarization charge density (dashed lines in Fig. 4) and the overall charge at each interface vanishes.

The excellent agreement between DFT and SP simulations allows us to rely on the much less expensive SP approach to compute the free carrier density over a wide range of system sizes and strain values at the same predictive level of DFT calculations. To reduce even further the computational cost, effective masses, band edges, polarizabilities and all necessary input parameters for SP simulations at arbitrary strain are obtained by fitting the corresponding DFT results over a finite set of strain values with low-order polynomials (see Supplemental Material). Results for the linear charge density of free electrons/holes as a function of strain and lateral size of the unstrained regions are reported in Fig. 5. As expected from the discussion above, the free-carrier density increases with increasing strain and with increasing system size. In particular, this allows us to obtain the critical width of the unstrained region beyond which the system becomes metallic as a function of strain $\epsilon$ in the strained region (black solid line). This critical width decreases with increasing strain, with a minimal value of $\sim 2.6$ nm for the largest value of $\epsilon$ considered here.

Up to now, we have always considered systems with sharp strain interfaces. This choice was mainly dictated by the need to restrain system sizes within the limits of DFT simulations and to simplify the description of the physical interpretation of the results. However, in a realistic setup the strain profile would be a continuous function of position. We want to show now that the phenomena reported above do not require sharp interfaces. Instead, the occurrence of polar discontinuities (and the magnitude of the charge density in the channels) is largely independent of the interface details. To prove this claim, we perform simulations of a setup having the strain profile illustrated in the inset of Fig. 6 using our SP solver, that is capable to consider any smooth strain profile. In particular, we consider a total system width $L$ in the range 70–80 Å, with the strain profile changing from a minimum of zero to a maximum of 8% over a range of length $S$. For illustration purposes, we interpolate linearly the strain between the two constant-strain regions, but we verified that the results would be only marginally affected by a change in the strain profile. The lengths of each of the two constant-strain regions are kept equal to $(L-2S)/2$ in each simulation, and we choose values of S from 0 to 20 Å (that is, up to interface regions very large compared to the system size, with $2S/L$ up to 50%). Although 2D materials can give rise to finite flexoelectric responses Yudin2013 ; Naumov2009 , we do not include them in simulations as they would only affect the profile of polarization charge density and not the overall amount of free carriers in the wires (especially in the asymptotic limit of large system sizes). In Fig. 6 we show a colour plot of the free-carrier density as a function of the geometrical parameters $L$ and $S$. We also show a contour line at constant free-carrier density at the value obtained for $L=70$ Å and $S=0$ Å, i.e. with sharp interfaces. The contour line makes it immediately evident that even for $S>0$ it is possible to obtain the same charge density that would be obtained for sharp interfaces (apart from the need of a slight adaptation of the total system length $L$, which is expected as the effective separation between the wires is reduced for increasing values of $S$). In the Supplemental Material we also show a comparison of the band edge profile and of the first confined states for two different values of $S$, showing that apart from a smoothing of the band profile, the electronic properties of the system are very similar. These results show the feasibility of the approach that we present here to induce polar discontinuities in 2D materials also in realistic situations where the strain might change gradually over several unit cells rather than abruptly.

Conclusions

We have shown that the local application of strain in piezoelectric 2D materials can give rise to polar discontinuities across the interfaces between strained and unstrained regions. Such interfaces could be realized by suspending the sample over an expanding gap or depositing it over neighbouring substrates with different expansion coefficients Duerloo2014 . By using density-functional theory simulations for the prototypical material SnSe, we have demonstrated the occurrence of an insulator-to-metal transition with the appearance of 1D wires of free electrons and holes along the interface. To further characterize the system, we have developed a Schrödinger–Poisson approach specifically designed for 2D materials that takes into account the subtleties associated with screening in low dimensions and that has been validated against first-principles results. The code implementing this solver is released open-source on GitHub schrpoisson and allowed us to compute the free-carrier density over a wide range of system sizes and strain values, and thus to assess the critical width at which the insulator-to-metal transition occurs as a function of strain and the effect of smooth strain profiles.

Acknowledgements

We would like to thank Nicola Marzari for useful discussions and support. This project has received partial funding from the Graphene Flagship through the European Union’s Horizon 2020 research and innovation programme under grant agreement No 696656 GrapheneCore1. We also acknowledge simulation time from the Swiss National Supercomputing Centre (CSCS) through project ID s580.

Supplementary Data

Potential generated by a periodic repetition of charged wires

We are interested in the potential generated by a periodic repetition of charged wires with linear charge density $\lambda$, separated by a distance $L$ along the $x$ direction and extending infinitely along $y$. Although more sophisticated derivations based on complex-plane representations of the problem, we prefer to adopt a more intuitive approach starting from the superposition principle. Indeed, the potential $\phi^{\rm PBC}_{\rm w}(x,z)$ can be simply computed as the sum of the potentials generated by each wire:

[TABLE]

Here is $K$ is an infinite constant that we set to zero by requiring that the potential on a wire generated by all other wires is finite and we made use of Euler’s formula $\sin(\zeta)=\zeta\prod_{n=1}^{\infty}\left(1-\frac{\zeta^{2}}{\pi^{2}n^{2}}\right)$. The second constant $c$ can also be set to zero in order to have that the potential far away from the wires is identical to the the one generated by a charged plane with planar density $\lambda/L$, i.e. $\phi^{\rm PBC}_{\rm w}(x,z\to\infty)=-2\pi\lambda/L|z|$. We thus finally have

[TABLE]