# Breakdown of optical phonons’ splitting in two-dimensional materials

###### Abstract

We investigate the long-wavelength dispersion of longitudinal and transverse optical phonon modes in polar two-dimensional materials, multilayers, and their heterostructures. Using analytical models and density-functional perturbation theory in a two-dimensional framework, we show that, at variance with the three-dimensional case, these modes are degenerate at the zone center but the macroscopic electric field associated with the longitudinal-optical modes gives rise to a finite slope at the zone center in their corresponding phonon dispersions. This slope increases linearly with the number of layers and it is determined solely by the Born effective charges of the material and the dielectric properties of the surrounding media. Screening from the environment can greatly reduce the slope splitting between the longitudinal and transverse optical modes and can be seen in the experimentally relevant case of boron nitride-graphene heterostructures. As the phonon momentum increases, the intrinsic screening properties of the two-dimensional material dictate the transition to a momentum-independent splitting similar to that of three-dimensional materials. These considerations are essential to understand electrical transport and optical coupling in two-dimensional systems.

###### pacs:

63.22.-m, 73.63.-b, 77.22.EjVan der Waals heterostructures Geim and Grigorieva (2013) will assuredly play a key role in future electronic and optoelectronic devices Fiori et al. (2014); Jariwala et al. (2014); Wang et al. (2012); Xia et al. (2014). Many of the potential candidates to build those next-generation devices are polar two-dimensional (2D) materials, including transition-metal dichalcogenides (TMDs) and hexagonal boron nitride (h-BN). A significant consequence of polarity in these materials is the generation of long-ranged electric fields by the polarization density associated with their longitudinal optical (LO) phonon modes. These fields strongly influence the transport properties of the monolayers as well as their heterostructures, via the Fröhlich electron-phonon interaction Verdi and Giustino (2015); Sjakste et al. (2015); Sohier et al. (2016). Importantly, they lead to additional dipole-dipole interaction terms affecting the dispersion characteristics of optical phonons. This leads to a splitting between the LO and transverse optical (TO) modes, driven by the long-ranged Coulomb interactions and electronic screening. Here, we show in detail how this LO-TO splitting is drastically affected by dimensionality. In 3D, it is well understood how this splitting is independent of the norm of the phonon momentum and how it lifts the degeneracy of the LO and TO phonon modes in the long wavelength limit. We show here that in 2D, the splitting depends on the norm of the phonon momentum, vanishes at the zone center and leads to a discontinuity in the derivative of the LO phonon dispersion. First-principles computations of this phenomenon are delicate independently of dimensionality, due to the long-ranged nature of the dipole-dipole interactions leading to electric fields that are not compatible with periodic boundary conditions. In 3D materials, however, the development of analytical models based on the Born effective charges and the dielectric tensor has allowed for the correct numerical treatment of the LO-TO splittingKunc (1982); Giannozzi et al. (1991); Baroni et al. (2001); Gonze et al. (1994); Gonze and Lee (1997); Kern et al. (1999). Similar efforts remain to be attained in 2D materials. A few theoretical and computational works have identified the main characteristics of the LO-TO splitting in 2D h-BN. Namely, it was pointed out that the splitting vanishes at the zone center Mele and Král (2002) while the slope of the LO dispersion is finite Sánchez-Portal and Hernández (2002); Michel and Verberck (2009, 2011). However, the peculiarity of screening in a 2D framework was not accounted for, and the proposed slopes were not quantitatively accurate. Furthermore, the empirical force constant models used lacked the generality and predictive power of first-principles simulations.

