Piezoelectricity and valley chern number in inhomogeneous hexagonal 2D crystals

Abstract

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 k·p model Hamiltonians. An analytical expression for the piezoelectric constant is obtained in terms of topological quantities, such as the 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₂. We further generalize 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 ≳10¹¹cm⁻² can be induced by realistic inhomogeneous strains, ϵ ≈ 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 third-rank piezoelectric tensor $\gamma _{ijk} \equiv \partial P_i/\partial \varepsilon _{jk}|_{{\mathbf{\varepsilon }} \to {\mathbf{0}}}$. The so-called modern theory of polarization exploits the properties of the Berry connection (BC) of the electronic wave-functions to quantify the change of polarization in an extended system.^3,4,5 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^6,7, as well as analytical approaches.^8,9

Inversion symmetry is broken in a large number of two-dimensional (2D) materials.¹⁰ This, together with their exceptional breaking strength and flexibility,¹¹ make them perfect platforms for strain engineering and, in particular, for piezoelectric applications.¹² Indeed, the isolation of monolayer and few-layer crystals of transition metal dichalcogenides (TMDs) or hexagonal Boron Nitride (h-BN) provides materials with symmetry properties that are different from their bulk counterparts. Bulk TMDs with the common formula MX₂ (M = Mo,W and X = S,Se) consist of stacked layers of MX₂ monolayers bonded by van der Waals forces, and have a center of symmetry located between the layers. Therefore, bulk TMDs are not piezoelectric. However, isolation of a monolayer of MX₂ from the bulk crystal removes the center of symmetry, leading to piezoelectricity, as reported experimentally.^13,14 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 $\sigma _{\mathrm{v}}:x \to - x$, and three-fold C₃ rotational symmetries, with the ${\hat{\boldsymbol x}}$ 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.^15,16,17

Realistic samples are often subject to non-uniform strain. This is particularly common in 2D crystals that are strained by controlled corrugation,¹⁸ deposition on substrates with nanodomes¹⁹ or nanopillars,^20,21 or because of the formation of bubbles due to trapped substances between the 2D crystal and the substrate.²² If the crystal is non-centrosymmetric, non-uniform strain can generate local charge density given by $\rho \left( {\boldsymbol{r}} \right) = - \nabla \cdot {\boldsymbol{P}}\left( {\boldsymbol{r}} \right)$, 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. Our theory relates the piezoelectric constant to the valley Chern number, and we show that piezoelectricity measurements can be used to access the valley degree of freedom in systems with a large gap like h-BN, for which non-local transport experiments are not viable. Explicit calculations for TMDs (MoS₂, MoSe₂, WS₂, and WSe₂) and h-BN are reported. Good agreement is found in comparison to existing ab initio calculations and experimental measurements of piezoelectric constant for MoS₂. We further studied the strain-induced polarization in undoped samples subject to non-uniform strain, like Gaussian and triangular bumps, bubbles, etc., and found that inhomogeneous deformations can induce large charge densities. TMDs are being extensively studied as platforms where the gap can be locally manipulated by strain,^19,23,24 and where strain gives rise to optical single-photon sources^20,21,25 (quantum emitters). Our theory can be used to determine the charge densities induced in these systems, as function of strain and device size.

