Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene

TL;DR

This paper develops a density functional perturbation theory tailored for two-dimensional heterostructures in field-effect setups, enabling detailed first-principles analysis of vibrational and electronic properties relevant for device design.

Contribution

The authors introduce a novel implementation of density functional perturbation theory that includes a truncated Coulomb interaction and simulates charging via field-effects for 2D heterostructures.

Findings

The method accurately computes vibrational properties and electron-phonon interactions in 2D heterostructures.

In doped graphene, flexural phonon coupling is strongly screened and negligible.

Contradicts previous reports by showing gate-induced coupling is not detectable in doped graphene.

Abstract

The ability to perform first-principles calculations of electronic and vibrational properties of two-dimensional heterostructures in a field-effect setup is crucial for the understanding and design of next-generation devices. We present here an implementation of density functional perturbation theories tailored for the case of two-dimensional heterostructures in field-effect configuration. Key ingredients are the inclusion of a truncated Coulomb interaction in the direction perpendicular to the slab and the possibility of simulating charging of the slab via field-effects. With this implementation we can access total energies, force and stress tensors, the vibrational properties and the electron-phonon interaction. We demonstrate the relevance of the method by studying flexural acoustic phonons and their coupling to electrons in graphene doped by field-effect. In particular, we show that while the electron-phonon coupling to those phonons can be significant in neutral graphene, it is strongly screened and negligible in doped graphene, in disagreement with other recent first-principles reports. Consequently, the gate-induced coupling with flexural acoustic modes would not be detectable in transport measurements on doped graphene.

Figures

Equations

$$n({\mathbf{r}}_p,z)$$

$${\psi }_{\mathbf{k},s}({\mathbf{r}}_p,z)$$

$$V_{\mathrm{KS}}^{2D}({\mathbf{r}}_p,z)=V_{\mathrm{ext}}^{2D}({\mathbf{r}}_p,z)+V_{\mathrm{H}}^{2D}({\mathbf{r}}_p,z)+V_{\mathrm{XC}}^{2D}({\mathbf{r}}_p,z).$$

$$f({\mathbf{r}}_p+{\mathbf{R}}_p,z)=f({\mathbf{r}}_p,z).$$

$$f({\mathbf{G}}_p,z)$$

$${\int }_{-t/2}^{t/2}⟨n{⟩}_p(z)dz$$

$$\int ⟨n{⟩}_p(z)dz=n_0=\sum _a\frac{Z_a}{S}+n_{\mathrm{dop}},$$

$$\int ⟨{\varrho }_{\mathrm{ion}}+{\varrho }_{\mathrm{elec}}{⟩}_p(z)dz$$

$${\varrho }_{\mathrm{gate}}(z)$$

$$V_{\mathrm{gate}}^{2D}(z)$$

$$\int ⟨{\varrho }_{\mathrm{ion}}+{\varrho }_{\mathrm{elec}}{⟩}_p(z)+{\varrho }_{\mathrm{gate}}(z)dz=0.$$

$$V_{\mathrm{ext}}^{2D}({\mathbf{r}}_p,z)=V_{\mathrm{gate}}^{2D}(z)+V_{\mathrm{ion}}^{2D}({\mathbf{r}}_p,z),$$

$$⟨V_{\mathrm{KS}}^{2D}{⟩}_p(+\infty )-⟨V_{\mathrm{KS}}^{2D}{⟩}_p(-\infty )=4\pi e^2{n}_{\mathrm{dop}}|z_{\mathrm{g}}|,$$

$$|{\mathbf{F}}_{\mathrm{gate}-\mathrm{material}}|=S\times 2\pi e^2{n}_{\mathrm{dop}}^2.$$

$$V_{\mathrm{barrier}}^{2D}(z)$$

$$V_{\mathrm{KS}}({\mathbf{r}}_p,z)$$

$$V_{\mathrm{KS}}({\mathbf{r}}_p,z+ic)$$

$${\bar{v}}_c(\mathbf{r})=\frac{e^2\theta (l_z-|z|)}{|\mathbf{r}|},$$

$$\bar{V}(\mathbf{r})=\int \frac{e\varrho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}\theta (l_z-|z|)d{\mathbf{r}}^{\prime }.$$

$$l_z=\frac{c}{2}.$$

$${\bar{V}}_{\mathrm{gate}}({\mathbf{r}}_p,z)$$

$$=V_{\mathrm{gate}}^{2D}({\mathbf{r}}_p,z)\text{ if }z\in [z_{\mathrm{g}}-\frac{c}{2};z_{\mathrm{g}}+\frac{c}{2}].$$

$$\bar{V}({\mathbf{G}}_p,G_z)$$

$$\begin{array}{rl}{\bar{v}}_c({\mathbf{G}}_p,G_z\ne 0)& =\frac{4\pi e^2}{|{\mathbf{G}}_p{|}^2+G_z^2}\times \left[1-e^{-|{\mathbf{G}}_p|l_z}\cos (G_z{l}_z)\right],\\ {\bar{v}}_c(\mathbf{G}=0)& =0.\end{array}$$

$${\bar{V}}_{\mathrm{KS}}({\mathbf{r}}_p,z)={\bar{V}}_{\mathrm{ext}}({\mathbf{r}}_p,z)+{\bar{V}}_{\mathrm{H}}({\mathbf{r}}_p,z)+V_{\mathrm{XC}}({\mathbf{r}}_p,z).$$

$${\bar{V}}_{\mathrm{ext}}(\mathbf{r})$$

$$v_a(|\mathbf{r}|)$$

$$v_a^{\mathrm{SR}}(|\mathbf{r}|)$$

$$v_a^{\mathrm{LR}}(|\mathbf{r}|)$$

$$v_a^{\mathrm{SR}}(\mathbf{G})$$

$$v_a^{\mathrm{LR}}(\mathbf{G})$$

$${\bar{v}}_a^{\mathrm{LR}}(\mathbf{G})$$

$${\bar{V}}_{\mathrm{ion}}^{\mathrm{loc}}(\mathbf{G})$$

Full text

I Introduction

Density Functional TheoryHohenberg and Kohn (1964); Kohn and Sham (1965); Ihm et al. (1979); Marting (2004); Engel and Dreizler (2011); Parr and Yang (1989) (DFT) based on plane-wave basis sets, Kohn-Sham equations and pseudopotentials, has proven to be a valuable tool to understand and predict electronic and structural properties of materials. DFT provides support in the process of understanding and controlling experimentally observed phenomenaFerrari et al. (2006). Inversely, it can be usedMounet et al. (2015) to identify interesting new two-dimensional (2D) compounds, thus encouraging their experimental study. The wide range of potential applications and fascinating phenomena offered by 2D materials would benefit from accurate DFT simulation in the 2D framework.

In general, the response of a material to a long wavelength periodic perturbation is highly dependent on dimensionality. We have recently shown the importance of working in the appropriate 2D framework for the computation of dielectric responsesSohier et al. (2015), optical phonons dispersionsSohier et al. (2017), and electron-phonon coupling in polar materialsSohier et al. (2016). Furthermore, a particularly relevant aspect of 2D materials is their sensitivity to external perturbations like external electric fields, as illustrated by the large gate-induced doping achievable for a 2D material in a field-effect transistor (FET) setupNovoselov et al. (2004); Bolotin et al. (2008); Efetov and Kim (2010). While the field effect is not part of the 2D material per se, its omnipresence in experimental setups and devices motivates the need to simulate it.

However, current implementations of DFT with three-dimensional periodic boundary conditions (3D PBC) are not adequate to the simulation of 2D materials doped in the FET setup. This is mainly due to two reasons. The first one pertains to 2D systems in general. In the response to long wavelength perturbations there is a spurious interaction between the system and the out-of-plane periodic images due to the $1/q^{2}$ behavior of the Fourier transform of the 3D Coulomb interaction. This effect becomes relevant for perturbations of in-plane momenta comparable or smaller than $\approx 2\pi/c$, where $c$ is the distance separating the periodic images. This range of momenta is relevant for electric transport. Additional long range interactions can arise if the slab presents a finite dipole along the z-axis.

The second reason relates to the treatment of charged systems. Charging of 2D materials is usually simulated via the use of a compensating jellium background. This approach is inappropriate as it represents a uniform doping of a 2D flake and it reproduces neither the strong voltage drop in proximity of its surfaceBrumme et al. (2014, 2015) nor the asymmetric nature of the FET configuration. This asymmetry is precisely the feature that is challenging to simulate because it breaks 3D PBC. Currently these aspects have not yet been taken into account in the calculation of vibrational properties via density functional perturbation theory.

Some methods have been proposed to deal with the FET setup at the DFT level for the calculation of total energies and forces. A dipole correction Neugebauer and Scheffler (1992); Bengtsson (1999); Meyer and Vanderbilt (2001); Brumme et al. (2014) can be used to recover 3D PBC in systems with an out-of-plane dipole moment. This method has been used to simulate chloronitrides Brumme et al. (2014) and transition-metal dichalcogenides Brumme et al. (2015) in a FET setup. Another approach to solve these issues is the effective screening medium (ESM) technique Otani and Sugino (2006). The Poisson equation is solved without 3D PBC, resulting in the correct potentials. The potentials are then modified where the electron density is negligibly small to allow their use in the KS equations with 3D PBC. This method has recently been used to simulate a graphene-based vertical field-effect tunneling transistor Wang and Cheng (2015) at the DFT level. However, linear response theory has not been developed for any of these methods. Here we solve this problem and develop a density functional perturbation theory approach tailored for 2D materials and heterostructures that includes the possibility of simulating vibrational properties in FET configuration.

We use our developments to study flexural phonons in field-effect doped graphene. In isolated graphene, those phonons do not couple to electrons due the mirror symmetry with respect to the graphene planeManes (2007). In a FET setup, however, this symmetry is broken and flexural phonons have recently been suggestedGunst et al. (2017) as a significant scattering mechanism in transport measurements. We recover the expected result that a symmetry-breaking electric field can activate a significant bare electron-phonon interaction with flexural phonons. However, we show that this interaction is strongly screened by the electrons and becomes negligible with respect to the coupling with in-plane phonons in doped graphene.