It should be pointed out that in most first-principles calculations of two-dimensional materials, the 2D system is repeated periodically while using a large interlayer distance. In such a setup, the splitting does not vanish at the zone center, due to the spurious long-range interactions with the periodic copies. The equality of LO and TO frequencies at zone center can be enforced by simply omitting the 3D splitting Zólyomi et al. (2014), but the dispersion at small momenta remains erroneous due to spurious interactions between the periodic images. Alternatively, the authors of Ref. Wirtz et al., 2003 study zone-border phonons so that LO phonon modes in neighbouring layers are out-of-phase. In this case, the electric fields generated by periodic images cancel out at long wavelengths. However, this approach also removes the effects of the electric field generated by a single isolated layer, and so it is not completely satisfactory. In this work we study in detail the LO-TO dispersions in 2D materials, that is their slope splitting at the zone center as well as the transition to a flatter dispersion for the LO mode at larger momenta. This is accomplished by formulating a detailed and quantitative analytical model that takes into account the subtle effects of screening in 2D, and by performing density-functional perturbation theory (DFPT) calculations of the LO-TO splitting in the appropriate open boundary 2D framework, using our implementation Sohier (2015); Sohier et al. (2016) of Coulomb cutoff techniques Ismail-Beigi (2006); Rozzi et al. (2006) in the relevant codes (PWSCF and PHONON) of the Quantum ESPRESSO distribution Giannozzi et al. (2009); Baroni et al. (2001). This latter approach unlocks the full potential and versatility of DFPT methods for 2D systems, and allows us to study LO-TO splitting in different monolayers, multilayers and heterostructures, providing the insight needed to interpret the results of various experimental setups and their effect on transport measurements.

In order to frame the discussion, we begin with introducing our
analytical model. We are interested in the dispersions of LO and TO phonons
near the point in 2D materials or layered 3D materials,
and we consider a phonon of in-plane momentum
. For simplicity, we use the
limit of the phonon displacements
and neglect the deviation from the strictly longitudinal
and transverse nature of the phonon modes as
momentum increases. The displacement of atom in the unit cell
is then given by ,
where is the mass of atom and is
the limit of the eigenvector of the dynamical matrix
corresponding to the LO mode, normalized over the unit cell.
The origin of the LO-TO splitting in polar materials
is the polarization density generated by the atomic
displacement pattern and
the associated electric fields.
The Fourier transform ^{1}^{1}1See supporting information for the
definition of the Fourier transform depending on the
dimensionality of the system of the polarization density is

(1) |

where is the elementary charge and is either the volume of a 3D-periodic system’s unit cell or the area of a 2D-periodic system’s unit cell . The tensor of Born effective charges associated with atom in the unit-cell is . Note that in the general case, treated in Appendix, it is not necessarily possible to define a pair of LO/TO modes belonging to the same irreducible representation at . For simplicity and clarity, we focus here on materials where this is possible, such as the commonly studied materials with hexagonal in-plane symmetry. Further, we assume in-plane isotropy with respect to long-wavelength perturbations, strictly in-plane phonon displacements and diagonal Born-effective-charges tensors. Within these assumptions, valid for all the materials mentioned in this work, the LO and TO modes would be mechanically similar in the long-wavelength limit. In polar materials, however, the emergence of long-ranged dipole-dipole interactions differentiates them. The divergence of the polarization, with Fourier transform , represents a polarization charge density. This polarization charge density is zero for the TO mode due to the orthogonality of the polarization and the direction of propagation. The LO mode, however, does induce a polarization charge density which, in turn, generates an electric field. The Born effective charges describe the atomic response to such a field and this leads to an additional restoring force on the atoms and an increase in energy cost for the displacement pattern of the LO mode with respect to the TO mode. The resulting relationship between the square of the frequencies of these modes is well-known in 3D bulk materials Baroni et al. (2001); Giannozzi et al. (1991); Gonze and Lee (1997). As the central analytical result of this paper, we generalize this relationship for both 2D and 3D layered materials as

(2) |

where and is a screened Coulomb interaction discussed below. In general, the second term on the right-hand side depends on the momentum direction via the Born effective charges and the screening. Here, we focus on in-plane isotropic materials for which all quantities depend only on the magnitude of for . We capture the fundamental role of dimensionality in the long wavelength limit by introducing the following simple model for the screened macroscopic Coulomb interaction

(3) |