Results

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 = \mathop {\sum}\nolimits_{jk} \gamma _{ijk}\varepsilon _{jk}$, where γ _ijk and ε _jk are the piezoelectric and the strain tensors, respectively. The quantity γ _ijk must respect the symmetries of the lattice, implying $\mathop {\sum}\nolimits_{i{\prime}j{\prime}k{\prime}} R_{ii{\prime}}^\dagger \gamma _{i{\prime}j{\prime}k{\prime}}R_{j{\prime}j}R_{k{\prime}k} = \gamma _{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₃ symmetries, we find that $\gamma _{xxx} = \gamma _{xyy} = \gamma _{yxy} = \gamma _{yyx} = 0$, and $\gamma _{yyy} = - \gamma _{yxx} = - \gamma _{xyx} = - \gamma _{xxy}$. The above symmetry properties allow us to write the piezoelectric-induced polarization as (See Supplementary Material for details).

$${\boldsymbol{P}} = \gamma _{yyy}{\cal A} \times \hat z,$$

(1)

where ${\cal A} = \left( {\varepsilon _{xx} - \varepsilon _{yy}} \right)\hat x - 2\varepsilon _{xy}{\hat{\mathbf y}}$. 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.²⁶ Notice also that, according to Eq. (1), the charge polarization is always perpendicular to ${\cal A}$, a result that has been reported in ref.⁹ 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 the electronic polarization of an insulator, P, can be calculated from the geometrical properties of the Bloch wave-functions,

$${\boldsymbol{P}}\left( {\cal A} \right) = - e\mathop {\sum}\limits_{\tau s} {\int}_{{\mathrm{B}}Z} \frac{{d{\boldsymbol{k}}}}{{\left( {2\pi } \right)^2}}{\boldsymbol{a}}_{\tau s}^{({\mathrm{v}})}\left( {{\boldsymbol{k}},\tau {\cal A}} \right),$$

(2)

where s = ± and τ = ± account for the spin and valley degrees of freedom, respectively. The valence-band BC reads ${\mathbf{a}}_{\tau s}^{({\mathrm{v}})}\left( {{\boldsymbol{k}},\tau {\cal A}} \right) = i\left\langle {u_{\tau s}^{({\mathrm{v}})}\left( {{\boldsymbol{k}},\tau {\cal A}} \right)\left| {\nabla _{\mathbf{k}}} \right|u_{\tau s}^{({\mathrm{v}})}\left( {{\boldsymbol{k}},\tau {\cal A}} \right)} \right\rangle$, where $\left| {u_{\tau s}^{({\mathrm{v}})}\left( {{\boldsymbol{k}},\tau {\cal A}} \right)} \right\rangle$ is the eigenvector of the system under strain. The piezoelectric constant is obtained as (See Supplementary Material for details).

$$\gamma = \frac{e}{2}\mathop {\sum}\limits_{\tau ,s} {\int}_{{\mathrm{B}}Z} \frac{{d{\boldsymbol{k}}}}{{\left( {2\pi } \right)^2}}\mathop {{lim}}\limits_{{\cal A} \to {\mathbf{0}}} \left[ {\frac{{\partial a_{\tau s,y}^{({\mathrm{v}})}\left( {{\boldsymbol{k}},\tau {\cal A}} \right)}}{{\partial A_x}} - \left( {x \leftrightarrow y} \right)} \right].$$

(3)

In the linear-response regime with respect to the strain gauge field, we can formally write the electronic Hamiltonian as $H_{\tau s}\left( {k,\tau {\cal A}} \right) \approx H_{\tau s}\left( k \right){\mathrm{ + }}\tau \mathop {\sum}\nolimits_\alpha A_\alpha \partial H_{\tau s}\left( {k,{\cal A}} \right)/\left. {\partial A_\alpha } \right|_{A \to 0} + O\left( {A^2} \right)$, where $H_{\tau s}\left( k \right)$ 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

$$\gamma = \frac{e}{{4\pi a_0}}\mathop {\sum}\limits_{\tau s} \tau \tilde C_{\tau s},$$

(4)

where a₀ is an effective lattice constant, $\tilde C_{\tau s} = {\int}_{{\mathrm{B}}Z} d{\boldsymbol{k}}\,\tilde \Omega _{\tau s}({\boldsymbol{k}})/2\pi$ has the form of the usual Chern number, and $\tilde \Omega _{\tau s}({\boldsymbol{k}})$ is formally similar to the Berry curvature (See Supplementary Material for details)

$$\begin{array}{l}\tilde \Omega _{\tau s}\left( {\boldsymbol{k}} \right) = \frac{i}{{\left[ {d_{\tau s}^{({\mathrm{c}}v)}\left( {\boldsymbol{k}} \right)} \right]^2}}\left\{ {\left\langle {u_{\tau s}^{\left( {\mathrm{v}} \right)}\left( {\boldsymbol{k}} \right)} \right|\tilde v_{\tau s,x}\left( {\boldsymbol{k}} \right)\left| {u_{\tau s}^{\left( {\mathrm{c}} \right)}\left( {\boldsymbol{k}} \right)} \right\rangle } \right.\cr \quad \quad \quad \, \times \left\langle {u_{\tau s}^{\left( {\mathrm{c}} \right)}\left( {\boldsymbol{k}} \right)} \right|v_{\tau s,y}\left( {\boldsymbol{k}} \right)\left| {u_{\tau s}^{\left( {\mathrm{v}} \right)}\left( {\boldsymbol{k}} \right)} \right\rangle - \left. {\left( {x \leftrightarrow y} \right)} \right\}.\end{array}$$

(5)

Here, $d_{\tau s}^{\left( {{\mathrm{cv}}} \right)}({\boldsymbol{k}}) = E_{\tau s}^{({\mathrm{c}})}({\boldsymbol{k}}) - E_{\tau s}^{({\mathrm{v}})}({\boldsymbol{k}})$ and $E_{\tau s}^{({\mathrm{c}}/v)}({\boldsymbol{k}})$ are the energy dispersion of the conduction/valence band. Notice that $v_{\tau s,\alpha }({\boldsymbol{k}}) = \partial H_{\tau s}({\boldsymbol{k}})/\partial k_\alpha$ is the standard band velocity, and the term $\tilde v_{\tau s,\alpha }({\boldsymbol{k}}) = a_0\partial H_{\tau s}({\boldsymbol{k}},{\cal A})/\partial A_\alpha |_{A \to 0}$ can be termed as “fictitious velocity”. For a generic two-band model for each (spin,valley) pair we can write ${\cal H}_{\tau s}({\boldsymbol{k}}) = \varepsilon _{\tau s}({\boldsymbol{k}})1 + {\boldsymbol{h}}_{\tau s}({\boldsymbol{k}}) \cdot {\mathbf{\sigma }}$, where ${\boldsymbol{h}}_{\tau s} = \left( {h_{\tau s,x},h_{\tau s,y},h_{\tau s,z}} \right)$ and ${\boldsymbol{\sigma }} = \left( {\sigma _x,\sigma _y,\sigma _z} \right)$ are Pauli matrices. The two-band model Hamiltonian of the strained crystal can be expressed as

$${\cal H}_{\tau s}({\boldsymbol{k}},{\cal A}) = \varepsilon _{\tau s}\left( {{\boldsymbol{k}} + \frac{{\eta _0{\cal A}}}{{2a_0}}} \right)1 + \mathop {\sum}\limits_i h_{\tau s,i}\left( {{\boldsymbol{k}} + \frac{{\eta _i{\cal A}}}{{2a_0}}} \right)\sigma _i,$$

(6)

where i = x, y, z and (η₀, η_i) are dimensionless Grüneisen parameters. Particularly, Eq. (6) has been obtained explicitly for monolayer h-BN and TMDs. After Eq. (6), we can evaluate the velocity as ${\boldsymbol{v}}_{\tau s}({\boldsymbol{k}}) = \nabla H_{\tau s}({\boldsymbol{k}})$, and the fictitious velocity as $\tilde v_{\tau s}({\boldsymbol{k}}) = \left\{ {\eta _0\nabla \varepsilon ({\boldsymbol{k}})1 + \mathop {\sum}\nolimits_i \eta _i\nabla h_{\tau s,i}({\boldsymbol{k}})\sigma _i} \right\}/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, $\Omega _{\tau s}({\boldsymbol{k}}) = 2\tilde \Omega _{\tau s}({\boldsymbol{k}})/\eta$ coincides with the conventional Berry curvature and $C_{\tau s} = 2\tilde C_{\tau s}/\eta$ 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 ${\boldsymbol{h}}_\tau ({\boldsymbol{k}}) = \left( {\hbar v\tau k_x,\hbar vk_y,\Delta } \right)$, where k = (k _x ,k _y ) is the relative momentum with respect to the K-point of the BZ, $\hbar v = 3t_0a_0/2$, where t₀~2.3 eV, Δ~5.97 eV, and $a = \sqrt 3 a_0 = 2.5{\mathrm{{\AA}}}$ are the nearest neighbor hopping, band gap, and lattice constant, respectively.^27,28 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²⁹

$$\begin{array}{l}\varepsilon _{\tau s}({\boldsymbol{k}}) = \frac{{\Delta _0 + \lambda _0\tau s}}{2} + \frac{{\hbar ^2|{\boldsymbol{k}}|^2}}{{4m_0}}\alpha ,\cr {\boldsymbol{h}}_{\tau s}({\boldsymbol{k}}) = \left( {t_0a_0\tau k_x,t_0a_0k_y,\frac{{\Delta + \lambda \tau s}}{2} + \frac{{\hbar ^2|{\boldsymbol{k}}|^2}}{{4m_0}}\beta } \right) .\end{array}$$

(7)

Here, m₀, the free electron mass, a₀, t₀, ∆₀, ∆, λ₀, λ, α, 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 1. We use Harrison’s prescription³⁰ for Slater–Koster parameters in our tight-binding modeling of strained TMDs.²⁹ Briefly, in this approach, the hopping between two lattice sites depends on the inter-site distance, r, through $\propto 1/r^{\ell _\mu + \ell _\nu + 1}$, where $\ell _{\mu (\nu )}$ stands for the angular momentum of two orbitals localized on two sites. In strained h-BN, we only have one independent particle-strain coupling, η _x  = η _y  = η~3.3⁹ leading to the simple relation $\tilde v_\tau = \eta v_\tau /2$. As a consequence, $\tilde C_\tau$ is proportional to the usual Chern number in the massive Dirac model, i.e., $C_\tau = \tau {\mathrm{sign}}[\Delta ]/2$, and we find

$$\gamma ^{{\mathrm{h}} - {\mathrm{BN}}} = \eta \frac{e}{{4\pi a_0}}C_{{\mathrm{v}}alley},$$

(8)

where the valley Chern number is defined by $C_{{\mathrm{valley}}} = \mathop {\sum}\nolimits_\tau \tau C_\tau = {\mathrm{sign}}[\Delta ]$.^31,32 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 non-local transport measurement in multi-terminal devices.^33,34,35 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³³ and bilayer graphene.^34,35 We notice that applying time-dependent strain can induce a synthetic valley-dependent electric field, which can derive charge current in topological systems, such as gapped graphene.^36,37 The numerical value of $\gamma ^{{\mathrm{h}} - {\mathrm{BN}}}$ obtained from our theory is given in Table 2, showing good agreement with ab initio calculations.

Table 1 k·p parameters of TMDs extracted from the low-energy projection of a tight-binding model.²⁹ The lattice parameters are taken from ref.⁶

Table 2 The numerical value of piezoelectric constant $\gamma [{\mathrm{1}}0^{ - 10}C/m]$ obtained here (k·p method), and the previously reported DFT (clamped-ion) and experimental results

The case of strained TMDs is more complex, since η _x  = η _y while they differ from η₀ 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 $\tilde C_{\tau 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 Supplementary Material for details)

$$\gamma ^{{\mathrm{TMDs}}} = (\eta _x + \eta _z)\frac{e}{{8\pi a_0}}\mathop {\sum}\limits_{s = \pm } \frac{{D_s(\eta _x,\eta _z) + 2\beta \Lambda _sC_s}}{{1 + 2\beta \Lambda _s}},$$

(9)

where $\Lambda _s = \hbar ^2(\Delta + s\lambda )/(4m_0t_0^2a_0^2)$,

$$C_s = \frac{{{\mathrm{sign}}[\Delta + s\lambda ] - {\mathrm{sign}}[\beta ]}}{2},$$

(10)

and

$$D_s(\eta _x,\eta _z) = \frac{{\eta _x{\mathrm{sign}}[\Delta + s\lambda ] - \eta _z{\mathrm{sign}}[\beta ]}}{{\eta _x + \eta _z}}.$$

(11)

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.^17,38 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

$$\gamma ^{{\mathrm{TMDs}}} = \frac{e}{{4\pi a_0}}\left( {\begin{array}{*{20}{c}} {\eta _x + \eta _z} & {{\mathrm{\beta < 0}}} \cr {\eta _x} & {{\mathrm{\beta = 0}}} \cr {\mathop {\sum}\limits_{s = \pm } \frac{{\eta _x - \eta _z}}{{2 + 4\beta \Lambda _s}}} & {{\mathrm{\beta > 0}}} \end{array}} \right..$$

(12)

The values of γ^TMDs obtained from our k·p method are shown in Table 2. 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.

Discussion

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.^{18,19,20,21,22,39,40,41,42} Neglecting, for long-wavelength deformations, the flexoelectric (i.e., a term accounting for the polarization induced by the strain gradient⁴³), we can generalize to the inhomogeneous strain case the linear-response relation for the piezoelectric tensor: $P_i({\boldsymbol{r}}) \approx \mathop {\sum}\nolimits_{jk} \gamma _{ijk}\varepsilon _{jk}({\boldsymbol{r}})$. Consequently, the induced charge density, following Eq. (1), reads $\rho ({\boldsymbol{r}}) = en({\boldsymbol{r}}) = - \nabla \cdot {\boldsymbol{P}}({\boldsymbol{r}}) \approx - \gamma \hat z \cdot \left[ {\nabla \times A({\boldsymbol{r}})} \right]$. 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 $\ell$. The variation of the strain leads to $\nabla \times {\cal A}({\boldsymbol{r}})\sim \Delta \varepsilon /\ell$, so that $n({\boldsymbol{r}})\sim (\gamma /e)\Delta \varepsilon /\ell \sim \Delta \varepsilon /(a_0\ell )$. The materials considered here can withstand large strains. In MoS₂ or h-BN bubbles,²² the variations in the strain can be of order Δε ≈ 0.02 over scales of $\ell \sim 300$ nm, leading to $n\sim 10^{11}{\mathrm{cm}}^{ - 2}$. Higher strain gradients, with maximum strains $\Delta \varepsilon \sim 0.1$ over short lengths, $\ell \sim 10 - 15$ nm have been reported in graphene bubbles on metallic surfaces.³⁹ Similar configurations will induce particle densities $n\sim 10^{13}{\mathrm{cm}}^{ - 2}$. We notice that, for larger values of applied strain, second-order piezoelectric effects might be relevant for an accurate estimation of the induced particle density. This is the case in Zinc-Blende⁴⁴ and wurtzite⁴⁵ semiconductors.

In order to illustrate further, the charge induced by non-uniform strains, we discuss in detail the case of MoS₂ and h-BN bubbles described in ref.²² 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 ${\mathbf{P}}({\mathbf{r}}) = p(r)[{\hat{\mathbf r}}sin(3\theta ) + {\hat{\mathbf \theta }}cos(3\theta )]$ where (See Supplementary Material for details.):

$$p(r) = \gamma \left\{ {\frac{{u(r)}}{r} - \frac{{\partial u(r)}}{{\partial r}} - \frac{1}{2}\left( {\frac{{\partial h(r)}}{{\partial r}}} \right)^2} \right\} = p_0\tilde p\left( {\frac{r}{R}} \right) .$$

(13)

Here, $p_0 = \gamma (1 + \nu )h_{{\mathrm{max}}}^2/R^2$, with v the Poisson’s ratio, and $\tilde p(x)$ is a universal function which does not depend on the material. The induced charge density is

$$\rho ({\boldsymbol{r}}) = \rho _0\tilde \rho \left( {\frac{r}{R}} \right)sin(3\theta )$$

(14)

where ρ₀ = p₀/R and, as before, the function $\tilde \rho (x)$ is universal. This analysis is consistent with our previous estimates, as $\Delta \varepsilon \sim h_{{\mathrm{max}}}^2/R^2$, and $\ell \sim R$. The charge distribution is shown in Fig. 1. 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⁻³ 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⁻¹. For γ~1 and h_max/R~0.1, we find $\rho \sim 10^7/R$ in units of electron charge × cm⁻². For R ≈ 1 μm, we obtain $\rho _0 \approx 10^{11}e \times {\mathrm{cm}}^{ - 2}$(See Supplementary Material for details). Finally, a number of schemes have been proposed to study gauge fields in graphene.^46,47,48 If the graphene layer in devices, such as quantum emitters,^20,21,25 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 \approx 10^{11}{\mathrm{cm}}^{ - 2}$.

Fig. 1

a Radial and out-of-plane displacement profiles corresponding to a pure-bending deformation with $\tilde h(x) \equiv h(x)/h_{{\mathrm{m}}ax} = (1 - x^2)^2\Theta (1 - x)$ and the dimensionless radial displacement, $\tilde u(x) \equiv u(x)/u_0$ with $u_0 = h_{{\mathrm{max}}}^2/R$, obtained from elasticity theory (See Supplementary Material for details). The horizontal red dashed line is given as a guide to the eye, showing the zero displacement level. b Bubble shape with a pure-bending deformation profile with $h_{{\mathrm{max}}}\sim 0.11R$. c Polarization field, $\tilde p(x)[{\hat{\boldsymbol r}}sin(3\theta ) + {\hat{\boldsymbol \theta }}cos(3\theta )]$, streamlines on top of the induced charge density, $\tilde \rho (x)sin(3\theta )$, colormap. The horizontal (vertical) colorbar corresponds to the streamline (colormap) plot. The dashed circle indicates the bubble’s boundary at r = R

In summary, we have performed a systematic study of the piezoelectric response of non-centrosymmetric 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 valleytronics^49,50, and piezotronics.¹² The particular cases of h-BN and TMDs (MoS₂, WS₂, MoSe₂, and WSe₂) 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 tremendous computational cost. We demonstrate that piezoelectric effect in inhomogenous crystals leads to the appearance of significant charge densities in the sample.

Data Availability

All data that support the findings of this study are available on request from the corresponding authors.