This paper is structured as follows. We first describe a model for 2D materials doped in the FET setup in Sec. II. We restrict ourselves to a description in terms of potentials. In Sec. III, we highlight the issues raised by the presence of periodic images to simulate 2D materials in the FET setup, and show how the 2D Coulomb cutoff technique can solve those issues. In Sec. IV, we detail the implementation of the 2D Coulomb cutoff for the potentials, total energy, forces, phonons and electron-phonon interactions in the Quantum ESPRESSO (QE) distributionGiannozzi et al. (2009). Finally, we exploit our implementation of the 2D Coulomb cutoff to study some properties of 2D materials specific to the FET setup. Namely, we focus on out-of-plane acoustic (ZA) phonons in a graphene FET. We show the emergence of a finite phonon frequency at $\mathbf{\Gamma}$ as well as a finite coupling to electrons for the ZA phonons, two quantities that are zero by symmetry for isolated graphene without electric field.

II Description of a 2D material doped in the FET setup

In this section we present our model for a 2D material doped in the FET setup, focusing on the potential of such a system. The central object is the 2D material itself. We consider a system with periodicity in the $\{x,y\}$ plane, defined as the infinite periodic repetition in the plane of a unit cell. The positions of the cells are $\mathbf{R}_{p}=m_{1}\mathbf{b}_{1}+m_{2}\mathbf{b}_{2}$, where $m_{1}$ and $m_{2}$ are two integers. The primitive lattice vectors $\mathbf{b}_{1},\mathbf{b}_{2}$ have coordinates in the $\{x,y\}$ plane. The $z$-component of $\mathbf{R}_{p}$ is a constant. The position of atom $a$ within the unit cell is labeled $\mathbf{d}_{a}$. The atomic internal coordinates $\mathbf{d}_{a}$ can have different $z$-components such that all atoms are not necessarily on the same plane, e.g. in the case of multilayered 2D materials. In reciprocal space, the crystal is described by reciprocal vectors $\mathbf{G}_{p}$, generated by two in-plane primitive reciprocal lattice vectors $\mathbf{b}^{*}_{1}$ and $\mathbf{b}^{*}_{2}$.

Within the DFT framework, the ground state properties of the system are determined by the ground-state electronic density $n(\mathbf{r}_{p},z)$, where we separate in-plane ( $\mathbf{r}_{p}$) and out-of-plane ($z$) space variables, as they clearly have different roles in a 2D system:

[TABLE]

The in-plane wave vector $\mathbf{k}$ and the band index $s$ define an electronic state. The Bloch wave functions $\psi_{\mathbf{k},s}$ are the solutions of the Kohn-ShamKohn and Sham (1965) (KS) equations. The KS potential of the 2D system is the sum of the external potential $V^{\rm{2D}}_{\rm{ext}}$ (which, for now, consists of the potential generated by the ions $V^{\rm{2D}}_{\rm{ion}}$), the Hartree potential $V^{\rm{2D}}_{\rm{H}}$, and the exchange-correlation potential $V^{\rm{2D}}_{\rm{XC}}(\mathbf{r}_{p},z)$:

[TABLE]

The above quantities possess the 2D-periodicity of the crystal. That is, for any $f$ 2D lattice-periodic function such as $n,V^{\rm{2D}}_{\rm{KS}},V^{\rm{2D}}_{\rm{ext}},V^{\rm{2D}}_{\rm{H}}$ or $V^{\rm{2D}}_{\rm{XC}}$, we can write

[TABLE]

The 2D Fourier transform of a 2D lattice periodic function reads

[TABLE]

where the integral is over the area of the unit cell $S$. In-plane averages are defined as $f(\mathbf{G}_{p}=0,z)=\langle f\rangle_{p}(z)$. In plane averages also extend in the out-of-plane direction. A relevant length scale for the out-of-plane extension of the 2D material would be the electronic density’s thickness t, defined such that:

[TABLE]

where $n_{0}\times S$ is the number of valence electrons per unit cell in the system, equal to the sum of the ionic charges $\sum_{a}Z_{a}$ in the neutral case. The total energy, forces, phonons and electron-phonon interactions of a neutral 2D material can be computed using the usual DFT formalismHohenberg and Kohn (1964); Kohn and Sham (1965); Ihm et al. (1979); Marting (2004); Engel and Dreizler (2011); Parr and Yang (1989) based on space integrals of products between the electronic density and various potentials. It is then sufficient to carry the out-of-plane integrals over a slab of thickness greater than $t$.

We now consider what we must do to simulate this 2D material doped in conditions emulating the FET setup. We consider a single-gate configuration, as shown in Fig. 1. A planar gate is placed parallel to 2D material and a voltage difference is applied between the two. An insulating material (gate dielectric) separates the 2D material and the gate, such that no current can flow between them and opposite surface charges accumulate on both sides. The key feature of the FET setup is its asymmetry in terms of electric field. Between the gate and the 2D material, the electric field is finite. On the other side of the 2D material, the electric field is zero. In the out-of-plane direction, it is essential that we simulate the correct 2D potentials in a region at least as large as the thickness $t$. We will not model every ions and electrons outside this region, in the gate-dielectric, substrate or gate. We rather propose ways to simulate the effects of those components on the 2D material.

The main purpose of the FET setup is to charge the 2D material. We consider an electron density such that:

[TABLE]

where $Z_{a}$ is the number of pseudo charges of atom $a$, $S$ is the surface of the 2D unit cell and $n_{\rm{dop}}\times S$ is the number of electrons added per unit cell. The total charge density of the 2D material is obtained by integrating the sum of the charge distributions associated with the ions and electrons:

[TABLE]

where the charge densities $\varrho$ are related to the corresponding electrons or ions density as $\varrho_{\rm{ion}}=en_{\rm{ion}}$ and $\varrho_{\rm{elec}}=-en$. In the FET setup of Fig. 1, the doping comes from the presence of accumulated counter-charges in the gate. We thus add a charged plane of opposite surface charge density $+en_{\rm{dop}}$ at $z_{\rm{g}}<-t/2$, playing the role of the gate:

[TABLE]

We now have a globally neutral system:

[TABLE]

The potential of the gate is included in the external potential:

[TABLE]

The resulting planar-averaged KS potential $V^{\rm{2D}}_{\rm{KS}}$ is the black line in Fig. 1, in the case of a hole-doping monolayer material. At this point, it has the features expected from a FET setup with vacuum as the insulating dielectric. The general characteristics of this potential are easily deduced from a simple parallel plate capacitor model: (i) outside the system, the electric field is zero and the potential is constant; (ii) between the 2D material and the gate, the electric field is constant and the potential is linear with a slope of $4\pi e^{2}n_{\rm{dop}}$; (iii) this electrostatic configuration translates into an out-of-plane dipolar moment which induces a shift in the KS potential:

[TABLE]

as represented in Fig. 1.

The other element to consider to have a minimal working model for the FET setup is the dielectric separating the gate and the material. Its necessity is obvious in the case of electron-doping. In that situation, the polarity of the system pictured in Fig. 1 is reversed. This means that the gate lies at a lower potential than the 2D material. In our simulations, there would then be some leaking of electrons towards the gate. This is not physical 111We could think about cold emission (an electron being emitted from a metal plate towards an other under a strong electric field), but in that case we would have to account for the work function of the gate. . In a FET setup, this is prevented by the presence of an insulator between the gate and the material. From a more mechanical point of view, the necessity of the dielectric is in fact more general. Indeed, both for hole- and electron-doping, there is an attractive force between the gate and the material, which is simply the electrostatic attraction between two oppositely charged plates:

[TABLE]

In this context, the dielectric provides a counteracting repulsive force. To emulate both the insulating and repulsive roles of the dielectric, we add a potential barrier in between the material and the gate:

[TABLE]

where $z_{\rm{g}}<z_{\rm{b}}<0$. This potential can be included in the external potential $V^{\rm{2D}}_{\rm{ext}}$. Adding such a barrier results in the potential represented by a blue line in Fig. 1. This barrier potential essentially forbids (or makes highly unlikely) the presence of electrons for $z<z_{\rm{b}}$, thus preventing electrons from leaking towards the gate. Since the electrons cannot go past the barrier, and since the ions are strongly attracted by the electrons, the barrier repulses the 2D material as a whole. As will be detailed later, the equilibrium position of the material with respect to the barrier can be determined by relaxation of the forces in the system.

In the following section, we explain how we deal with the periodic images to obtain the KS potential we just described in a plane-wave DFT code wit 3D PBC. Then, we will detail the modifications implemented to compute the total energy, forces, phonons and electron-phonon interactions for a 2D material doped in the FET setup.

III Treatment of the periodic images

Ab initio calculations based on plane-wave basis sets require periodic boundary conditions along the three dimensions (3D PBC). In this framework, periodic images of the 2D system are present in the out-of-plane direction. Our goal is for each periodic image to be strictly equivalent to the 2D system presented in the previous section, at least within a certain “physical region” around the 2D material (for example, within the boundaries of Fig. 1). In this section, we detail the issues that arise from the use of 3D PBC for the simulation of doped systems, systems with out-of-plane dipolar moment, and systems perturbed at long wavelengths. We then show how the Coulomb cutoff technique can solve those issues.

III.1 Inadequacy of 3D PBC

The 3D-periodic system obtained by adding translated copies of the 2D system generates potentials that are different from the ones described in the previous section. This comes from interactions between periodic images, due to the combination of their potentials while satisfying PBC. The sum of the KS potential from each periodic image can be written:

[TABLE]

where $i$ is an integer, and $c$ is the distance between the periodic images. $V^{\rm{2D}}_{\rm{KS}}$ is the potential of the 2D system, while $V_{\rm{KS}}$ is the one simulated in DFT with 3D PBC. In addition to the 2D PBC of Eq. (4) that $V^{\rm{2D}}_{\rm{KS}}$ already fulfills, $V_{\rm{KS}}$ has to fulfill the PBC in the third direction:

[TABLE]

We first consider a doped 2D material. Away from the direct vicinity of the material, this system behaves like a monopole, with $\lim_{|z|\to\infty}\langle V^{\rm{2D}}_{\rm{KS}}\rangle_{p}(z)=\infty$, and $V_{\rm{KS}}$ is obviously ill-defined. As mentioned before the standard method in current plane-wave DFT packages amounts to the use of a jellium background. Each slab is then globally neutral, containing the doped material and a uniform distribution of compensating charges. In between the periodic images, the resulting potential is quadratic in $z$, with extrema at mid-distance between layers, as shown in Fig. 2. This potential does fulfill the PBC and doesn’t diverge. However, it is quite different from the potential one would expect for a charged, isolated 2D material. Indeed, away from the direct vicinity of the materials, one would expect to recover a linear potential, similar to that generated by an isolated monopole.