where, in addition to the different powers at which appears in the Coulomb interaction, screening is considered differently in 2D and 3D. In the 3D layered materials discussed here, screening can be described by the in-plane dielectric constant of the bulk . The splitting (second term in Eq (2)) is then independent of momentum, with an expression that depends notably on the effective charges and the dielectric constant. In 2D, screening can be described by Keldysh (1979); Wehling et al. (2011); Cudazzo et al. (2011); Berkelbach et al. (2013); Steinhoff et al. (2014); Sohier et al. (2016), where the constant describes the dielectric properties of the environment. In the case of two semi-infinite dielectrics on each side of the 2D material, with relative permittivity and , we have . In the case of an isolated monolayer . The effective screening length describes the screening properties of the 2D material itself. It can be approximated as where is the thickness of the 2D material (see Ref. Sohier et al., 2016 for details). Due to the effect of dimensionality on the screened Coulomb interaction, the 2D splitting now depends on momentum. Namely, Eq. (2) can be written in the shorthand form where is a constant depending on the effective charges and the masses. For , the splitting is linear in and screened solely by the surrounding medium. As momentum decreases, the electric field lines associated with the polarization charge density spread more and more in the surrounding medium, leading to a vanishing dipole-dipole interaction and thus a vanishing splitting. At , although the splitting is zero, the slope of the LO dispersion is finite and discontinuous. It has a positive value in all directions. For , one recovers the 3D case. Indeed, in this limit, the electric field associated with the polarization density is confined within the thickness of the monolayer. The dipoles interact within the layer and it makes no difference whether the monolayer is isolated or surrounded by other monolayers (as in a layered 3D material). The material-specific effective screening length determines the transition between the two regimes, and can be estimated from first-principles Wehling et al. (2011); Cudazzo et al. (2011); Berkelbach et al. (2013); Steinhoff et al. (2014); Sohier et al. (2016).

Monolayer | (eV Å ) | (Å) | (cm) |
---|---|---|---|

h-BN | |||

MoS | |||

MoSe | |||

MoTe | |||

WS | |||

WSe | |||

InS | |||

InSe |

To complement these analytical results we performed DFPT calculations in 2D open boundary conditions of the of LO and TO dispersions in a selection of monolayers, as shown in Fig. 1. These are in excellent agreement with the 2D analytical model of Eq. (2), using the parameters obtained independently from first principles and reported in Table 1. Note that at larger momenta the LO dispersion displays some material-specific behaviour, as seen more clearly for MoS. Contrary to the 3D case, the computation of the frequency of the LO mode is not problematic at exactly, since the splitting is zero. To obtain the correct behaviour at small but finite momentum, however, several issues arise. In standard plane-wave DFPT, spurious interactions with the artificial periodic images yield erroneous results in the long-wavelength limit: when the phonon momentum is smaller than the inverse of the distance between periodic images, the atoms of one monolayer feel the polarization field of the neighbouring monolayers. In addition, the periodic images lead to spurious screeningKozinsky and Marzari (2006); Sohier et al. (2015). No matter the amount of vacuum inserted in between the periodic images, for sufficiently small momenta one will always end up with the response of a 3D-periodic system, i.e. a non-vanishing splitting, as shown in Fig. 2. This issue is solved by using our implementation of the 2D Coulomb cutoff technique Sohier (2015); Sohier et al. (2016), as shown in Fig. 1. In addition, a technical difficulty arises when using Fourier interpolation schemes Gonze and Lee (1997); Baroni et al. (2001); Giannozzi et al. (1991), that are otherwise very useful to obtain phonon dispersions on fine momentum grids at minimal computation cost.

In fact, the discontinuity in the first derivative of the LO dispersion, due to the long-ranged dipole-dipole interactions, leads to interatomic force constants (IFC) decaying with power-law in real-space Piscanec et al. (2004). Since the Fourier interpolation scheme relies on finite-ranged IFCs, it is not suited to capture correctly LO-TO splitting. A similar issue is present in 3D, with a discontinuity in the value of the frequency rather than in its derivative. In 3D, the established solution Kunc (1982); Giannozzi et al. (1991); Baroni et al. (2001); Gonze et al. (1994); Gonze and Lee (1997); Kern et al. (1999) is to construct a model of the long-ranged dipole-dipole interactions in reciprocal space, such that the corresponding contribution to the dynamical matrices can be excluded from the interpolation process. We apply the approach detailed in Ref. Gonze and Lee, 1997 to the 2D case, simply replacing the 3D screened Coulomb interaction with the 2D one (for the case of isotropic materials it is given in Eq. (3), for the anisotropic materials it is in the Appendix). In the long-wavelength limit, this amounts to excluding from the interpolation the following contribution to the dynamical matrix :

(4) |

