Piezoelectricity and valley Chern number in inhomogeneous hexagonal 2D crystals

Conversion of mechanical forces to electric signal is possible in non-centrosymmetric materials due to linear piezoelectricity. The extraordinary mechanical properties of two-dimensional materials and their high crystallinity make them exceptional platforms to study and exploit the piezoelectric effect. Here, the piezoelectric response of non-centrosymmetric hexagonal two-dimensional crystals is studied using the modern theory of polarization and ${\bm k} \cdot {\bm p}$ model Hamiltonians. An analytical expression for the piezoelectric constant is obtained in terms of topological quantities such as the {\it valley Chern number}. The theory is applied to semiconducting transition metal dichalcogenides and hexagonal Boron Nitride. We find good agreement with available experimental measurements for MoS$_2$. We further generalise the theory to study the polarization of samples subjected to inhomogeneous strain (e.g.~nanobubbles). We obtain a simple expression in terms of the strain tensor, and show that charge densities $\gtrsim 10^{11} {\rm cm}^{-2}$ can be induced by realistic inhomogeneous strains, $\epsilon \approx 0.01 - 0.03$.

Introduction.-Piezoelectricity is a property of crystals with broken inversion symmetry, which allows conversion of mechanical to electric energy [1,2]. When subjected to an external strain field ε, piezoelectric crystals acquire a polarization P that is described by the thirdrank piezoelectric tensor γ ijk ≡ ∂P i /∂ε jk | ε→0 . The socalled modern theory of polarization exploits the properties of the Berry connection (BC) of the electronic wavefunctions to quantify the change of polarization in an extended system [3][4][5][6]. For crystalline insulators, the BC is obtained in terms of the Bloch orbitals, and the polarization can thus be calculated as an integral of the BC on whole Brillouin zone. This quantum mechanical description of the polarization has been used to calculate the piezoelectric constant of a number of crystals from ab initio [7,8] as well as analytical approaches [9,10].
Inversion symmetry is broken in a large number of two-dimensional (2D) materials [11]. This, together with their exceptional breaking strength and flexibility [12,13], make them perfect platforms for strain engineering and, in particular, for piezoelectric applications [14]. Indeed, the isolation of monolayer and few-layer crystals of transition metal dichalcogenides (TMDs) or hexagonal Boron Nitride (h-BN) [15] provides materials with symmetry properties that are different from their bulk counterparts. Bulk TMDs with the common formula M X 2 (M = Mo,W and X = S,Se) consist of stacked layers of M X 2 monolayers bonded by van der Waals forces, and have a centre of symmetry located between the layers. Therefore bulk TMDs are not piezoelectric. However, isolation of a monolayer of M X 2 from the bulk crystal removes the center of symmetry, leading to piezoelectricity, as reported experimentally [16,17]. Similarly, h-BN consists of a honeycomb lattice with different elements in the two sublattices of the unit cell, making this material piezoelectric as well. Both monolayer h-BN and TMDs are hexagonal crystals that belong to the D 3h point symmetry group which contains two main symmetry elements: mirror reflection σ v : x → −x, and threefold C 3 rotational symmetries, with thex axis along the zigzag direction. They present a direct band gap at the two inequivalent K and K points of the Brillouin zone (BZ), and their low-energy electronic excitations are well described by massive Dirac-like Hamiltonians [18][19][20].
Realistic samples are often subject to non-uniform strain. This is particularly common in 2D crystals which are strained by controlled corrugation [21,22], deposition on substrates with nanodomes [23] or nanopillars [24,25], or because of the formation of bubbles due to trapped substances between the 2D crystal and the substrate [26]. If the crystal is non-centrosymmetric, nonuniform strain can be a source of carrier doping, with charge density given by ρ(r) = −∇ · P (r) [27], due to local variations of the polarization. This issue will be one of the main focuses of this work.
Using a generic k · p model Hamiltonian, we derive analytical expressions for the piezoelectric coefficients of hexagonal 2D crystals. Explicit calculations for TMDs (MoS 2 , MoSe 2 , WS 2 and WSe 2 ) and h-BN are reported. Good agreement is found in comparison to existing ab initio calculations and experimental measurements of piezoelectric constant. We further study the strain-induced polarization in undoped samples subject to non-uniform strain, like Gaussian and triangular bumps, bubbles, etc., finding that inhomogeneous deformations can induce large charge densities. TMDs are being extensively studied as platforms where the gap can be locally manipulated by strain [23,[28][29][30], and where strain give rises to optical single-photon sources [24,25,31] (quantum emitters). Our theory can be used to determine the charge densities induced in these systems, as function of strain and device size.
General formulation.-Let us consider a 2D crystal subject to a uniform static strain field. In the linear response regime, the induced-polarization due to the piezoelectric effect is given by P i = jk γ ijk ε jk where γ ijk and ε jk are the piezoelectric and the strain tensors, respec-tively. The quantity γ ijk must respect the symmetries of the lattice, implying i j k R † ii γ i j k R j j R k k = γ ijk where R accounts for a point group symmetry element. For instance, for a 2D system with D 3h symmetry lying in the xy-plane, after considering σ v and C 3 symmetries, we find that γ xxx = γ xyy = γ yxy = γ yyx = 0, and γ yyy = −γ yxx = −γ xyx = −γ yxx . The above symmetry properties allow us to write the piezoelectric-induced polarization as (see Appendix A) where A = (ε xx − ε yy )x − 2ε xyŷ . It is worth noting that this expression is formally equivalent to the gauge field that describes the effect of strains on the electronic structure of graphene [32]. Notice also that, according to Eq. (1), the charge polarization is always perpendicular to A, a result that has been reported in Ref. [10]. Finally, Eq.
(1) also implies that the trace of the strain tensor, V = ε xx + ε yy , does not contribute to the polarization. From now on, we use γ to indicate γ yyy . According to the modern theory of polarization [3][4][5][6], the electronic polarization of an insulator, P , can be calculated from the geometrical properties of the Bloch wave-functions, where s = ± and τ = ± account for the spin and valley degrees of freedom, respectively. The valence-band BC reads a is the eigenvector of the system under strain. The piezoelectric constant is obtained as (see (3) In the linear response regime with respect to the strain gauge field, we can formally write the electronic Hamiltonian as H τ s (k, τ A) ≈ H τ s (k) + τ α A α ∂H τ s (k, A)/∂A α | A→0 +O(A 2 ), where H τ s (k) is the unstrained Hamiltonian. Moreover, the valence band eigenvector can be evaluated by using first-order perturbation theory. After some straightforward calculations, we obtain the following expression for the piezoelectric coefficient where a 0 is an effective lattice constant, C τ s = BZ dk Ω τ s (k)/2π has the form of the usual Chern number, and Ω τ s (k) is formally similar to the Berry curvature (see Appendix B) τ s (k), and E (c/v) τ s (k) are the energy dispersion of the conduction/valence band. Notice that v τ s,α (k) = ∂H τ s (k)/∂k α is the standard band velocity, and the term v τ s,α (k) = a 0 ∂H τ s (k, A)/∂A α | A→0 can be termed as "fictitious velocity".
For a generic two-band model for each (spin,valley) pair we can write H τ s (k) = τ s (k)1 + h τ s (k) · σ, where h τ s = (h τ s,x , h τ s,y , h τ s,z ) and σ = (σ x , σ y , σ z ) are Pauli matrices. The two-band model Hamiltonian of the strained crystal can be expressed as where i = x, y, z and (η 0 , η i ) are dimensionless parameters accounting for the strength of particle-strain coupling. Particularly, Eq. (6) has been obtained explicitly for graphene, monolayer h-BN and monolayer TMDs. After Eq. (6), we can evaluate the velocity as v τ s (k) = ∇H τ s (k), and the fictitious velocity as v τ s (k) = {η 0 ∇ (k)1 + i η i ∇h τ s,i (k)σ i } /2. Notice that, for the simplest graphene-like case, η i=0,x,y,z = η, the fictitious velocity is proportional to the velocity. In this case, therefore Ω τ s (k) = 2 Ω τ s (k)/η coincides with the conventional Berry curvature and C τ s = 2 C τ s /η is the Chern number.
Piezoelectric constant of h-BN and TMDs.-In the following, we apply the developed theory to calculate the piezoelectric constant of two paradigmatic families of 2D crystals with D 3h symmetry: h-BN and TMDs. The effective k · p Hamiltonian of h-BN in the "sublattice" space is given by τ (k) = 0 and h τ (k) = ( vτ k x , vk y , ∆), where k = (k x , k y ) is the relative momentum with respect to the K-point of the BZ, v = 3t 0 a 0 /2, where t 0 ∼ 2.3 eV, ∆ ∼ 5.97 eV, and a = √ 3a 0 = 2.5Å, are the nearest neighbor hopping, band gap and lattice constant, respectively [10,33,34]. The spin degree of freedom leads to a double degeneracy of the states and therefore we drop the subindex "s" in the h-BN Hamiltonian. On the other hand, the effective k · p model for monolayer TMDs (ignoring trigonal warping effects) in "band" space is [35] τ s (k) = Here, m 0 , the free electron mass, a 0 , t 0 , ∆ 0 , ∆, λ 0 , λ, α, and β are strain-independent parameters that can be obtained in terms of the Slater-Koster parameters of the original tight-binding Hamiltonian. Numerical values for the different monolayer TMDs considered in this work are given in Table I.
In strained h-BN, we only have one independent particle-strain coupling, η x = η y = η ∼ 3.3 [10] leading to the simple relation v τ = ηv τ /2. As a consequence, C τ I: k · p parameters of TMDs extracted from the lowenergy projection of a tight-binding model [35,36]. The lattice parameters are taken from Ref. [7].

MoS2
MoSe2 WS2 WSe2 where the valley Chern number is defined by . This result differs from the one reported in Ref. [10], which contains a high-energy cutoff that does not appear in our derivation. Therefore, measurements of the piezoelectric constant can be used as direct tools to analyze the valley Chern number. Topological valley currents have been recently detected through nonlocal transport measurement in multi-terminal devices [37][38][39]. Here, we propose that piezoelectricity measurements can be used to access the valley degree of freedom in systems like h-BN, whose large gap impedes non-local transport experiments like those performed in graphene superlattices [37] and bilayer graphene [38,39]. The numerical value of γ h-BN obtained from our theory is given in Table. II, showing good agreement with ab initio calculations.
The case of strained TMDs is more complex, since η x = η y while they differ from η 0 and η z . Contrary to the simpler h-BN case, the fictitious velocity in TMDs is not proportional to the velocity and consequently we cannot use the simplified relation with the usual Chern number. However, we still can evaluate C τ s explicitly from the TMD's Hamiltonian. After a straightforward calculation, we arrive at the following analytical expression for the piezoelectric constant of TMDs (see Appendix C) where Λ s = 2 (∆ + sλ)/(4m 0 t 2 0 a 2 0 ) , and Notice that C s is the usual K-valley (τ = +) Chern number of monolayer TMDs with spin s = ±. Intriguingly, depending on the relative sign of β and ∆ ± λ, we either have a topological (C s = ±1) or a trivial (C s = 0) phase in each valley [20,40]. This topological property is protected as long as inter-valley scattering is suppressed. Since ∆ ± λ > 0 for the case of interest here, one can simplify Eq. (9). We find The values of γ TMDs obtained from our k · p method are shown in Table II. Again, in spite of the simplicity of our model, the results that we find are in good agreement with existing ab-initio and experimental results, strengthening the validity of our approach and providing microscopic insight into piezoelectricity in 2D crystals. -2.9 ± 0.5 ---Effect of inhomogeneous strain.-In the following we consider the polarization induced in samples subjected to inhomogeneous strain. This is a highly relevant problem due to the large number of recent experiments in which 2D crystals are subjected to a non-uniform strain profile [21][22][23][24][25][26][41][42][43][44][45]. Neglecting, for long-wavelength deformations, the flexoelectric (i.e. a term accounting for the polarization induced by the strain gradient [46]), we can generalize to the inhomogeneous strain case the linear-response relation for the piezoelectric tensor: P i (r) ≈ jk γ ijk ε jk (r). Consequently the induced charge density, following Eq. (1), reads ρ(r) = en(r) = −∇ · P (r) ≈ −γẑ · [∇ × A(r)]. The dependence of ρ(r) on the strain tensor is the same as the dependence of the strain-induced pseudomagnetic field acting on electrons in graphene. Unlike this case, the induced charge density has the same sign in the two valleys.
We can obtain simple estimates of the charge density induced by a variation in strain, ∆ε, over a length . The variation of the strain leads to ∇ × A(r) ∼ ∆ε/ , so that n(r) ∼ (γ/e)∆ε/ ∼ ∆ε/(a 0 ). The materials considered here can withstand large strains. In MoS 2 or h-BN bubbles [26], the variations in the strain can be of order ∆ ≈ 0.02 over scales of ∼ 300 nm, leading to n ∼ 10 11 cm −2 . Higher strain gradients, with maximum strains ∆ ∼ 0.1 over short lengths, ∼ 10 − 15 nm have been reported in graphene bubbles on metallic surfaces [41]. Similar configurations will induce carrier densities n ∼ 10 13 cm −2 .
In order to illustrate further the charge induced by non-uniform strains, we discuss in detail the case of MoS 2 and h-BN bubbles described in [26]. The shape of these bubbles is determined by the competition between the elastic energy of the 2D material and the van der Waals attraction to the substrate. We consider bubbles with radial symmetry. The shape and internal strains are defined by the in plane and out of plane displacements, u(r), h(r). The form of these functions are universal, and determined by the ratio h max /R, where h max is the height of the bubble and R is its radius. The polarization vector for this case is given by P (r) = p(r)[r sin(3θ)+θ cos(3θ)] where (see Appendix D) Here, p 0 = γ(1 + ν)h 2 max /R 2 , with ν the Poisson's ratio, and p(x) is a universal function which does not depend on the material. The induced charge density is where ρ 0 = p 0 /R and, as before, the function ρ(x) is universal. This analysis is consistent with our previous estimates, as ∆ε ∼ h 2 max /R 2 , and ∼ R. The charge distribution is shown in Fig. 3. The charge density reflects the trigonal symmetry of the lattice, and, as a result, it vanishes at the apex of the bubble, r = 0. We assume that the piezoelectric layer slides and relaxes outside the region where it is detached from the substrate. As a result, the charge density decays as r −3 outside the bubble. Note that the aspect ratio, h max /R is independent of the size of the bubble, so that the size dependence of the induced charge density is R −1 . For γ ∼ 1 and h max /R ∼ 0.1, we find ρ ∼ 10 7 /R in units of electron charge × cm −2 . For R ≈ 1µm, we obtain ρ 0 ≈ 10 11 e × cm −2 . For more details, see Appendix D. Finally, a number of schemes have been proposed to study gauge fields in graphene [47][48][49]. If the graphene layer in devices, such as quantum emitters [24,25,31], was encapsulated in h-BN, the pseudomagnetic fields in graphene and the charge density induced in the h-BN layers are proportional. For example, a pseudomagnetic field of one Tesla in graphene implies n ≈ 10 11 cm −2 .
Conclusions.-In summary, we have performed a systematic study of the piezoelectric response of noncentrosymmetric hexagonal 2D crystals. Starting from a general k · p model Hamiltonian, we have obtained a closed analytical expression for the piezoelectric constant in terms of the valley Chern number, bridging valleytron-ics valleytronics [50,51] and piezotronics [14]. The particular cases of h-BN and TMDs (MoS 2 , WS 2 , MoSe 2 and WSe 2 ) have been studied. The validity of the theory has been proven by the good quantitative agreement found between the piezoelectric constant obtained from our method, and that calculated from ab initio approaches and experimental measurements.
We finally generalize the theory to study samples subjected to inhomogeneous strain, which is a case of great experimental interest and which cannot be studied with standard DFT methods due to the computational cost. We demonstrate that piezoelectric effect in inhomogenous crystals leads to the appearance of significant carrier densities in the sample.

Appendix A: Symmetry consideration
Within the linear response theory, the strain-induced polarization reads where P and ε are the polarization vector and strain tensor respectively, and γ ijk stands for the piezoelectric tensor element. For the rank-3 tensor γ, we can obtain the following symmetry relation [52] i j k where R stands for the matrix representation of a symmetry operator. Considering three-fold rotational, C 3 , and vertical mirror, Notice that Owing to the mirror symmetry, all tensor elements with an odd number of x Cartesian index are identically zero. Considering the symmetry relations given in Eq. (A3), we can rewrite Eq. (A1) as follows where According to this relation the charge polarization is always perpendicular to the vector A.

Appendix B: Piezoelectric constant
The electronic polarization is given by [3][4][5][6] where τ A is the fictitious gauge field at τ -valley and the superscript (v) stands for the valence band index. The piezoelectric constant reads We can use the following decomposition On the other hand, the first term in Eq. (B3) does not respect such anti-symmetric relation under x ↔ y exchange. Therefore, the integral of the first symmetric term in the above relation on whole BZ is zero based on the D 3h symmetry. In this regard, we reach the following relation which is identical to Eq. (3) of the main text.
Using the definition of Berry connection, we can obtain We linearize the strain-electron interaction part: Notice that a 0 is an effective lattice constant and the fictitious velocity is defined as follows Within first order perturbation theory, the Bloch eigen-function reads where the sum runs over band indices m. For a two-bands model, we have We use the following identity [53] u (c) τ s (k) where the velocity operator is v τ s,α (k) = ∂H τ s (k) ∂k α .
(B11) Therefore, we reach From this we obtain the following relation for the piezoelectric constant.

Appendix C: Piezoelectric constant of TMDs
The two-band model Hamiltonian of strained monolayer TMDs can be formally written as follows [35] H τ s (k, where i = x, y, z and (η 0 , η i ) are dimensionless parameters accounting for the strength of particle-strain coupling of each term. Notice that τ /s indicate of the valley/spin index. The eigenvalues of unstrained Hamiltonian, H τ s (k, 0), read where |h τ s (k)| = h τ s,x (k) 2 + h τ s,y (k) 2 + h τ s,z (k) 2 .
The corresponding eigenvectors are Note that +/− sign corresponds to the conduction(c)/valence(v) band. According to the definition of velocity and fictitious velocity, we have Notice that the off-diagonal strain tensor elements are zero due to symmetry, i.e. ε rθ = ε θr = 0 [54,55]. Using the polar coordinates, r = (r, θ), it is straightforward to find: The electronic polarization is given by Therefore, the induced charge density reads In order to calculate the polarization and induced charge density, we first need to evaluate the strain tensor. For this purpose, we define strain tensor per area through the stress tensor, Σ αβ , and it can be described by the Hook's law [54,55] Σ rr (r) = E 1 − ν 2 [ε rr (r) + νε θθ (r)] , where E and ν are the Young's modulus and Poisson's ratio, respectively. At this stage, we need to obtain Σ rr(θθ) .
We define the following Green's function equation Therefore, we have where the source function follows and the general solution is which contains two unknown parameters c 1 and c 2 which can be obtained by considering continuity of the strain tensor. It is easy to prove that where r > = max(r, r ) and r < = min(r, r ). By plugging the above Green's function in Eq. (D11), we find the following solution for the radial displacement Using the dimensionless variables given in Table III, we can write where From the knowledge of the radial displacement, we can calculate the strain tensor elements by using Eq. (D1). Eventually, the radial profile of polarization and induced charge density can be evaluated from the following relations and The rest of this section is devoted to three particular cases of experimental interest, namely bubbles with different shapes as Gaussian, parabolic-like [26] and pure-bending profiles. The parabolic-like bubble has been discussed in Ref. [26]. Here, we briefly explain the pure-bending case for which the bending contribution to the mechanical free energy, ∝ ∇ 2 h(r) 2 , dominates the elastic energy, ∝ |∇h(r)| 2 . In this approximation, the bubble profile is given is the flexural rigidity. Considering a uniform force (per area) profile as f z (r) = f Θ(R − r), with Θ(x) the Heaviside function, and implementing h(R) = ∂h(r)/∂r| r→R = 0 as the boundary condition, the out-of-plane displacement reads h(r) = h max h(x) with h(x) = Θ(1 − x) 1 − x 2 2 and h max = f R 4 /64D [54,55]. Straightforward calculations lead to explicit expressions for the displacement vector, { u(x), h(x)}, the strain tensor, { ε rr (x), ε θθ (x)}, the radial profile of electronic polarization, p(x), and the radial profile of induced charge density, ρ(x), for each of the three bubble profiles discussed in the text, as shown in Tables IV, V and VI. All of these mathematical functions are shown in Fig. 3. IV: Functional expressions for the displacement vector, strain tensor, radial profile of electronic polarization and radial profile of induced charge density, for a 2D crystal deformed following a Gaussian bubble profile.

Quantity
Functional expression . We set ν = 0 in these plots.