Now we consider a 2D system with a global dipolar moment in the out-of-plane direction and a $V^{\rm{2D}}_{\rm{KS}}$ potential like in Fig. 1. Here, each periodic image is globally neutral. However, the potential $\langle V_{\rm{KS}}\rangle_{p}(z)$ would experience a shift with each periodic images, eventually diverging. Imposing PBC forbids this kind of situation. Instead, it leads to a combination of additional electric field and re-organization of the charge so that the total average electric field in one slab is zeroNeugebauer and Scheffler (1992); Bengtsson (1999); Meyer and Vanderbilt (2001). Here again, we loose the equivalence with $V^{\rm{2D}}_{\rm{KS}}$.

Finally, 3D PBC are very problematic when the system is perturbed at small wave vector. If a 2D charge density is modulated according to an in-plane wave vector $\mathbf{q}$, it generates a potential decreasing as $e^{-|\mathbf{q}||z|}$ in the out-of-plane direction. At small wave vector, the extent of the potential induced by the electron density is thus very large. When it is of the order of the distance between periodic images, there is some spurious interactions. This issue is critical when simulating the screening propertiesSohier et al. (2015) of the material as well as its response to phonon perturbationsSohier et al. (2016, 2017).

III.2 Isolate the layers with 2D Coulomb cutoff

To reduce interactions between periodic images, a naive solution is to increase the distance between them. However, in a plane-waves framework, the cost of the calculations increasing linearly with the distance, this method can be very expensive, especially in the particular cases presented before. Furthermore, this method inevitably fails for DFPT in the long wavelength limit, as spurious interactions are bound to affect the response of the material at small enough wave vectors.

One solution for the FET setup and for systems with an out-of-plane dipolar moment in general is to add a dipole correction to catch up the potential shift Neugebauer and Scheffler (1992); Bengtsson (1999); Meyer and Vanderbilt (2001); Brumme et al. (2014). However, the dipole correction has to be recalculated self-consistently at each iteration, and this method has not been extended to the DFPT framework.

Here, we tackle 3D PBC issues by using the Coulomb cutoff techniqueJarvis et al. (1997); Rozzi et al. (2006); Ismail-Beigi (2006), successfully used by other codesMarini et al. (2009); Gonze et al. (2016) in different contexts. The general idea is to cut all the potentials off between the periodic images. In effect, all physical links between periodic images are severed because the potential generated by one periodic image does not reach the others. Each slab is effectively isolated. There is no physical 3D-periodic system anymore. There is a 2D-periodic system, copied and repeated in the third dimension in order to build potentials that mathematically fulfill 3D PBC.

Each long-range potential ($V\equiv V_{\rm{ion}},V_{\rm{H}},V_{\rm{gate}}$) in the original 3D code is generated by a certain distribution of charges via the Coulomb interaction $v_{c}(\mathbf{r})=\frac{e^{2}}{|\mathbf{r}|}$. To build the corresponding cutoff potentials in the code ($\bar{V}\equiv\bar{V}_{\rm{ion}},\bar{V}_{\rm{H}},\bar{V}_{\rm{gate}}$), we use the following cutoff Coulomb interaction:

[TABLE]

where $\mathbf{r}\equiv(\mathbf{r}_{p},z)$ is a generic 3-dimensional space variable . An arbitrary charge density $\varrho$ then generates the following potential:

[TABLE]

Roughly speaking, considering a single charged plane, we generate its potential only within a certain slab of thickness $2l_{z}$ centered on the charge distribution. Within this slab, we have that $\bar{V}(\mathbf{r})=V^{\rm{2D}}(\mathbf{r})$. Outside of this slab, the potential is zero. Each periodic image of each charge distribution ($\varrho_{\rm{ion}},\varrho_{\rm{elec}},\varrho_{\rm{gate}}$) generates its own potential within its own slab. To fulfill 3D PBC, the simpler way is to cut off midway between the periodic images:

[TABLE]

Since the potentials $V^{\rm{2D}}_{\rm{ion}}$, $V^{\rm{2D}}_{\rm{H}}$, and $V^{\rm{2D}}_{\rm{gate}}$ are symmetric with respect to the plane of the associated subsystem (ions, electrons, gate), they have the same value on both sides of the corresponding slab. $\bar{V}_{\rm{ion}}$, $\bar{V}_{\rm{H}}$, and $\bar{V}_{\rm{gate}}$ are each continuous and periodic, and so is their sum $\bar{V}_{\rm{KS}}$. However, since the slabs of each subsystem do not coincide, the KS potential is only physical within the overlap of the subsystems’ slabs.

This overlap region defines a “physical region”, as illustrated in Fig. 3 , where all the potentials make sense. Outside of this region, there are some spurious unphysical variations of the KS potential. Those spurious variations are a necessary consequence of fulfilling 3D PBC. Let us consider the example Fig. 3 in more details. The simplest subsystem is the gate because $\varrho_{\rm{gate}}$ is infinitely thin in the out-of-plane direction. Within the slab $z\in[z_{\rm{g}}-\frac{c}{2};z_{\rm{g}}+\frac{c}{2}]$, we see the potential generated by the gate at $z_{\rm{g}}$.

[TABLE]

For $z$ outside of this interval, we see the potential generated by the neighboring periodic images of the gate, which has no physical sense with respect to the 2D system represented in Fig. 3. For the electrons, the charge distribution $\varrho_{\rm{elec}}$ is spread in the out-of-plane direction. Each infinitesimal slice of electronic density with surface charge density $\varrho_{\rm{elec}}(z)dz$ generates its contribution to the Hartree potential only within a certain slab. The Hartree potential is physical only within the overlap of all those slabs. If the electrons are centered around a position $z_{\rm{e}}$, that would be $z\in[z_{\rm{e}}-\frac{c}{2}+\frac{t}{2};z_{\rm{e}}+\frac{c}{2}-\frac{t}{2}]$. The ions are in a similar situation, but the charge distribution is much less spread. The difference in the spreading of $\varrho_{\rm{elec}}$ and $\varrho_{\rm{ion}}$ leads to the bumps we can observe for $\bar{V}_{\rm{ion}}+\bar{V}_{\rm{H}}$ , at $z\approx\pm 18.5$ a.u. in Fig. 3. Thus, the unphysical variations of the KS potential outside the physical region are due to the addition of potentials generated by incomplete subsystems or different periodic images.

Nevertheless, everything happens as in the previous section within the physical region associated to the KS potential. To simulate the system, we just need to make sure that the 2D material lies in this physical region. We will need 3D Fourier transforms $\bar{V}(\mathbf{G}_{p},G_{z})$, easily related to the 2D Fourier transform of $V^{2D}(\mathbf{G}_{p},z)$, Eq. 5, as follows:

[TABLE]

IV Implementation

In this section we detail the implementation of the Coulomb cutoff for the potentials, and show how it affects the energies, forces, phonons and EPC in the code. Most of the potentials, or at least their long-range parts, are calculated in reciprocal space. We thus need the Fourier transform of Eq. 18, as defined in Ref. Rozzi et al., 2006:

[TABLE]

The choice of the $\mathbf{G}=0$ value is just a convention since every potential is defined up to a constant. Here, we choose the same convention as in the original 3D code, such that the average of a potential over the unit cell is zero. A more detailed justification about this choice, especially its implications in terms of energy, can be found in App. D. For clarity, we will often need to describe the implementation of the original 3D code first in order to identify what needs to be modified. We use different notations to distinguish the potentials that stay as implemented in the original 3D code (noted $V$) and those that are modified with the implementation of the 2D Coulomb cutoff (noted $\bar{V}$). For other quantities (energy, forces, phonons and EPC), such distinction in the notation is not necessary. Indeed, their definition essentially does not change, it is the potential that is used to compute them that changes.

IV.1 KS Potential

The KS potential is the sum of the external potential, the Hartree potential, and the exchange-correlation potential:

[TABLE]

The exchange-correlation potential is short-range and does not need to be cut off. Thus, it will be ignored in the following. We note, however, that a great majority of commonly used functionals are derived in the framework of the 3D electron gas. While their relevance in the context of 2D materials is obviously questionable, the development of new functionals for 2D materials is out of scope for the present work. The implementation presented here is valid for all the usual types of pseudopotentials (norm-conserving, ultra-soft, projector-augmented wave functions). Indeed, we only modify the long-range parts of the potentials, which are independent of the pseudopotential type.

The external potential is the sum of the ionic, gate and barrier potentials:

[TABLE]

IV.1.1 Ionic Potential

The ionic potential is separated in local and non-local parts $\bar{V}_{\rm{ion}}=\bar{V}^{\rm{loc}}_{\rm{ion}}+V^{\rm{non-loc}}_{\rm{ion}}$. The non-local part is short-range. It does not need to be cut off and is ignored here. We need to compute the Fourier transform of the cutoff local potential $\bar{V}^{\rm{loc}}_{\rm{ion}}(\mathbf{G})$. It is computed from the the pseudopotentials, which are separated in short-range and long-range parts. We first describe this separation as it is done in the original code, identify what must be modified, then present the implementation of the cutoff.

In the original 3D code, the local part of the pseudopotential is a radial function in real space $v_{a}(|\mathbf{r}|)$ associated to each type of atom. It is separated in short-range (SR) and long-range (LR) parts:

[TABLE]

where $\rm{erf}(\sqrt{\eta}|\mathbf{r}|)$ is the error function with $\eta$ as a tuning parameter (see App. C for more details on that parameter). The pseudopotential $v_{a}^{\rm{SR}}(\mathbf{r})$ is indeed short-range because it always behaves as $-\frac{Z_{a}e^{2}\rm{erf}(\sqrt{\eta}|\mathbf{r}|)}{|\mathbf{r}|}$ for $|\mathbf{r}|$ large enough. In particular, we have that $v_{a}^{\rm{SR}}(\mathbf{r})=0$ for $|\mathbf{r}|\geq r_{\rm{SR}}$. The Fourier transform of the SR part is calculated via numerical integration, while the LR part is analytic. The SR part, specific to each atom, is Fourier transformed on a finite sphere:

[TABLE]

where $\Omega=S\times c$ is the volume of the unit cell. The potential $v_{a}^{\rm{SR}}$ does not need to be cut off as long as $r_{\rm{SR}}<l_{z}$, which is easily satisfied. The Fourier transform of the LR part $v_{a}^{\rm{LR}}(\mathbf{G})$ is easily found analytically, since $v_{a}^{\rm{LR}}(|\mathbf{r}|)$ is the potential generated by a Gaussian distribution of charges:

[TABLE]

This SR/LR separation is implemented in the original 3D code to enable the restriction of numerical Fourier transforms to a finite region of space. The original code also relies heavily on the rotational invariance of the radial pseudopotentials to define the arrays containing their Fourier transforms.

In our 2D implementation, we replace the analytic LR part of the pseudopotential by its cutoff version:

[TABLE]

The SR/LR separation turns out to be very convenient to implement the Coulomb cutoff. However, since the Coulomb cutoff breaks the rotational invariance, it cannot be implemented as a simple modification of the existing array. A separate array for the cutoff LR part is calculated in a separate routine. It is then added to the SR part when constructing the local part of the ionic potential:

[TABLE]

IV.1.2 Hartree Potential

The Hartree potential is relatively easy to cut off. It is computed in reciprocal space from the electronic density:

[TABLE]

IV.1.3 Gate Potential

The gate potential must be added for a doped system. In practice, the potential of the gate is added in real space to the external potential. We define directly in real-space the saw-tooth potential generated by $\varrho_{\rm{gate}}$, Eq. (9) via the cutoff Coulomb interaction Eq. (24). Within the interval $z\in[z_{\rm{g}}-\frac{c}{2},z_{\rm{g}}+\frac{c}{2}]$, it is defined as:

[TABLE]

where the constant term is due to the definition of $\bar{v}_{c}(\mathbf{G}=0)$. It sets the out-of-plane average of the potential to zero. A second gate can be added to provide more flexibility and to simulate the combination of bottom and top gate. We simply define two separate gate potentials (index ”bot” for bottom gate and ”top” for top gate):

[TABLE]

and add them together to form the total gate potential: $\bar{V}_{\rm{gate}}(z)=\bar{V}^{\rm{bot}}_{\rm{gate}}(z)+\bar{V}^{\rm{top}}_{\rm{gate}}(z)$ where the charges should be such that the whole system is neutral: $n_{\rm{top}}+n_{\rm{bot}}=n_{\rm{dop}}$.

IV.1.4 Barrier Potential

The barrier potential is necessary to relax the forces in the system and to prevent electrons from leaking towards the gate. It is also used to prevent electrons from going outside the physical region. Indeed, the variations of the potential outside the physical region can lead to the presence of potential wells. Placing a barrier potential outside the physical region ensures that no unphysical leaking occurs. The barrier is added in real space along with the gate. In practice, this barrier consists in the addition of a constant to the external potential within a certain region in the out-of-plane direction. For $z\in[-\frac{c}{2},+\frac{c}{2}]$, it is defined as:

[TABLE]

The borders of the barrier at $z_{\rm{b}1}$ and $z_{\rm{b}2}$ are smoothed via a linear transition from $V_{\rm{b}}$ to [math] on a small distance. The implementation of the gate and the barrier was adapted from a previous modification of the code, discussed in Ref. Brumme et al., 2014.

IV.1.5 Verifications

To check the consistency of our modifications on the potentials, we can first simulate the potentials of a neutral and non-polar 2D system, without gate or barrier. The corresponding ionic, Hartree and KS potentials are plotted with and without the 2D Coulomb cutoff in Fig. 4. With 3D PBC, setting the $\mathbf{G}=0$ value of the ionic or Hartree potential to zero is equivalent to the inclusion of a compensating jellium background. The potentials we observe then correspond to either ions or electrons bathed in the associated jellium. This leads to a quadratic behavior in $z$ between the periodic images.

When the 2D Coulomb cutoff is applied, we recover the linear behavior in $z$. Setting the $\mathbf{G}=0$ value of the ionic or Hartree potential to zero leads to a simple shift. For such a neutral and non-polar system, the KS potentials with and without cutoff coincide up to a constant within the physical region. This constant comes from the fact that both KS potential average to zero but the cutoff KS potential has bumps outside the physical region while the other does not.

Let us now simulate the KS potential of a hole-doped 2D material as shown in Fig. 5. Using the original code with 3D PBC, we obtain the potential of the material bathed in a jellium compensating for the added charge (or missing electrons). In that case, the KS potential is quadratic, with a varying slope and thus a varying electric field. The electric field is symmetric with respect to the plane of the 2D material. It vanishes midway between the periodic images, on the left and right borders of Fig. 5. If we use the 2D Coulomb cutoff without adding a gate, we obtain the potential that would be generated by the doped 2D material in vacuum, as in Eq. (8), within the physical region. If we add a compensating charged plane to simulate the gate, we obtain the configuration of the FET setup, with a finite electric field on the left of the 2D material, and zero electric field on the right.

Finally, we simulate the KS potential of an electron-doped system to show the necessity of the barrier potential in Fig. 6. Without the barriers, some potential wells appear on both sides of the 2D material. On the left, this is due to the presence of the positively charged gate. On the right, this is due to the unphysical variations of the KS potential outside the physical region. Electrons leak towards those potential wells, which can be inferred here from the slopes of the KS potential in the vicinity of the 2D material. Compared to what we should obtain in the FET setup, the slope of the KS potential on the left of the 2D material is too small while the slope on the right is not zero. This is due to the Hartree potential contribution from the electrons that leaked in the potential wells. This is not what we want to simulate. The addition of a potential barrier prevents the electrons from leaking towards the barrier or outside the physical region, and we find the right slopes (or electric field) in the vicinity of the material.

IV.2 Total Energy

The total energy per unit cell associated with the system is:

[TABLE]

It is the sum of the kinetic energy of the electrons, the energy of the electrons in the external potential, the Hartree energy, the exchange-correlation energy, the ion-ion interaction energy, the energy of the ions in the potential of the gate(s), and finally, the self-interaction energy of the gate(s). The terms $E_{\rm{kin}}$ and $E_{\rm{XC}}$ are short-range. They are computed as in the standard 3D code. The computation of the remaining terms is detailed in the following. A general definition of $E_{\rm{ext}}$, $E_{\rm{H}}$ and $E_{\rm{i-i}}$ can be found in App. A. Those definitions stay valid provided one uses the cutoff potentials $\bar{V}_{\rm{ext}}$, $\bar{V}_{\rm{H}}$ and $\bar{\Phi}$. Unlike the potentials, the total energy is not defined up to a constant. The choice of the $\bar{v}_{c}(\mathbf{G}=0)$ can affect the value of the energy contributions, but it should not affect the total energy. In the following, choosing $\bar{v}_{c}(\mathbf{G}=0)=0$ means that the $\mathbf{G}=0$ terms of the long range contributions to the energy will be put to zero in the 2D framework. This is justified in App. D.

IV.2.1 External Energy

The external energy $E_{\rm{ext}}$ is calculated via the sum of the eigenvalues of the KS system. Once the external potential is defined as in Eq. 26, it is used to solve the KS system and the sum of its eigenvalues includes the correct cutoff energy contribution. The external energy decomposes into the energies of the electrons in the potentials of: i) the ions $E^{\rm{ion}}_{\rm{ext}}$, ii) the gate $E^{\rm{gate}}_{\rm{ext}}$ and iii) the barrier $E^{\rm{barrier}}_{\rm{ext}}$. We use sufficiently sharp and high barrier potentials to write that $E^{\rm{barrier}}_{\rm{ext}}\approx 0$ because there are (almost) no electrons where there is a potential barrier. The other contributions $E^{\rm{ion}}_{\rm{ext}}$ and $E^{\rm{gate}}_{\rm{ext}}$ are non-zero, but we have no further modifications to make.

IV.2.2 Hartree Energy

The Hartree energy is easily written in reciprocal space as:

[TABLE]

and is computed in practice by replacing $\bar{V}_{\rm{H}}(\mathbf{G})$ by its expression (Eq. (34)):

[TABLE]

IV.2.3 Ion-ion Interaction Energy

The ion-ion interaction energy $E_{\rm{i-i}}$ is computed using the ion-ion interaction potential $\bar{\Phi}$. The computation is based on the Ewald summation technique Ewald (1921), which involves a separation into SR and LR parts $E_{\rm{i-i}}=E^{\rm{SR}}_{\rm{i-i}}+E^{\rm{LR}}_{\rm{i-i}}$. Much like for the ionic potential, we do not need to modify the SR part. Here again, we start by presenting what is done in the original code, identify what we must modify, then present the implementation of the cutoff.

In the original 3D code, following the Ewald summation technique,the ion-ion interaction potential $\Phi$ is separated in SR and LR part as follows:

[TABLE]

where the prime in the first sum excludes the case $\{\mathbf{R^{\prime}=R},a^{\prime}=a\}$ if $\mathbf{r}=\mathbf{R}+\mathbf{d}_{a}$ and $\Phi^{\rm{self}}$ subtracts that term from the second sum. The constant $\eta_{\rm{ew}}$ tunes the SR/LR separation (see App. C for more details). The SR part of the ion-ion interaction potential $\Phi^{\rm{SR}}$ is dealt with in real space and does not need to be modified as long as ${\rm{erf}}(\sqrt{\eta_{\rm{ew}}}l_{z})\approx 1$ (easily satisfied). $\Phi^{\rm{self}}$ is simply the value of $\Phi^{\rm{LR}}(\mathbf{r})$ for $\mathbf{r}-\mathbf{R}^{\prime}-\mathbf{d}_{a^{\prime}}=\bm{0}$. As such $\Phi^{\rm{self}}$ is also short-range. We include the corresponding energy contributions in $E^{\rm{SR}}_{\rm{i-i}}$. Those contributions do not need to be cut off. The contribution of the LR potential $\Phi^{\rm{LR}}$ to the energy is computed in reciprocal space and needs to be be modified.

In our implementation, we replace the Fourier transform of the LR part of the ion-ion interaction potential by its cutoff version:

[TABLE]

The LR contribution to the ion-ion interaction energy is then computed in reciprocal space as follows:

[TABLE]

where $n_{\rm{ion}}(\mathbf{G})=\frac{1}{\Omega}\sum_{a}Z_{a}e^{i\mathbf{G}\cdot\mathbf{d}_{a}}$ is the Fourier transform of the distribution of ions ($\varrho_{\rm{ion}}=en_{\rm{ion}}$).

IV.2.4 Other Energies

The other energies to account for are the energy of the ions in the potential of the gates $E_{\rm{g-i}}$, and the self interaction of the gates $E_{\rm{g-g}}$.

[TABLE]

where $d_{a,z}$ is the $z$-component of $\mathbf{d}_{a}$. We consider the most general case of a double-gate setup. Note that those contributions to the energy have a manifest dependency on the cutoff distance $l_{z}$. The total energy, of course, should not depend on $l_{z}$. As detailed in App. D, the $l_{z}$-dependent terms in the expression above will cancel with corresponding terms in $E_{\rm{ext}}$ and $E_{\rm{H}}$.

IV.2.5 Verifications

In the absence of doping, gate and barrier, the total energy is:

[TABLE]

We can thus check the consistency of the implementation of $E^{\rm{ion}}_{\rm{ext}}$, $E_{\rm{H}}$ and $E_{\rm{i-i}}$ in a neutral system. We first compute the total energy of the neutral, non-polar system of Fig. 4 with and without cutoff. We should obtain the same result as there is no issue with the periodic images in that case. We checked that the difference is below numerical precision. We can then use a neutral system with an out-of-plane dipolar moment such that interactions between periodic images do play a role without the 2D Coulomb cutoff. We use graphene with hydrogen atoms on top of half of the carbon atoms, see Fig. 7. The effect of the Coulomb cutoff is clear on both the KS potential and total energy. The KS potential of the system without cutoff illustrates the comments of Sec. III. Namely, imposing 3D PBC leads to the compensation of the out-of-plane dipolar moment by an external electric field, visible here via the finite slope of the KS potential away from the material. When we use the 2D Coulomb cutoff, we observe the right behavior, with a potential shift and no external electric field. The energy of the system simulated without cutoff tends to the one with cutoff at large distances between periodic images. With the cutoff, the energy is independent of the distance. There is a lower limit to the distance between periodic images, which is when the boundaries of the physical region are too close to the material. Still, the minimal distance we can use in our implementation of the code is negligible with respect to what we would have to use without cutoff. In the case of Fig. 7, we see that the distance between the periodic images would have to be roughly five times larger without the cutoff to obtain the same total energy as with the cutoff within $10^{-4}$ Ry. The computational cost would also be 5 times larger. A way to get the right total energy in this kind of polar material is to simulate the mirror image of the system within the unit cell. We checked that this leads to the same energy as what we find with the 2D Coulomb cutoff. Still, adding a mirror image of the system rather than the cutoff leads to a drastic increase of the computational cost.

IV.3 Forces

The forces on the ions are found by computing the derivative of the total energy with respect to a displacement $\mathbf{u}_{a,i}$ of atom $a$ in direction $i$. Only the terms related to an interaction in which the ions are involved remain. The force acting on ion $a$ in direction $i$ is written:

[TABLE]

where the first term is the force on the ion from the electrons, the second is from the other ions, the third is from the gate and the last from the barrier. The notation $\frac{\partial}{\partial\mathbf{u}_{a,i}}$ implies taking the derivative at zero displacement $\mathbf{u}_{a,i}=0$. All the quantities involved are known once the self-consistent calculation is completed. We will only detail the terms for which we need to apply the 2D Coulomb cutoff. The first term is calculated by computing the derivative of the ionic potential, separated in local and non-local parts. The derivative of the local part is found by using the Fourier transform of the pseudopotentials:

[TABLE]

The effect of the derivative in reciprocal space is to bring down a factor $-iG_{i}$ from the exponential. The corresponding force is then calculated in reciprocal space:

[TABLE]

The gate and the barrier have indirect contributions to this term. Indeed, they have an effect on $n(\mathbf{r})$, via their presence in the self-consistent KS potential.

The second term in Eq. (54) is the force from the other ions. It is found by derivation of the ion-ion interaction energy. We only treat the LR contribution, because it is the only one that needs to be cut off:

[TABLE]

The third term is the direct contribution of the electrostatic force applied by the gates to the ions. Depending on doping, it can be repulsive or attractive:

[TABLE]

where we consider the most general case of the double-gate setup.

The barrier applies no direct force on the ions $\mathbf{F}^{\rm{b-i}}_{a,i}=\mathbf{0}$. The ions effectively never see the barrier potential (they could if the barrier was smoother). However, barriers can act on the ground-state electronic density $n(\mathbf{r})$ which in turn acts on the ions. For example, if the system is too close from a barrier, the repulsive effect of the barrier will show in the self-consistent cycles, shifting the electrons away from the barrier. The first term in Eq. (54) will thus include a force that tends to push the ions away from the barrier, with the electrons. The barrier is then essential to relax the forces, arriving at an equilibrium between the attraction from the gate and repulsion from the barrier. This is illustrated by the energy and forces of the gate-barrier-material system represented in Fig. 8. For large material-barrier distances, the force tends to the attraction from the gate, and the total energy is linear. In that case, we have checked that the total force on the ions is the force between two charged plates, Eq. (14).

[TABLE]

For small material-barrier distances, as the system gets too close to the barrier, there is a sharper increase of the force and the total energy.

IV.4 Stresses for neutral 2D materials

Stresses were not implemented in the FET setup. We focus here on planar stress for neutral materials without external electric field, as it used in the relaxation of cell parameters for isolated 2D materials. The following thus applies for Cartesian coordinates $i,j\in[x,y]$. The stress is the first derivative of the energy with respect to strainNielsen and Martin (1985a, b)

[TABLE]

For the full expressions of stresses, we refer the reader to Refs. Nielsen and Martin, 1985a, b, noting that the authors’ definition of stress differs from QE’s and ours by a minus sign. We only need to cutoff the terms involving long-range potentials. The expressions of those terms in the 3D framework are given in App. A, along with a minimal outline of their derivation. They are modified as follows.

The key difference with the standard 3D code comes from the cutoff Coulomb interaction Eq. 24, and its differentiation with respect to strain:

[TABLE]

The Hartree contribution to the stress tensor reads

[TABLE]

The contribution from the long-range part of the local part of the ion-electron interaction is

[TABLE]

Finally, the long-range contribution from the ion-ion interaction gives:

[TABLE]

IV.5 Phonons and EPC

To calculate the phonon dispersion and electron-phonon interactions, as reminded in App. B, we need to compute the response of the electronic density to a phonon perturbation. In essence, the linear response of the system involves derivatives of the previous potentials and energies. Once the previous framework is set up, we just have to apply the Coulomb cutoff to the derivatives consistently. The following applies to insulators, semiconductors and metals. In the latter case, Fermi-surface effects, and notably the shift of the Fermi level that may arise with phonons at zero momentum, are treated as in the standard code, see Ref. Baroni et al., 2001.

IV.5.1 Dynamical matrix

In practice, the dynamical matrix is given by the following integrals on the unit cell:

[TABLE]

where the ”$\mathcal{V}$” notation indicates that we are using the lattice periodic part of the potential, see App. B. The first term can readily be computed from the quantities obtained in the ground-state calculation. It is computed in reciprocal space as:

[TABLE]

where $\Re(x)$ gives the real part of $x$. The last term comes from the second derivative of the ion-ion interaction $E_{\rm{i-i}}$ with respect to a phonon displacement. The contribution from $E^{\rm{SR}}_{\rm{i-i}}$ does not change. The LR part $E^{\rm{LR}}_{\rm{i-i}}$ yields the following contribution to the dynamical matrix:

[TABLE]

The second term in Eq. (70) is computed via numerical integration over the unit cell in real space. The quantities inside the integral are computed during the calculation of the electronic density response to the perturbed KS potential, presented in the following.

IV.5.2 Perturbed KS potential

The linear electronic density response is found by solving a self-consistent system involving the effective perturbation, that is the derivative of the KS potential with respect to a phonon displacement $\frac{\partial\mathcal{\bar{V}}_{\rm{KS}}(\mathbf{r}_{p},z)}{\partial\mathbf{u}_{a,i}(\mathbf{q})}$ , Eq. (108) (the notation ”$\mathcal{\bar{V}}$” indicates lattice periodic functions, see App. B). The first term is the perturbation of the external potential. The phonons only bring a direct perturbation to the potentials in which the ions are involved. This means the perturbed external potential contains only the contribution from the ionic potential. The Fourier transform of the derivative of the local part of the ionic potential has non-zero components at wave vectors $\mathbf{q}+\mathbf{G}$:

[TABLE]

where the Fourier components of the long-range part of the local pseudopotential are similar to Eq. (32):

[TABLE]

The perturbed ionic potential of Eq. (73) (along with the non-local part that is computed as in the original 3D code), is Fourier transformed and inserted in the second term of the dynamical matrix Eq. (70).

The remaining long-range potential to cut off is the Hartree potential generated by the density response, computed in reciprocal space:

[TABLE]

The density response, solution of the self-consistent system corresponding to the effective perturbation $\frac{\partial\mathcal{V}_{\rm{KS}}(\mathbf{r}_{p},z)}{\partial\mathbf{u}_{a,i}(\mathbf{q})}$ is inserted in the second term of the dynamical matrix Eq. (70).

IV.5.3 Born Effective charges and LO-TO splitting

In polar materials, the long-wavelength behavior of longitudinal optical (LO) modes depends strongly on dimensionality. Indeed, the displacement patterns of LO phonons are asociated with dipoles that interact with each others via long-range Coulomb interactions. These dipole-dipole interactions lead to an extra term in the energy of the LO mode with respect to the tranverse optical (TO) mode, thus leading to the so-called LO-TO splitting. In 2D, as shown in Ref. Sohier et al., 2017, the splitting vanishes in the zero momentum limit, but the dispersion of the LO mode displays a finite slope at the $\mathbf{\Gamma}$ point. The implementation of the 2D cutoff in DFPT as detailed above guarantees the correct treatment of the LO-TO splitting. A key quantity for this phenomenon is the tensor of Born effective charges. Notably, it gives the values of the finite slope of the LO dispersion at zero momentum. It can by computed either via the forces induced by an electric field perturbation, or via the polarization induced by atomic displacements. In both cases, the quantities involved (forces, perturbed KS potential) are already corrected as detailed above.