where are Cartesian coordinates. In principle, the long-wavelength bare dipole-dipole interaction above is sufficient to treat the non-analyticity of the dynamical matrix giving rise to the power-law decay of IFCs. In practice, we also include screening at finite momenta via the use of the screened 2D Coulomb interaction. This helps convergence with respect to the phonon momentum grid from which the dispersions are interpolated. Indeed, screening brings some analytical but potentially sharp variations to the LO dispersion. Fig. 2 shows that this interpolation scheme is successful in reproducing the correct long-wavelength behavior. To highlight the importance of adapting DFT and DFPT to 2D materials, we also show how standard 3D periodic boundary conditions calculations fail in describing the long-wavelength behavior of polar-optical phonons in 2D.

For a comprehensive understanding of 2D systems we also study and propose a simple model of LO-TO splitting in multilayers. We take here h-BN as an example. Van-der-Waals interactions are neglected. These might affect the absolute value of the LO and TO frequencies by influencing interlayer distances but they would not change the LO-TO splitting itself since it is an electrostatic effect depending on in-plane Born effective charges and screening. We assume that, aside from thicker slabs, the dielectric properties of the multilayers remain unchanged. The underlying assumption is that the perturbing field and the material’s response to it are uniform over the multilayer. Furthermore, we focus on the the highest LO branch, with an in-phase LO mode in each monolayer, noted here as “hLO”. We find the frequency of this mode to be

(5) |

where the index “hTO” designates the highest TO mode (see Appendix for details). Both the strength of the splitting and the screening length are multiplied by a factor with respect to the monolayer. We compare the model of Eq. (5) with our DFPT calculations in Fig. 3, showing excellent agreement for the slope at and its evolution. Moving away from , the splitting is dictated by the screening of the 2D material and the model is in very good agreement considering the assumptions made about the dielectric properties of the multilayer. The model of Eq. (5) only treats the highest LO mode of the multilayer, for which the polarization densities are in-phase and the splitting effect from each layer adds up constructively. The lower LO modes, noted “out-of-phase LO” in Fig. 3, behave quite differently. In particular, they do not display a finite slope at zone center. Indeed, the electric fields generated by the polarization densities in different layers cancel each other in the long-wavelength limit, due to the atomic displacements being out-of-phase. The same argument explains the absence of splitting at zone center for the lower TO mode in bulk h-BN with AB stacking.

The expression in Eq. (5) offers further insights
on the 2D/3D transition.
As the number of layers increases, the slope of the LO mode increases
and the range of momenta over which the dispersion is linear decreases.
For increasing the slope of the highest LO mode becomes very large,
but the range of momenta for which the frequency increases closes up around
^{2}^{2}2Note that if one could probe phonons closer and closer
to as , the slope (or group velocity) would approach
and surpass the speed of light, and the non-relativistic electrostatic framework
used here would break down.
. Eventually, we arrive to a situation where all physically
relevant momenta are such that .
In this bulk limit (), the difference in the squared frequencies
approaches a constant value, giving rise to a finite LO-TO splitting close to
,
which can be estimated to be by considering
in Eq. (5)^{3}^{3}3More rigorously, one should rather
consider the limit
and , keeping the product of the two constant.
In BN, and using Table 1, we obtain
eV. This value is within
of the bulk splitting obtain in 3D DFPT (see Fig. 3),
and in excellent agreement with experimental results Rokuta et al. (1997); Geick et al. (1966)
in bulk h-BN. Such agreement points to a relatively easy way to estimate
this quantity from experiments or bulk calculations.

Last, in order to assess the role of environmental screening on LO-TO phonons we consider van-der-Waals heterostructures. These considerations are also relevant for transport properties of heterostructures, as the LO dispersion slope is representative of the strength of the long-ranged polarization fields that can scatter electrons remotely in all layers. In practice, it is essential to consider that the slope of the LO dispersion at the zone center is divided by the effective dielectric constant of the environment . Isolated monolayers are not always easily fabricated Rokuta et al. (1997); Serrano et al. (2007), and future devices will rely on Van-der-Waals heterostructures, in which the polar 2D material is surrounded by a variety of other 2D layers. We study the effect of the dielectric environment by simulating h-BN on top of monolayer and bilayer graphene (see Fig. 4). The calculations are performed with an electronic smearing equivalent to room temperature. This smearing is small enough to consider monolayer graphene as neutral, with a constant static screening function Kozinsky and Marzari (2006); Shung (1986); Gorbar et al. (2002); Ando (2006); Wunsch et al. (2006); Barlas et al. (2007); Wang and Chakraborty (2007); Hwang and Das Sarma (2007), namely Kozinsky and Marzari (2006); Sohier et al. (2015). At small momenta, the graphene layer then behaves like an insulating bulk dielectric and the slope of the LO dispersion is divided by . Bilayer graphene has a larger density of states at the intersection of the valence and conduction bands. In this case the smearing is enough for the bilayer graphene to exhibit a metallic behaviour in the long wavelength limit. The slope of the LO dispersion vanishes, as the polarization field from h-BN is completely screened by the electrons of bilayer graphene. Finally, note that in the general case of heterostructures containing other polar materials, the LO-TO splitting will be the result of a complex interplay between the polarization densities and screening properties of the various layers.

In conclusion, we use our implementation of the 2D Coulomb cutoff within density-functional perturbation theory to study the long wavelength limit of polar-optical phonons in a two-dimensional framework and complement this with physical, analytical models of the interactions and their screening. At the zone center, although the splitting of LO and TO phonon modes vanishes, a discontinuity appears in the slope of the LO phonon dispersion, and we provide a model to evaluate this slope in various situations. For isolated 2D materials, the slope can be estimated directly from the Born effective charges. For a multilayer, we find that the slope of the highest LO mode is proportional to the number of layers. In the general case, the slope also depends on the dielectric environment of the 2D material. In the experimentally relevant case of h-BN/graphene heterostructures, the slope is reduced, and can even vanish when screened by the metallic behaviour of the electrons. Last, screening from the electrons of the 2D material occurs only for phonon wavelengths smaller than an effective screening length. A wavelength-independent splitting similar to bulk 3D materials is then recovered.

Acknowledgements – This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 696656 â GrapheneCore1. T.S., M.G., and N.M. acknowledge financial support from the Swiss National Science Foundation (SNSF – project number 200021-143636). M.C. acknowledges support from the Agence Nationale de la Recherche under the reference no ANR-13-IS10-0003-01. Computer facilities were provided by the PRACE project and CINES, IDRIS and CEA TGCC (Grant EDARI No. 2016091202).

## Appendix A Computational details

Phonon calculations are performed within density functional perturbation
theory using a modified version of the PWSCF and PHONON codes of the Quantum ESPRESSO distribution. In particular, the modification includes the implementation
of the 2D Coulomb cutoff for the calculation of total energy, forces, bands
and the linear response of the system to a phonon perturbation.
We use pseudopotentials from the Standard Solid-State Pseudopotentials (SSSP) library ^{4}^{4}4see http://www.materialscloud.org/sssp (accuracy version),
with the exception of BN, for which we use
ultrasoft pseudopotentials within the local density approximation.
We use electron-momentum grids, except when graphene is involved in which case the grid is .
Calculations were performed using the AiiDA materials informatics platform Pizzi et al. (2016).

## Appendix B LO-TO splitting in the non-isotropic case

The expression for the polarization density induced by a given phonon mode in a single layer is:

(6) |

where , are respectively in plane and out-of-plane space variables, and is the out-of-plane profile of the monolayer, homogeneous to an inverse distance and normalized to unity (). Like in the main text, we use the limit of the atomic displacement pattern . The polarization density induces a potential with the same periodicity, given by Poisson equation

(7) |

where is a dielectric tensor depending on . Although we preserve generality in the notation of phonon mode , note that a field is effectively generated only if the divergence of the polarization is non-zero (which is true for optical modes with a longitudinal component). To solve Eq. 7, it is convenient to work in reciprocal space. In the context of 2D materials, since there is a periodicity only in the plane, we work with in-plane Fourier Transforms. However, the third dimension must be treated with care to account for screening properly. In the following, the subtleties of dimensionality and the treatment of the out-of-plane direction will be hidden in the screening. The in-plane Fourier Transform of is defined as follows. In a 2D-periodic framework, we integrate over the out-of-plane variable . The potential induced by a slab of polarization density extends much further than its source in the out-of-plane direction, as it decays as . This calls for a slightly different definition of the in-plane Fourier Transform in the 2D framework. Since we are only interested in the value of the potential within the material, we define the in-plane Fourier Transform of the potential as in 2D. Note that is normalized to unity. In the 3D-periodic framework, we work with the usual average over the unit cell and , where is the size of the unit cell in the out-of-plane direction.

Given those definitions, the Poisson equation is solved by:

(8) |