IV.5.4 Fourier interpolation of phonon dispersions

Dynamical matrices can be Fourier interpolatedGonze and Lee (1997); Baroni et al. (2001) to obtain phonons on dense grids at minimal computational cost. The Fourier interpolation in itself is carried out as in the standard 3D code. In polar materials, however, the LO-TO splitting corresponds to a discontinuity either in the zeroth (in 3D) or first (in 2D) order derivative of the phonon dispersion. Some non-analytic terms arise at long wavelengths due to the long- range nature of the dipole-dipole interactions. Those non-analytic terms must be modeled and excluded from the interpolation process. Since they depend on dimensionality, the interpolation requires a different treatment in 3D and 2D. The 2D treatment is implemented as detailed in Ref Sohier et al., 2017.

IV.5.5 EPC

We have all the quantities necessary to compute the EPC:

[TABLE]

The EPC matrix elements are screened via the induced part of the effective KS perturbation (Hartree and exchange-correlation). The Hartree part of the screening is then that of a 2D material. The gate and the barrier have no direct effect in the KS perturbation. Note that they are absent from this section. However, they broke the symmetry of the ground-state. In particular, the electronic distribution is not centered on the ions’ plane anymore. We will study the consequences of their presence in the following section.

V Application to Graphene FET setup

In this section we exploit our implementation of DFPT for gated 2D systems to simulate some predicted peculiarities of the FET setup. For isolated graphene without any external electric field, it can be shown that the flexural ZA phonons disperse quadratically and their energy is zero in the long-wavelength limitKatsnelson (2012). Based on the mirror symmetry with respect to the graphene plane, one can further show that ZA phonons do not couple linearly to electronsManes (2007). Those characteristics do not hold for graphene in the FET setup. First, the phonon dispersion changes due to the presence of a substrate and a gate-dielectric. Second, the presence of an electric field breaks the mirror symmetry with respect to the graphene plane, making linear coupling to electrons possible. Those FET-specific effects have not been studied in the context of DFPT. The electron-phonon coupling with flexural phonons in gated graphene was recently studied by first-principles and suggested to be a significant scattering mechanismGunst et al. (2017). However, in this work, the calculations performed do not completely include the effect of metallic screening on the electron-phonon coupling. Indeed, at the two lowest doping considered, the phonon momentum allowed for by $11\times 11$ supercell is too large with respect to the size of the Fermi surface. Furthermore, the method used in Ref. Gunst et al., 2017 assumes the Fourier transform of the derivative of the self-consistent potential to be phonon-momentum independent. We will show that this is not the case in doped graphene.

We perform DFPT calculations on graphene doped in the FET setup. We simulate the main consequences of the presence of the substrate and gate dielectric by placing two barriers at $z_{b1}$ and $z_{b2}=-z_{b1}$, such that $z_{b2}-z_{b1}\approx 5.3$ Å. Compared to graphite, it corresponds to a graphene-barrier distance that is smaller than the distance between two adjacent graphene atomic planes, but larger than the distance separating the tails of the electronic densities associated with those planes. In practice, the distances between graphene and the substrate or the dielectric depend on the details of the interactions between those materials. Here, we simply make a choice. As we will see, the results of the phonon calculations point to a rather conservative choice for the graphene-barrier distance. We plot the KS potential of the system in the upper panel of Fig. 9. We will consider three setups:

• •

“constant field”: we use two oppositely charged gates. There is a finite and constant electric field between the gates, but graphene stays neutral, i.e. $\varepsilon_{F}=0$ eV with respect to the Dirac point. In this setup, the mirror symmetry associated to the graphene plane is broken and electronic screening from graphene is minimized. Barriers are present.

• •

“one gate”: we use a single bottom gate ($z_{\rm{bot}}<0$), with a charge equal and opposite to that of the graphene sheet. The electric field is finite on the side of the bottom gate, but zero on the side of the top gate. In this setup the mirror symmetry associated to the graphene plane is also broken but the graphene is electron-doped such that $\varepsilon_{F}\approx 0.7$ eV. This implies a stronger metallic screening from $\pi^{*}$ electrons. Barriers are present.

• •

“isolated”: this setup is simply for comparison. We simulate isolated, neutral graphene, without gates and without barriers.

Linear response calculations are performed within the Quantum ESPRESSO distributionGiannozzi et al. (2009), using the 2D cutoff and FET setup implementation described in this work. We use a norm conserving pseudopotential within the local density approximation Perdew and Zunger (1981) (LDA). A dense k-point grid ($96\times 96\times 1$) is chosen to sample the Fermi surface of graphene and account for screening effects. We use a $0.01$ Ry Methfessel-Paxton smearing function for the electronic integration and a $65$ Ry kinetic energy cutoff. We use the relaxed equilibrium structure of isolated graphene in all setups. We neglect the change in lattice parameter due to doping, which was calculated to be under $0.1\%$. In the out-of-plane direction, the graphene sheet is fixed midway between the barriers, where the repulsive forces from the barriers cancel out. We neglect the comparatively small attractive forces from the charged gates.

V.1 Acoustic out-of-plane (ZA) phonons

We first show the emergence of a finite ZA phonon frequency at $\mathbf{\Gamma}$ when graphene is enclosed between two barriers. We calculate the acoustic phonons in this system and get the dispersion in the lower panel of Fig. 9. The dispersions obtained for isolated graphene (no gate, no barrier) are also shown. The dispersion of the in-plane modes is essentially unaffected by the presence of gates and barriers. In contrast, a shift in the dispersion of the ZA phonons is observed for the “one gate” and “constant field” setups. The shift being similar for both setups, it can be attributed to the presence of barriers rather than the electric field configuration. We observe in this case a rather flat dispersion, with $\omega_{\rm{ZA}}(\mathbf{\Gamma})\approx 35$ cm*-1*. When a 2D material is enclosed between two potential barriers, the ZA phonon dispersion loses its quadratic behavior in the long wavelength limit. Instead, it goes to a finite value at $\mathbf{\Gamma}$. The closer the barriers, the more confined is the 2D material and the larger is $\omega_{\rm{ZA}}(\mathbf{\Gamma})$. For graphene, a relevant reference for that value might be the ZO’ mode of graphite, in which neighboring layers are in out-of-phase ZA modes. It is often foundNicklow et al. (1972); Oshima et al. (1988); Mounet and Marzari (2005) to have a $\mathbf{\Gamma}$ frequency close to $100$ cm*-1*. The relatively small value of $\omega_{\rm{ZA}}(\mathbf{\Gamma})$ found here would thus indicate that the chosen graphene-barrier distance is rather conservative in the sense that the effect is most likely underestimated. Such a situation is preferred here, in order to find an upper bound for the strength of scattering by ZA phonons.

V.2 Gate-induced coupling to ZA phonons

We now demonstrate the emergence of a finite coupling to linear order between the electrons and out-of-plane acoustic ZA phonons, due to the electric field breaking the mirror symmetry with respect to the graphene plane. More importantly, we unravel the critical impact of screening on this coupling. We consider scattering of electrons on an iso-energetic line at $\varepsilon=0.7$ eV in the $\pi^{*}$ band. In the “one gate” setup, this corresponds to the Fermi surface of graphene. Thus, the results will be representative of the scattering involved in electronic transport. We use the same iso-energetic line in the “constant field” setup, although the line does not represent the Fermi surface since the graphene layer is neutral. In this situation, the results are not relevant for electronic transport. They correspond to the relaxation of a single electron excited at an energy of $\varepsilon=0.7$ eV. The motivation behind comparing scattering on the same iso-energetic lines is to observe the effect of electronic screening.

We fix the initial state $|\mathbf{k}\rangle$ and define a set of scattered states $|\mathbf{k}+\mathbf{q}\rangle$ on the iso-energetic line that we assume circular, as represented in Fig. 10. This implies that we neglect trigonal warping and assume elastic scattering, that is $\varepsilon_{\mathbf{k}+\mathbf{q}}=\varepsilon_{\mathbf{k}}\pm\hbar\omega_{\mathbf{q},\nu}\approx\varepsilon_{\mathbf{k}}$. We cover only half the line, the other half being equivalent by symmetry. We plot the corresponding EPC as a function of the norm of momenta. Keep in mind that increasing norm then corresponds to increasing values of the angle between initial and scattered states $\theta_{\mathbf{k}+\mathbf{q},\mathbf{k}}$.

We plot the corresponding values of $g_{\nu}(\mathbf{q})=g_{\mathbf{k},\pi^{*},\mathbf{k}+\mathbf{q},\pi^{*},\nu}$ for $\nu=$ZA, TA, LA in the upper panel of Fig. 10. The in-plane acoustic modes are plotted for comparison. We clearly observe the emergence of a non-zero value for $g_{\rm{ZA}}$, as expected. In both neutral and doped cases, the value of the coupling vanishes for the large phonon momenta corresponding to backscattering. In the neutral case, $g_{\rm{ZA}}$ is comparable to the coupling with in-plane acoustic phonons. However, while the coupling to in-plane phonons is essentially doping-independent, $g_{\rm{ZA}}$ is much smaller in doped graphene. Indeed, the coupling to in-plane phonon is dominated by gauge fieldsSohier et al. (2014), which do not affect the local charge density and thus are not screened. In contrast, the gate-induced coupling to ZA phonon acts as a deformation potential, that is a periodic modulation of charge density in the underlying lattice potential in which the free electrons move. In terms of Hamiltonian, the perturbation is diagonal and proportional to identity in the Dirac spinor basis. As such, this perturbation is screened by graphene’s electrons. In the doped case, the gate-induced coupling $g_{\rm{ZA}}$ is strongly screened by metallic graphene and it becomes negligible. In Eq. (2) of Ref. Gunst et al., 2017, the authors considered a deformation potential (called ”field-induced coupling constant”) that depends on gate-voltage but is independent of momentum $\mathbf{q}$. The screening from the conduction electrons of graphene does in fact bring a strong momentum dependencySohier et al. (2015) to this quantity.