where is the screened Coulomb interaction. Assuming in-plane isotropy for the dielectric properties of the material, we can write it as in the main text:

(9) |

with , and is the in-plane dielectric constant of the bulk. See the appendix of Ref. Sohier et al., 2016 for more details. In the anisotropic case, and become tensors and we use the following model in the long-wavelength limit:

(10) |

Each component of can be obtained separately with the same method as in the isotropic case Sohier et al. (2016).

The associated electric field is:

(11) |

The corresponding force on atoms is:

(12) | ||||

(13) |

This bring the following additional term to the dynamical matrix:

(14) | ||||

(15) |

Where the index indicates that this is only the contribution related to Born-effective charges. After selecting the eigenvalue by multiplying left and right by (the limit of) a generic eigenvector , we obtain that the frequency is increased by:

(16) | ||||

(17) | ||||

(18) |

with or depending on the dimensionality. Here we have made no assumption on the dependency of phonon eigenvectors or Born effective charges on the direction of the phonon momentum. Using the LO phonon eigenvector in a material with hexagonal in-plane symmetry like in the main text, we have . In more general cases of 2D materials with lower symmetry, it is not always possible to define a pair of LO/TO modes that belong to the same irreducible representation at . In this case we would have

## Appendix C Evolution with the number of layers

To compute the splitting of the highest LO branch (hLO) in multilayers, we first assume that the effective charges are unchanged. We then notice that, with respect to the single layer, the squared term on the right-hand side of Eq. 18 is multiplied by a factor , where the numerator comes from the sum over the atoms while the denominator comes from the normalization of the phonon eigenvectors.

For the screening, we have the same Coulomb interaction, but with and because is proportional to the thickness of the 2D materialSohier et al. (2016). We thus have:

(19) | ||||

Where the squared term on the right-hand side is summed only on the atoms of one layer. Thus, it has the same value as in the monolayer case and all the consequences of having several layers is explicitly conveyed by the presence of the factors.

## References

- Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013), arXiv:1307.6718 .
- Fiori et al. (2014) G. Fiori, F. Bonaccorso, G. Iannaccone, T. Palacios, D. Neumaier, A. Seabaugh, S. K. Banerjee, and L. Colombo, Nature Nanotechnology 9, 768 (2014).
- Jariwala et al. (2014) D. Jariwala, V. K. Sangwan, L. J. Lauhon, T. J. Marks, and M. C. Hersam, ACS Nano 8, 1102 (2014), arXiv:1402.0047 .
- Wang et al. (2012) Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Nature Nanotechnology 7, 699 (2012).
- Xia et al. (2014) F. Xia, H. Wang, D. Xiao, M. Dubey, and A. Ramasubramaniam, Nature Photonics 8, 899 (2014).
- Verdi and Giustino (2015) C. Verdi and F. Giustino, Physical Review Letters 115, 176401 (2015).
- Sjakste et al. (2015) J. Sjakste, N. Vast, M. Calandra, and F. Mauri, Physical Review B - Condensed Matter and Materials Physics 92, 054307 (2015), arXiv:1508.06172 .
- Sohier et al. (2016) T. Sohier, M. Calandra, and F. Mauri, Physical Review B 94, 085415 (2016).
- Kunc (1982) K. Kunc, Physical Review Letters 48, 406 (1982).
- Giannozzi et al. (1991) P. Giannozzi, S. De Gironcoli, P. Pavone, S. Baroni, S. D. Gironcoli, P. Pavone, S. Baroni, S. De Gironcoli, P. Pavone, and S. Baroni, Physical Review B 43, 7231 (1991).
- Baroni et al. (2001) S. Baroni, S. De Gironcoli, A. Dal Corso, and P. Giannozzi, Reviews of Modern Physics 73, 515 (2001), arXiv:0012092v1 [arXiv:cond-mat] .
- Gonze et al. (1994) X. Gonze, J.-C. Charlier, D. C. Allan, and M. P. Teter, Physical Review B 50, 13035 (1994).
- Gonze and Lee (1997) X. Gonze and C. Lee, Physical Review B 55, 10355 (1997).
- Kern et al. (1999) G. Kern, G. Kresse, and J. Hafner, Physical Review B 59, 8551 (1999).
- Mele and Král (2002) E. J. Mele and P. Král, Physical Review Letters 88, 056803 (2002).
- Sánchez-Portal and Hernández (2002) D. Sánchez-Portal and E. Hernández, Physical Review B 66, 235415 (2002).
- Michel and Verberck (2009) K. H. Michel and B. Verberck, Physical Review B 80, 224301 (2009).
- Michel and Verberck (2011) K. H. Michel and B. Verberck, Physical Review B 83, 115328 (2011).
- Zólyomi et al. (2014) V. Zólyomi, N. D. Drummond, and V. I. Fal’ko, Physical Review B 89, 205416 (2014).
- Wirtz et al. (2003) L. Wirtz, A. Rubio, R. A. de la Concha, and A. Loiseau, Physical Review B 68, 045425 (2003).
- Sohier (2015) T. Sohier, Electrons and phonons in graphene : electron-phonon coupling , screening and transport in the field effect setup, Phd thesis, Universite Pierre et Marie Curie, Paris VI (2015).
- Ismail-Beigi (2006) S. Ismail-Beigi, Physical Review B 73, 233103 (2006).
- Rozzi et al. (2006) C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Physical Review B 73, 205119 (2006).
- Giannozzi et al. (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009), arXiv:0906.2569 .
- (25) See supporting information for the definition of the Fourier transform depending on the dimensionality of the system.
- Keldysh (1979) L. V. Keldysh, Journal of Experimental and Theoretical Physics Letters 29, 658 (1979).
- Wehling et al. (2011) T. O. Wehling, E. Sasioglu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Physical Review Letters 106, 236805 (2011), arXiv:1101.4007 .
- Cudazzo et al. (2011) P. Cudazzo, I. V. Tokatly, and A. Rubio, Physical Review B 84, 085406 (2011), arXiv:1104.3346 .
- Berkelbach et al. (2013) T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Physical Review B 88, 045318 (2013), arXiv:1305.4972 .
- Steinhoff et al. (2014) A. Steinhoff, M. Rösner, F. Jahnke, T. O. Wehling, and C. Gies, Nano Letters 14, 3743 (2014).
- Kozinsky and Marzari (2006) B. Kozinsky and N. Marzari, Physical Review Letters 96, 166801 (2006).
- Sohier et al. (2015) T. Sohier, M. Calandra, and F. Mauri, Physical Review B 91, 165428 (2015), arXiv:1503.02530v1 .
- Piscanec et al. (2004) S. Piscanec, M. Lazzeri, F. Mauri, A. C. Ferrari, and J. Robertson, Physical Review Letters 93, 185503 (2004).
- (34) Note that if one could probe phonons closer and closer to as , the slope (or group velocity) would approach and surpass the speed of light, and the non-relativistic electrostatic framework used here would break down.
- (35) More rigorously, one should rather consider the limit and , keeping the product of the two constant.
- Rokuta et al. (1997) E. Rokuta, Y. Hasegawa, K. Suzuki, Y. Gamou, C. Oshima, and A. Nagashima, Physical Review Letters 79, 4609 (1997).
- Geick et al. (1966) R. Geick, C. H. Perry, and G. Rupprecht, Physical Review 146, 543 (1966).
- Serrano et al. (2007) J. Serrano, A. Bosak, R. Arenal, M. Krisch, K. Watanabe, T. Taniguchi, H. Kanda, A. Rubio, and L. Wirtz, Physical Review Letters 98, 095503 (2007).
- Shung (1986) K. Shung, Physical Review B 34, 979 (1986).
- Gorbar et al. (2002) E. V. Gorbar, V. P. Gusynin, V. A. Miransky, and I. A. Shovkovy, Physical Review B 66, 045108 (2002).
- Ando (2006) T. Ando, Journal of the Physical Society of Japan 75, 074716 (2006).
- Wunsch et al. (2006) B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006).
- Barlas et al. (2007) Y. Barlas, T. Pereg-Barnea, M. Polini, R. Asgari, and A. H. MacDonald, Physical Review Letters 98, 236601 (2007).
- Wang and Chakraborty (2007) X.-F. Wang and T. Chakraborty, Physical Review B 75, 033408 (2007).
- Hwang and Das Sarma (2007) E. H. Hwang and S. Das Sarma, Physical Review B 75, 205418 (2007).
- (46) See http://www.materialscloud.org/sssp.
- Pizzi et al. (2016) G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, Computational Materials Science 111, 218 (2016).