In Ref. Gunst et al., 2017, it is argued that despite relatively weak coupling, the high occupation of ZA phonons lead to a considerable scattering probability. In the lower panel of Fig. 10, we study the transport-relevant quantity

[TABLE]

where $N_{\mathbf{q},\nu}$ is the phonon occupation at room temperature (Böse-Einstein distribution with $T=300$ K) and the angular term conveys the fact that backscattering is more detrimental to electronic transport. In the framework of the relaxation time approximation and elastic processes, the integral of this term over the Fermi surface corresponds to the scattering rate. In the neutral case, we see in the lower panel of Fig. 10 that despite vanishing backscattering and the relatively small coupling overall, $\mathbf{\Gamma}_{\mathbf{k},\mathbf{k}+\mathbf{q},\rm{ZA}}$ is comparable to the other acoustic mode, thanks to a relatively large occupation of ZA phonons (note that as mentioned before, the phonon frequency is probably a lower bound so the occupation and the scattering rate are upper bounds). This makes field-induced scattering by ZA phonons potentially important for carrier relaxation in neutral graphene under a constant electric field. In the more relevant case of single-gated and doped graphene, the scattering from the ZA phonon is screened and negligible. The doping level considered here is rather large. However, for lower, more experimentally realistic doping levels, the coupling would be similar or smaller. Indeed, the electric field and thus the field-induced bare coupling would be smaller. Since the screening scales with the dimension of the Fermi surface Sohier et al. (2015), similar screening would be obtained for scattering around the Fermi surface at any doping. This general trend is verified in Fig. 11, where we use two different doping and compare the electron-phonon coupling as a function of the momentum rescaled by the size of the Fermi surface $|\mathbf{q}|/2k_{F}$.

Finally, in addition to screening effects, the flatness of the ZA dispersion also plays a role in decreasing the coupling with respect to the quadratic dispersion of isolated graphene (see definition Eq. 77). In experimental setups, the ZA dispersion will depend on the details of the interactions with the subtrate and the gate dielectric. In our simulations, the barriers act as an approximation to these interactions. However, as mentionned above, we expect experimental values of $\omega_{\rm{ZA}}(\mathbf{\Gamma})$ to be comparable or larger than the simulated value. Thus, we expect the values of the coupling found in the highly doped case to be an upper bound for scattering around the Fermi surface of graphene at any finite doping level. This coupling would thus be undetectable in transport measurements. Those results are in stark contrast with the recent first-principles study of this effectGunst et al. (2017). The main reasons for this discrepency are the dispersion of flexural phonons and electronic screening. Due to the flat dispersion of the flexural phonons, we find overall smaller couplings to flexural phonons than in Ref. Gunst et al., 2017. For the relaxation of photoexcited carriers in neutral graphene (”constant field” setup), we find that the coupling is still comparable to the coupling with in-plane phonons. For carrier transport in doped graphene, however, the coupling to flexural phonons is screened and becomes negligible. Electronic screening plays a key role in electron scattering. It is a highly dimensionality-dependent quantity that requires the correct 2D framework. It can be quite difficult to model and predict its effect on complex mechanisms like electron-phonon coupling, even more so in complex systems like gated 2D heterostructures. As illustrated here, adequate DFPT methods are then an invaluable tool.

VI Conclusion

Manipulating the electronic properties of heterostructures via the field effect is key to many future usage of 2D materials. We first set the framework for the simulation of charged heterostructures within the field-effect setup. We then show that various issues arise within the standard three-dimensional periodic boundary conditions for those systems, and that the truncation of the Coulomb interaction in the out-of-plane direction is a simple and efficient solution. We detail the implementation of the two-dimensional Coulomb cutoff and the field effect setup within the Quantum ESPRESSO distribution, for ground state and linear response calculations. The most basic changes concern the construction of potentials equivalent to the those generated by an isolated two-dimensional system. Changes are then made accordingly throughout the code to compute physical quantities properly defined in the two-dimensional framework. This includes total energies, forces, stresses, phonons and electron-phonon interactions. We demonstrate the relevance of the implementation by studying flexural (or out-of-plane acoustic) phonons for graphene in the field effect setup. Our results show the emergence of a finite phonon frequency in the long wavelength limit, as well as a finite coupling to electrons. However, electronic screening makes the coupling to flexural phonons negligible with respect to the coupling to in-plane phonons. This implies that the phenomenon is undetectable in transport measurements for graphene at finite doping.

VII acknowledgements

The authors acknowledge support from the European Union Horizon 2020 research and innovation program under Grant agreement No. 696656-GrapheneCore1 and from Agence Nationale de la Recherche under the reference n. ANR-13-IS10-0003-01. Part of the calculations were carried out at IDRIS, CINES, CEA and BSC TGCC.

Appendix A DFT

In this section we introduce the quantities that can be calculated in DFT. We give the most straightforward expressions for a generic system in 3D space. The formulas might be applied to 3D-periodic or 2D-periodic materials. The following description is obviously very far from exhaustive. The aim is simply to set the notations and provide definitions for the quantities mentioned in the main text. For more details, we refer the reader to the literature, for example Refs. Marting, 2004; Engel and Dreizler, 2011; Parr and Yang, 1989. The lattice vectors are noted $\mathbf{R}$, the internal coordinate of atom $a$ is $\mathbf{d}_{a}$. The reciprocal lattice vectors $\mathbf{G}$ are such that $e^{i\mathbf{G}\cdot\mathbf{R}}=1$.

A.1 Potentials

The central potential in the self-consistent process is the Kohn-ShamKohn and Sham (1965) potential. It contains:

• •

$V_{\rm{ext}}$: the principal contribution to the external potential $V_{\rm{ext}}$ is the potential generated by the ions $V_{\rm{ion}}$ calculated via the pseudopotentials. For the purpose of this paper, we consider only the local part of $V_{\rm{ion}}$, written as:

[TABLE]

where $v_{a}$ is the pseudopotential associated to atom $a$. We can put other contributions into $V_{\rm{ext}}$, but in this appendix $V_{\rm{ext}}=V_{\rm{ion}}$.

• •

$V_{\rm{H}}$: the Hartree potential is given by:

[TABLE]

Here, and whenever there is no specified interval, the integrals are carried out over the entire space spanned by the corresponding variable.

• •

$V_{\rm{XC}}$: the exchange-correlation potential is based on the local density approximation Perdew and Zunger (1981) (LDA).

We also have to mention the potential $\Phi(\mathbf{r})$ which is the potential generated by the ions used to compute ion-ion interactions. Although it is generated by the same source as $V_{\rm{ion}}$, it is not calculated via the pseudopotentials. To generate this potential, the ions are modeled by a collection of point charges. This potential is defined as:

[TABLE]

where the prime on the sum excludes the $\{\mathbf{R}=\mathbf{R}^{\prime},a=a^{\prime}\}$ case if $\mathbf{r}=\mathbf{R}+\mathbf{d}_{a}$.

A.2 Total Energy

What we call the total energy of the system is the clamped-ions energy or the Born-Oppenheimer energy surface Born and Oppenheimer (1927):

[TABLE]

where

[TABLE]

and $\mathcal{E}_{\rm{XC}}$ is the exchange-correlation energy. For any of the above quantity $\mathcal{E}$, one can define the corresponding energy per unit cell $E=\mathcal{E}/N$ where $N$ is the number of unit cells. This is often more useful in practice, since the system is infinite, as are the energies $\mathcal{E}$.

A.3 Forces

To calculate the force on atom $a$ in direction $i$, we compute the derivative of the total energy per unit cell with respect to a displacement $\mathbf{u}_{a,i}$ of this atom in this direction, and take the value at $\mathbf{u}_{a,i}=0$. Using the Hellmann-Feynman theorem Feynman (1939); Hellmann (1937), the force acting on ion $a$, in direction $i$ is given by:

[TABLE]

where the integral is carried over the volume of the unit cell $\Omega$. Here and in the following, the notation $\frac{\partial}{\partial\mathbf{u}_{a,i}}$ represents the value of the derivative at zero displacement. The first term is the contribution from the electrons, the second from the ions. The forces can be computed as soon as we have solved the ground state, since there are only known quantities and their derivatives.

A.4 Stresses

The stress is the first derivative of the energy with respect to the strain tensor $\overleftrightarrow{\epsilon}$ Nielsen and Martin (1985a, b):

[TABLE]

In practice, it is derived from the total energy by applying the scaling procedure $\mathbf{r}^{\prime}=(1+\overleftrightarrow{\epsilon})\mathbf{r}$. We will focus here on the contributions from long-range potentials, the only ones involved in the 2D Coulomb cutoff process. We give an outline of the general derivation, below is the treatment of the Hartree contribution:

[TABLE]

Since $\Omega n(\mathbf{G})$ is the number of electrons per unit cell, it is invariant under strain. We are thus left with the derivative of the Hartree potential:

[TABLE]

Starting from the fact that $\Omega n(\mathbf{G})$ is invariant and knowing that the volume transforms as $\Omega^{\prime}\approx(1+\sum_{i}\epsilon_{ii})\Omega$, one finds that $\frac{\partial n(\mathbf{G})}{\partial\epsilon_{ij}}=-\delta_{ij}n(\mathbf{G})$. We are left with the derivative of the Coulomb interaction. Since the scaling procedure gives $\mathbf{G}^{\prime}=(1-\overleftrightarrow{\epsilon})\mathbf{G}$ in reciprocal space, we have that $\frac{\partial G_{l}}{\partial\epsilon_{ij}}=-\delta_{li}G_{j}$. The derivative of the Coulomb interaction is then computed using the chain rule:

[TABLE]

with the 3D Coulomb interaction as $v_{c}(\mathbf{G})=\frac{4\pi e^{2}}{\mathbf{G}^{2}}$. The Hartree contribution is thus:

[TABLE]

Similar derivations yield the contribution from the long-range part of the local ionic potential. By differentiation of (see also Eq. 27):

[TABLE]

one gets:

[TABLE]

Finally, the long-range part of the ion-ion interaction is written:

[TABLE]

and gives the following stress:

[TABLE]

Appendix B DFPT

DFPT enables the computation of the linear response of the ground state to given perturbation. We focus here on phonon perturbations, as implemented in the phonon code of the Quantum ESPRESSO distribution. Again, we only give minimal description to settle the notation and define the quantities referred to in the main text. The reader might refer to the literature for more details, for example Refs. Gonze and Lee, 1997; Baroni et al., 2001.

B.1 Phonons

A phonon perturbation of momentum $\mathbf{q}$ is represented by a collection of displacements $\mathbf{u}_{a,i}(\mathbf{R})$ of atom $a$ in Cartesian direction $i$:

[TABLE]

where $\mathbf{u}_{a,i}(\mathbf{q})$ is the Fourier transform of $\mathbf{u}_{a,i}(\mathbf{R})$. The phonon frequencies are obtained from the second derivative of the total energy of the crystal $\mathcal{E}_{\rm{tot}}$ (not the energy of a unit cell) via the matrix of the interatomic force constants defined as Giannozzi et al. (1991); Gonze and Lee (1997); Baroni et al. (2001):

[TABLE]

In this particular context, it does not make sense to talk about energy per unit cell. Indeed, the energy is not lattice periodic because of the phonon perturbation. There are two contributions, one from the electrons, one from the ions:

[TABLE]

The first term of Eq. (103) and Eq. (105) are simply the second derivatives of quantities already computed in DFT. The second term of Eq. (103), however, contains the linear response of the electronic density to a phonon perturbation. This quantity can be calculated within DFPT. A phonon perturbation translates into a periodic perturbation of the potential generated by the ions, that is a periodic perturbation of $V_{\rm{ext}}$:

[TABLE]

where we now work with the (single-component) Fourier transform of the phonon perturbation $\mathbf{u}_{a,i}(\mathbf{q})$. The phonon perturbation triggers the linear response of the electronic density:

[TABLE]

which is found by solving a new set of equations, involving the linear perturbation to the KS potential:

[TABLE]

where the $\mathcal{V}$ notation indicates that we take the lattice-periodic part of the perturbations:

[TABLE]

The last two terms of Eq. (108) are generated by the density response Eq. (107). We thus have a new self-consistent system to solve. Once self-consistency is reached, we can calculate the dynamical matrix $\mathcal{D}$ which is the Fourier transform of the matrix of the force constants:

[TABLE]

where $M_{a}$ is the mass of atom $a$, and we have used translational invariance to express the matrix of the force constant as a function of the generic lattice vector $\mathbf{R}$. The eigenvalue problem:

[TABLE]

gives the frequencies $\omega_{\mathbf{q},\nu}$ ($\omega^{2}_{\mathbf{q},\nu}$ being the eigenvalues) and eigenvectors $\mathbf{e}_{\mathbf{q},\nu}$ of mode $\nu$ at momentum $\mathbf{q}$.

B.2 EPC

The electron-phonon interaction matrix elements are obtained from the derivative of the KS potential as follows:

[TABLE]

Appendix C Long-range / short range separations

There are two short-range / long-range (SR/LR) separations performed in the code. One for the pseudopotentials (Eq. 27) and one for the computation of the ion-ion interaction, Eq. 42. They are done for different reasons. The first is done to enable the computation of the Fourier transform of the pseudopotentials, the second is done to optimize the convergence of the real/reciprocal space computation of the Ewald sums. To each is associated a tuning parameter, $\eta$ or $\eta_{\rm{ew}}$. Since in both case, the long-range/short-range contributions are always put back together before computing any physical quantity, the tuning parameters can be chosen separately. As in the original code, we set $\eta=1$ for the pseudopotentials while the Ewald splitting parameter $\eta_{\rm{ew}}$ used in the computation of the ion-ion interaction is chosen depending on the plane-wave cutoff.

Appendix D Treatment of the $\mathbf{G}=0$ singularities in the 2D code

D.1 $\mathbf{G}=0$ value of the Coulomb interaction

The treatment of the $\mathbf{G}=0$ terms with a 2D Coulomb cutoff is developed in Ref. Rozzi et al., 2006, where the authors show that one should use the following value for the $\mathbf{G}\to 0$ limit of the truncated Coulomb interaction:

[TABLE]

In the following, we quickly justify this recommendation. Note that we use the notation $\mathbf{G}\to 0$ to distinguish this value from the value $\bar{v}_{c}(\mathbf{G}=0)=0$ used in our implementation. The potential $V(\mathbf{r}_{p},z)$ generated by a generic 2D distribution $m(\mathbf{r}_{p},z)$ via the cutoff Coulomb interaction is written:

[TABLE]

By changing variables and exploiting the in-plane periodicity of $m(\mathbf{r}_{p},z)$, it can be shown that the planar average of the potential $V(\mathbf{G}_{p}=0,z)$ reads:

[TABLE]

This can be written as :

[TABLE]

with

[TABLE]

The first term of the above equation is the one that gives the diverging behavior in the potential of a charged plane. However, this term vanishes as soon as the 2D system is globally neutral within the cutoff because it does not depend on $z$. If we replace $m$ by a globally neutral distribution $n_{\rm{tot}}$ that would be the sum of the distributions of the electrons, ions and gate, we get:

[TABLE]

We can thus drop this term. The definition of the $\bar{v}_{c}(\mathbf{G}\to 0)$ is then found by Fourier transform of the remaining term along the third direction (as in Eq. (23)):

[TABLE]

D.2 Implementation

We now show why we can further simplify the process and use $\bar{v}_{c}(\mathbf{G}=0)=0$ in our implementation. The exchange-correlation and barrier contributions to the potentials and energies are ignored here because they bring no divergence. Since we are considering cutoff quantities, the following concerns the long-range part of the potentials and the corresponding contributions to energy when the long-range/short-range separation is done. In order to simplify the argument, we do not make the distinction in the notation. In the following, the ”tilde” quantities are those defined using the value $\bar{v}_{c}(\mathbf{G}\to 0)=-2e^{2}\pi l_{z}^{2}$ recommended in Ref. Rozzi et al. (2006). Here is how we define the potentials in our implementation:

[TABLE]

Defined this way, the $\mathbf{G}=0$ value of our potentials is zero (at least for the long-range part in the case of $\bar{V}_{\rm{ion}}$ and $\bar{\Phi}$). Note that if we sum $\bar{V}_{\rm{H}}$, $\bar{V}_{\rm{ion}}$ and $\bar{V}_{\rm{gate}}$, we find that $\bar{V}_{\rm{KS}}=\tilde{V}_{\rm{KS}}$, which is essential. The potentials give the following energies:

[TABLE]

where, once again, we have that all the $\mathbf{G}=0$ contributions to the energy are zero, and that $E_{\rm{tot}}=\tilde{E}_{\rm{tot}}$, if we sum all the contributions. The process described above is equivalent to setting $\bar{v}_{c}(\mathbf{G}=0)=0$, and it gives the same KS potential and total energy as using $\bar{v}_{c}(\mathbf{G}\to 0)=-2\pi e^{2}l_{z}^{2}$. It is also very close to the process used in the original 3D code, which allows us to minimize the changes.

The $\mathbf{G}=0$ values of the cutoff potentials can thus be set to zero. This is pretty straightforward to apply when the whole potential is cutoff, without prior SR/LR separation, that is in the case of $\bar{V}_{\rm{H}}$ and $\bar{V}_{\rm{gate}}$. When the SR/LR separation is done, for $\bar{V}_{\rm{ion}}$ and $\bar{\Phi}$, it becomes more subtle, because the short-range parts are not cutoff and the corresponding $\mathbf{G}=0$ term are finite. Let’s examine what must be done, always trying to minimize the changes with respect to the 3D code. The ionic potential’s $\mathbf{G}=0$ term in the 3D code is

[TABLE]

This term is non-divergent, non-zero and short-range. It is computed by numerical integration in a finite sphere. It can be referred to as the ”$\alpha$” term, leading to the so-called ”$\alpha$Z” energy contribution. It is a combination of two divergent terms. Following the overall strategy of the 2D cutoff implementation, the first term should be separated in SR/LR part and then the LR part replaced by its cut off counterpart. The second term would be directly replaced by its cut off counterpart. Let us follow this process and determine how the ”$\alpha$” term should be corrected in 2D. We first separate the 3D ionic potential in SR and LR parts. The SR part depends on the pseudopotential and it is left unchanged. The LR part of the 3D potential $\tilde{V}^{3D,LR}_{\rm{ion}}(\mathbf{G}=0)$ can be separated in two terms: (i) a divergent part and (ii) a finite, nonsingular partMartyna and Tuckerman (1999)

[TABLE]

The divergent part (i) cancels out with $v_{c}(\mathbf{G}\to 0)n_{\rm{ion}}(\mathbf{G}=0)$ in 3D. It’s cut off counterpart similarly cancels with $\bar{v}_{c}(\mathbf{G}\to 0)n_{\rm{ion}}(\mathbf{G}=0)$ in 2D. The changes to be made to the ”$\alpha$” term would thus be subtracting and adding zero, i.e, doing nothing. The term (ii) in its 3D form, $\tilde{V}^{3D,LR,NS}_{\rm{ion}}(0)$ should be subtracted and replaced by the $\mathbf{G}=0$ value of its cut off counterpart. The latter being zero, the only correction to make to the 3D ”$\alpha$” term is to subtract $\frac{\pi e^{2}}{\eta\Omega}\sum_{a}Z_{a}$. In the 3D Ewald summation, the $\mathbf{G}=0$ term is:

[TABLE]

Similarly to the ”$\alpha$” term, the only thing to do is to subtract $-\frac{\pi e^{2}}{\eta_{\rm{ew}}\Omega}\sum_{a}Z_{a}$. Note that it looks like we are adding some $\eta$-dependent terms to $E^{\rm{ion}}_{\rm{ext}}$ and $E_{\rm{i-i}}$. This is actually not the case. The $\eta$-dependent terms we add in the $\mathbf{G}=0$ terms cancel out with equal and opposite contributions in the $\mathbf{G}\neq 0$ terms. We are defining the $\mathbf{G}=0$ term that corresponds to the cutoff potential used for the $\mathbf{G}\neq 0$ terms. $E^{\rm{ion}}_{\rm{ext}}$ and $E_{\rm{i-i}}$ are independent of $\eta$ or $\eta_{\rm{ew}}$, and those parameters can be still be chosen independently for the two SR/LR separation processes as said in App. C.