Bond-Orbital-Resolved Piezoelectricity in Sp2-Hybridized Monolayer Semiconductors

Sp2-hybridized monolayer semiconductors (e.g., planar group III-V and IV-IV binary compounds) with inversion symmetry breaking (ISB) display piezoelectricity governed by their σ- and π-bond electrons. Here, we studied their bond-orbital-resolved electronic piezoelectricity (i.e., the σ- and π-piezoelectricity). We formulated a tight-binding piezoelectric model to reveal the different variations of σ- and π-piezoelectricity with the ISB strength (Δ). As Δ varied from positive to negative, the former decreased continuously, but the latter increased piecewise and jumped at Δ=0 due to the criticality of the π-electrons’ ground-state geometry near this quantum phase-transition point. This led to a piezoelectricity predominated by the π-electrons for a small |Δ|. By constructing an analytical model, we clarified the microscopic mechanisms underlying the anomalous π-piezoelectricity and its subtle relations with the valley Hall effect. The validation of our models was justified by applying them to the typical sp2 monolayers including hexagonal silicon carbide, Boron-X (X = N, P, As, Ab), and a BN-doped graphene superlattice.


Introduction
Piezoelectricity, a linear electromechanical intercoupling in non-centrosymmetric dielectrics, is an intriguing subject in solid state physics and plays an important role in various technological applications. Piezoelectric polarization of a crystal contains ionic and electronic contributions. The former stems from the internal ionic displacement caused by strain, while the latter is described by the strain-induced geometric (or Berry) phase shift of the ground-state wavefunction [1][2][3]. Specifically, the electronic piezoelectric coefficient is determined by the Brillouin zone (BZ) integration of the piezoelectric Berry curvature (PBC) defined in the parameter space containing strain [4][5][6]. The PBC formally resembles the momentum space Berry curvature (MBC) [7], which is crucial for understanding many geometric phase phenomena such as the quantum Hall effect [8] and quantum phase transition (QPT) [9][10][11], etc. Moreover, encoded with the information of the ground state, electronic piezoelectricity may also serve in detecting the QPT characterized by a Berry phase jump [12][13][14]. It is therefore worth investigating the mechanisms of electronic piezoelectricity and its correlation with other geometric phase phenomena.
Technically, nanomaterials with excellent piezoelectricity are pursued for next-generation electronic devices. Recent breakthroughs in two-dimensional (2D) piezoelectricity offer new opportunities in this respect. Many 2D materials [15][16][17][18][19][20][21][22][23][24][25][26] such as 2D transition metal dichalcogenides [16], graphitic carbon nitride [18], and hexagonal boron nitride [19], have been theoretically predicted or experimentally confirmed to be piezoelectric. Among them, sp 2 -hybridized monolayers; e.g., planar group III-V [19][20][21][22][23][24] and IV-IV [24][25][26] binary compounds, have attracted increasing attention due to their bright application prospects in piezotronics and theoretical importance in understanding 2D piezoelectricity. Sp 2 piezoelectrics exhibit superior electro-mechanical properties that are deeply rooted in their peculiar sp 2 hybridization: the strong σ-bonds construct a robust honeycomb lattice that can withstand a large deformation [27]; the π-orbitals form the gapped low-energy dispersion with two non-equivalent valleys, mimicking the behavior of massive Dirac fermions [28][29][30], and an elastic deformation can cause the two valleys to drift along opposite directions, then emerging as pseudo-gauge fields [31]. Due to the inversion symmetry breaking (ISB), they are piezoelectric. The ISB strength and the bond polarity of the sp 2 piezoelectrics are determined by the electronegativity difference between their two constituent atoms and are tunable by altering the combination of atoms [27]. Physically, their electronic piezoelectricity is contributed by both σand π-bond electrons. We thus naturally asked (i) what respective roles the two contributions play in the piezoelectricity and (ii) how they vary with the ISB strength or bond polarity. The ISB also endows the sp 2 monolayers with a non-zero MBC, resulting in the valley Hall effect [32]. Given the geometric essence of electronic piezoelectricity, it is unclear how the non-zero MBC manifests itself in the piezoelectricity of sp 2 monolayers. In addition to intrinsic piezoelectricity, an engineered piezoelectricity can also be induced by ISB in centrosymmetric sp 2 crystals such as graphene [33]. Several schemes of this kind for piezoelectric engineering in sp 2 materials have been proposed [34][35][36], but further guidelines beyond this ISB criterion are still lacking.
In this work, by combining the tight-binding (TB) approximation and the modern theory of polarization, we systematically explored the bond-orbital-resolved piezoelectricity of the πand σ-electrons (i.e., the πand σ-piezoelectricity) in sp 2 semiconductors. Firstly, by elaborating a TB piezoelectric model for a prototypical sp 2 crystal, we revealed the different variations in the σand π-piezoelectricity with the ISB-strength (∆). As ∆ varied from positive to negative, the former decreased continuously; whereas the latter increased piecewise and showed a discontinuity at ∆ = 0 due to the criticality of the π-electrons' ground-state geometry during the QPT. Thus, the electronic piezoelectricity near ∆ = 0 was predominated by the π-electrons. Then, the microscopic mechanisms of anomalous π-piezoelectricity and its subtle relation with the valley Hall effect were clarified using an analytical model. Finally, we validated our models by applying them to typical sp 2 materials: hexagonal silicon carbide (h-SiC), Boron-X (h-BX, X = N, P, As, Ab), and a BNdoped graphene (BNG) superlattice. The piezoelectric engineering of the sp 2 materials and relevant strategies are also discussed.

Strain-Dependent TB Hamiltonian for Sp 2 Piezoelectric Crystals
We began with the nearest-neighbor TB Hamiltonian for unstrained sp 2 piezoelectrics shown in Figure 1a, which, by neglecting the spin-orbit coupling, can be written as [37][38][39]: where the operator c † i A ,µ (c i B ,µ ) creates (annihilates) an electron in the atomic orbital µ ∈ s, p x , p y , p z located at the A(B) sublattice site i A (i B = i A + δ I ) with the on-site energy A(B) µ , and t I µ,µ , (I = 1, 2, 3) represents the hopping from orbital µ to orbital µ along the three nearest-neighbor vectors δ 1 = (0, −1)a, δ 2 = ( √ 3, 1)a/2, and δ 3 = (− √ 3, 1)a/2, with a being the unstrained bond length. The on-site energy differences ∆ µ = A µ − B µ reflect the ISB strength and acted as the control parameters in our following piezoelectric model. Under an in-plane strain ε, the nearest-neighbor vectors transform as δ I = δ I + ε · δ I if assuming the clamped-ion approximation [5]. Accordingly, within the Slater-Koster framework [37], the strain-dependent hopping t I µ,µ (ε) can be expressed in terms of the two-center integrals V χ (ε) = V 0 χ exp(1 − β χ δ I /a) as shown in Table 1, where V 0 χ (χ = ssσ, spσ, ppσ, or ppπ) is the two-center integral for unstrained lattice and β χ is the electron-lattice cou- pling parameter [40]. Then, the strain-dependent TB Hamiltonian can be straightforwardly constructed by replacing the hopping terms in Equation (1) with the strain-modified ones (t I µ,µ (ε)).
reflect the ISB strength and acted as the control parameters in our following piezoelectric model. Under an in-plane strain ε, the nearest-neighbor vectors transform as = + I I I ′ ⋅ δ δ ε δ if assuming the clamped-ion approximation [5]. Accordingly, within the Slater-Koster framework [37], the strain-dependent hopping , ( ) I t μ μ′ ε can be expressed in terms of the two-center integrals , or ppπ ) is the two-center integral for unstrained lattice and χ β is the electron-lattice coupling parameter [40]. Then, the strain-dependent TB Hamiltonian can be straightforwardly constructed by replacing the hopping terms in Equation (1) [37]. The direction cosines of vector pointed from the μ -orbital site to the μ′ -orbital site is denoted by  Table 1 vanish. Accordingly, the eight bands of sp 2 piezoelectrics are decoupled into the π-and σ-bands. This can be seen explicitly in the k-space strain-dependent Hamiltonian, which in the basis { , , , , , takes the following block-diagonal form: (2) Figure 1. (a) Undeformed crystal structure of the sp 2 piezoelectrics with three nearest-neighbor vectors δ I=1,2,3 of magnitude a =|δ I | connecting its inequivalent atoms A (blue) and B (red). (b) The corresponding BZs, which are demarcated into two triangular subzones by the Γ-M mirror lines (red) for each valley centered at K and K points. Table 1. The strain-modified nearest-neighbor hopping t I µ,µ (ε) expressed in terms of Slater-Koster two-center integrals V χ (ε) [37]. The direction cosines of vector pointed from the µ-orbital site to the µ -orbital site is denoted byn I µ,µ = (n x , n y , n z ). Other elements can be found by permuting indices.
For a uniform strain reserving the in-plane reflection symmetry; i.e., n z = 0, the π-bond orbital p z does not couple with the orbitals forming the σ-bond µ ∈ s, p x , p y because the corresponding hoppings t I µ ,p z (ε) given in Table 1 vanish. Accordingly, the eight bands of sp 2 piezoelectrics are decoupled into the πand σ-bands. This can be seen explicitly in the k-space strain-dependent Hamiltonian, which in the basis s A , p A y , p A x , s B , p B y , p B x ⊕ p A z , p B z takes the following block-diagonal form: The element of the 6 × 6 matrix H σ (k, ε) for σ-bands reads: where µ m(n) ∈ s A , p A y , p A x , s B , p B y , p B x with m or n ranging from 1 to 6; and 1 = A s , 2 = 3 = A p , 4 = B s , and 5 = 6 = B p are the on-site energies. The remaining 2 × 2 matrix for π-bands is given by: where f (k, ε) = ∑ 3 I=1 t I (ε)e −ik·δ I with t I (ε) = t I p z ,p z (ε) and t I (0) = t for simplicity.

General Formulas for Electronic Piezoelectricity
Piezoelectricity can be regarded as "unquantized" charge pumping driven by an adiabatic lattice deformation ε(t) evolving along an open-path t ∈ [0, T]. The polarization difference accumulated during this evolution is related to the piezoelectric current J via the continuity equation ∆P = T 0 dt J(t). For the leading order, the piezoelectric response is characterized by the rank-3 piezoelectric tensor properly defined as e ijk = (∂ J i /∂ . ε jk ) ε→0 [3].
For the clamped-ion model in the present work, only the piezoelectric current of electrons was relevant in this definition, so henceforth, e ijk will refer to the electronic or clamped-ion piezoelectric tensor unless otherwise specified. In the framework of semiclassical wave packet dynamics, the piezoelectric current is given by [4][5][6][7]: where e is the charge of the electron, the factor of 2 accounts for the spin degeneracy, and Ω n is the PBC built from the strained Bloch eigenstate |u n (k, ε) of the nth occupied valence bands. Consequently, the piezoelectric tensor, by definition, is obtained as [5,6]: with Equation (6) presents a geometric interpretation of the electronic piezoelectricity; the PBC formally resembles the MBC Ω n (k) = i[ ∂ k x u n ∂ k y u n − ∂ k y u n ∂ k x u n ] in the TKNN formula σ H = (e 2 / )∑ n BZ dkΩ n (k)/(2π) 2 [8]. However, they differ from each other in two aspects. Firstly, the MBC is defined on the compact BZ torus (k x , k y ), while the PBC is defined in a non-compact space (k i , ε jk (t)) as topologically different from the BZ torus. This results in the difference between their k-integral σ n H and e n ijk for the isolated nth band: σ n H must be quantized [8], but e n ijk is not necessarily quantized. Secondly, in the presence of time-reversal symmetry (TRS) or inversion symmetry (IS), their k-space distributions show contrary parities. In systems with TRS, the MBC is an odd function of k, leading to a vanishing σ n H [7]; in contrast, the PBC is an even function of k, thus admitting a non-zero e n ijk [4]. For systems with IS, the MBC is an even function of k, while the PBC is an odd function of k, the latter of which requires that e n ijk = 0 and, as expected, excludes the piezoelectricity in such systems. Despite these differences, there can be subtle relations between them when the local patches of (k i , ε jk (t)) and (k x , k y ) can be linearly mapped into each other (see Section 3.2).

Details of DFT Computation
To fit the TB parameters for the sp 2 crystals considered in Section 3.3.1, their band structures were calculated within the density functional theory (DFT) using the VASP package [41]. In the calculation, the exchange-correlation effects were treated at a GGA-PBE level, the electron-ion interactions were described by the projector augmented planewave method, and the energy cutoff for basis-set expansion was chosen to be 600 eV. For structure optimization, the total energy was convergent within 10 −7 eV, and a k-point mesh of 24 × 24 × 1 was used. To exclude the interactions between the neighboring layer, a vacuum layer with a thickness of 20 Å was applied. The processes of fitting TB parameters for these unstrained sp 2 crystals are presented in Section S4 of the Supplementary Materials (Supplementary Materials Section S4).

Results and Discussions
Since piezoelectricity results from ISB, one might intuitively anticipate that it is positively related to the strength of ISB and thus will decrease to zero when the ISB strength decays [24]. However, we will show in this section that this is not the real case for the piezoelectricity in sp 2 crystals.

Bond-Orbital-Resolved Piezoelectricity in Generic Sp 2 Crystals
Armed with the above analysis, we will now use the established formulas to study the electronic piezoelectricity of generic sp 2 crystals. The decoupling between σand π-bands allows splitting of the electronic piezoelectric tensor e ijk into the bond-orbital-resolved con- i,jk (k)/(2π) 2 . If indexing the σand π-valence (conduction) bands by n = 1, 2, 3 (m = 4, 5, 6) and n = 7 (m = 8), respectively, the corresponding PBCs, for convenience of calculation, can be written in the following forms [7]: where E 0 n(m) (k) and |u 0 n(m) (k) are the eigenenergy and eigenstate of the strain-independent Hamiltonian H(k), respectively; v ; and c.c. denotes the complex conjugates. The D 3h symmetry of sp 2 crystals requires that e ijk has only one independent component; e.g., e 222 [15]. Hence, in the following, we will only calculate e σ 222 and e π 222 . Before diving into the calculation, we will briefly explain the parameter setting. Despite the TB parameters for realistic sp 2 crystals that rely on the material details, their similarity in band structures allows for a unified treatment of their piezoelectricity. Instead of modeling any special sp 2 crystal, here we are mainly concerned with the general variation trends in e σ(π) 222 with the ISB strength parametrized by the on-site energy differences ∆ µ=s,p . Roughly speaking, the ∆ µ are positively related to the electronegativity difference between A and B atoms [42]; for a qualitative investigation, it was adequate to assume ∆ p = ∆ s = ∆.
What we wanted to determine was how e σ(π) 222 varied with ∆. To focus on this main motif, we further assumed that the hopping terms did not change with ∆. In our TB calculations, we adopted the typical TB parameters V 0 χ and µ fitted for unstrained graphene in [39] and β χ fitted for strained graphene in [40], then modified the onsite energies as In other words, we took the "ISB-graphene" as a prototype of piezoelectric sp 2 crystals in view of their similar band structure. In doing this, we did not expect our qualitative results to depend on such a parameter choice.
The calculated e σ 222 and e π 222 as functions of ∆ are plotted in Figure 2a. As the figure shows, the magnitude of e π 222 was overall much larger than that of e σ 222 in a quite large range of ∆. This coincided with the previous DFT prediction for the typical sp 2 crystal h-BN that the π-electrons would dominate its electronic piezoelectricity [43]. More interestingly, e σ 222 and e π 222 showed opposite variation trends with the ISB strength: as |∆| decayed, e σ 222 decreased and approached zero as intuitively expected, but e π 222 increased anomalously and finally saturated at a giant finite value. This prominent difference between them implied that sp 2 crystals with a very small ISB would have a quite strong (rather than weak) piezoelectricity contributed almost entirely by the π-electrons, since e σ 222 is neglectable near ∆ = 0. Remarkably, in contrast to the continuous variation in e σ 222 , e π 222 showed an abrupt jump when crossing the peculiar point at ∆ = 0. Since piezoelectricity is just the strain-induced geometric phase shift of the ground-sate wavefunction, such discontinuity reflected the non-analyticity in the geometry of the π-electron ground-state [13].
| | Δ decayed, 222 e decreased and approached zero as intuitively expected, but 222 e increased anomalously and finally saturated at a giant finite value. This prominent difference between them implied that sp 2 crystals with a very small ISB would have a quite strong (rather than weak) piezoelectricity contributed almost entirely by the π-electrons, since 222 e σ is neglectable near 0 Δ = . Remarkably, in contrast to the continuous variation in 222 e σ , 222 e π showed an abrupt jump when crossing the peculiar point at 0 Δ = . Since piezoelectricity is just the strain-induced geometric phase shift of the ground-sate wavefunction, such discontinuity reflected the non-analyticity in the geometry of the π-electron ground-state [13].  In the context of QPT [9][10][11], the non-analyticity of the ground-state geometric phase accompanying gap-closing is a hallmark of the quantum criticality. The critical points marked by gap-closing divide the Hamiltonian parameter space into several regions. Two ground states lying in the same region can be connected by an adiabatic (well-gapped) parameter path along which the observable ground-state properties vary smoothly [9]. However, when the system's Hamiltonian varies across the critical point that separates two different regions, it will undergo a QPT that is characterized by the "critical behavior"; i.e., an abrupt change in properties related to the ground-state geometry such as the Hall conductance [8] and the Born effective charge [13]. The distinct variation trends in e σ 222 and e π 222 in Figure 2a can also be understood within the framework of QPT. To see this more clearly, let us focus on the expression of the PBC in Equations (7) and (8). The denominators 2,22 tends to increase as the occupied E 0 n and unoccupied E 0 m energy levels get close to each other during variation in ∆ and will show singularity when the energy gap between them is closed at certain points ∆. As a result, the corresponding 222 would show critical behavior near this point. To trace the possible singularity in the PBC, in Figure 2b, we plotted the energy gap between the lowest unoccupied and highest occupied levels for the σ-bands and π-bands; i.e., It can be seen that the E σ g for σ-bands retained a finite value over the given interval in the ∆. Therefore, any Hamiltonian H σ (∆) lying in this range was adiabatically connectable to the non-piezoelectric H σ (0) without the appearance of singularity in Ω σ 2,22 . Thereby, e σ 222 should vary smoothly through zero with ∆. Noting that E σ g E π g and hence Ω σ

2,22
Ω π 2,22 over a quite large range of ∆, it is reasonable that e σ 222 would have a relatively minor magnitude compared to the e π 222 shown in Figure 2a.
However, the energy gap E π g = ∆ for π-bands would close at ∆ = 0, across which a QPT occurs between quantum valley Hall states marked by the opposite valley Chern numbers C v = sgn(∆) [12]. When ∆ approaches this critical point, the PBC (Supplementary Materials Section S1) and the MBC [29] of the π-valence band will diverge. The trigonometry terms in the above square brackets are finite (∼ 1) for any k, while at the K or K points k D = (±4π/3 √ 3a, 0), the energy of the π-conduction band in the denominators is E 0 22 and Ω π will diverge as ∝ C v · ∆ −2 for a decreasing |∆|, and then sharply peak around k D like the Dirac function, as shown in Figure 2c. For an sp 2 crystal with TRS, Ω π 2,22 [Ω π ] is an even [odd] function of k and satisfies Ω π 22 [Ω π ](∆). In the |∆|→ 0 limit, the even symmetric and Dirac-function-like distribution of Ω π 2,22 ensures that its k-integral is non-vanishing, while Ω π 2,22 (−∆) = −Ω π 2,22 (∆) requires that e π 222 (0 − ) = −e π 222 (0 + ), thus leading to the abrupt sign-changing of e π 222 when ∆ varies from 0 − to 0 + . Therefore, it was the singularity in Ω π 2,22 near the QPT point that gave rise to the critical behavior of e π 222 shown in Figure 2a. On the other hand, because Ω π (−k) = −Ω π (k), its usual k-integral given by the TKNN formula always vanishes and cannot show any criticality near the critical point at ∆ = 0. To capture the non-analyticity in the ground-state geometry of π-electrons during the QPT, an alternative way is to define an auxiliary quantity called the valley Hall conductance (VHC) [44]: where h = 2π , and the subscripts K and K denote the triangular domain in Figure 1b as delimited by the red Γ-M lines (along which Ω π (k) = 0). It should be noted that in a TB model for π-bands, σ π VH is not quantized for any finite non-zero ∆ [44]. As illustrated in Figure 2d, σ π VH showed quite similar variation trends to those of e π 222 . The discontinuity of σ π VH at the QPT point can also be explained by essentially the same argument employed for that of e π 222 . This similarity between the critical behaviors of σ π VH and e π 222 during the QPT motivated us to further explore their intimate relations in the following.

Valley Model for the Anomalous π-Piezoelectricity
Having discerned the dominant role of π-electrons in the piezoelectricity of sp 2 crystals, we now turn to closely examining the microscopic mechanism underneath the anomalous π-piezoelectricity and its subtle relation to the valley Hall effect as suggested by the TB calculation.

Correlation between π-Piezoelectric Coefficient and VHC
Owning to the sharply peaked distribution of the PBC in the K and K valleys, the main physics of the π-piezoelectricity in sp 2 crystals can be captured by an analytical valley model based on the low-energy effective Hamiltonians. When expanding H π (k, ε) in Equation (4) at k D = (±4π/3 √ 3a, 0), we obtain the massive Dirac Hamiltonian: where v F = 3at/2 is the Fermi velocity; σ 1 , σ 2 , and σ 3 are the Pauli matrices; τ = ± refers to the K or K valley (when appearing in the sub-or superscript of a variable τ = K or K ); and q = (q 1 , q 2 ) is the crystal momentum measured from the k D point. The strain ε emerges as a pseudo-gauge field given by A i = β ppπ /2a∑ jk γ ijk ε jk ; here, the non-zero elements of the rank-3 tensor γ ijk are restricted by the D 3h symmetry to obey −γ 111 = γ 122 = γ 212 = γ 221 = 1 [31]. In real sp 2 crystals, the lattice constant a offers a natural high-energy cutoff of q =|q|≤ Λ ∼ 1/a for the above Hamiltonian, beyond which the Dirac approximation is inapplicable [30]. By starting from H τ (q, ε), repeating the derivation procedure of Equation (9), and exploiting the linear mapping relationship ∂ ε jk → ∂ q i , we find that the PBC of the τ-valley valence band is (see Supplementary Materials Section S2): where is just the MBC for the τ-valley valence band [32], and il is the antisymmetric Levi-Civita tensor. Then, by adding up the integrals of Ω τ i,jk (q) over the K and K valleys, the π-piezoelectric tensor can be obtained as: Thus, in the valley model (labeled with the superscript "v"), e v ijk is directly related to the VHC σ v VH = σ K H − σ K H with σ τ H = (e 2 /h) q≤Λ dqΩ τ (q)/2π. When integrating out the crystal momentum q straightforwardly, we arrive at: where α =|∆/t| measures the polarity of π-bonds [24], and the energy cutoff is set as Λ = √ S BZ /2π to preserve the total number of states in the first BZ with an area of S BZ = 8π 2 /3 √ 3a 2 [30]. When substituting Equation (15) into Equation (14), we finally obtain the explicit expression for the piezoelectric coefficient: Here, the factor −∑ l il γ ljk ensures that the components of e v ijk satisfy the D 3h symmetry; i.e., e v 211 = e v 112 = e v 121 = −e v 222 and e v 111 = e v 122 = e v 212 = e v 221 = 0 [15]. Roughly, β ppπ ∼ l p + l p + 1 is determined by the angular momentum of the involved p z orbitals l p = l p = 1 [42] and therefore can be treated as constant in the following.
We will now discuss the mechanisms of the π-piezoelectricity revealed by Equations (14)- (16). Firstly, e v ijk is an inverse proportional function of the bond length a, and thus it tends to diminish with an elongating a; secondly, as inherited from σ v VH , e v ijk is also a decreasing function of the bond polarity α. Hereinafter, these two trends will be referred to as "the bond length and polarity mechanisms", respectively. When α 1, the VHC is nearly quantized as σ v VH ≈ C v · e 2 /h, and thus the π-piezoelectric coefficient can be approximated as e v 222 ≈ C v · eβ ppπ /2πa. However, this α-independent formula becomes physically unreasonable and would overvalue the magnitude of e v 222 for partially ionic sp 2 crystals such as h-BN (α > 2, see Table 2) [21]. Indeed, recent ab initio studies showed that the piezoelectricity in ISB graphene would decay significantly with an increasing effective α (or ∆) [25].
The valley Chern number C v in Equation (16) reflects the topological aspects of the π-piezoelectricity in sp 2 crystals. Although σ v VH itself is not quantized for a finite ∆ = 0, its jump of when crossing the critical point at ∆ = 0 is exactly quantized. Actually, as pointed out in [11], N 3 is a well-defined topological invariant, and in the spirit of the bulk-edge correspondence, was related to the number of kink states or the zero-energy modes [45] that existed at the interface between the two pieces of sp 2 crystals with an opposite sgn(∆) shown in Figure 3. In this sense, the QPT between ground states of different C v is regarded as marginally topological [46]. In sync with σ v VH , the jump of e v 222 across ∆ = 0 should also be quantized (in units of eβ ppπ /πa) according to Equation (14) and thereby can serve as an alternative signal for probing such topological QPT in sp 2 crystals [12,14]. its jump of when crossing the critical point at 0 Δ = is exactly quantized. Actually, as pointed out in [11], 3 N is a well-defined topological invariant, and in the spirit of the bulk-edge correspondence, was related to the number of kink states or the zero-energy modes [45] that existed at the interface between the two pieces of sp 2 crystals with an opposite sgn( ) Δ shown in Figure 3. In this sense, the QPT between ground states of different v C is regarded as marginally topological [46]. In sync with VH

π-Piezoelectricity as a Hall-Type Response to Pseudo-Electric-Field
From a dynamic viewpoint, the anomalous π-piezoelectric response can be intuitively illustrated by drawing an analogy with the valley Hall effect. Consider now that the strain ( ) t ε is adiabatically time-dependent. The resultant pseudo-gauge potential that couples to the τ -valley electrons through the Peierls substitution will act as an electric field with opposite directions for non-equivalent valleys [47].
This pseudo-electric field is distinct from the real electric field and hence no Hall voltage will be measured, albeit one can define a valley [44]. On the other hand, the piezoelectric current contributed by the τ -valley takes the following form (Supplementary Materials Section S3):

π-Piezoelectricity as a Hall-Type Response to Pseudo-Electric-Field
From a dynamic viewpoint, the anomalous π-piezoelectric response can be intuitively illustrated by drawing an analogy with the valley Hall effect. Consider now that the strain ε(t) is adiabatically time-dependent. The resultant pseudo-gauge potential that couples to the τ-valley electrons through the Peierls substitution q → q − τ A(t) will act as an electric fieldẼ [47]. This pseudo-electric field is distinct from the real electric field E = − . A(t) yielded by the varying magnetic potential A(t), which couples to both valleys with the same signs. It is well known [32] that in the valley Hall effect illustrated in Figure 4a, a longitudinally applied real electric field E will drive the valley-contrasting transverse Hall current of J τ (t) =σ τ H × E (here,σ τ H = σ τ Hẑ = ±σ v VHẑ /2 = ±σ v VH /2). Since J K = −J K , the charge current J K + J K = 0 and hence no Hall voltage will be measured, albeit one can define a valley current J v = J K − J K =σ v VH × E [44]. On the other hand, the piezoelectric current contributed by the τ-valley takes the following form (Supplementary Materials Section S3):

A(t) with opposite directions for non-equivalent valleys
which allows an intuitive reinterpretation of J τ as a Hall-type current driven by the pseudo-electric fieldẼ τ . Differing from the canceling-out by J τ in the charge-neutral valley Hall effect, the strain-induced J τ can add up to result in a non-zero charge current Its accumulation during t ∈ [0, T] yields a bulk polarization of and a potential difference between two edges of the sample shown in Figure 4b, which in an analogous sense can be regarded as the "Hall voltage". To summarize, in a dynamic picture, the anomalous π-piezoelectric effect is a charge-non-neutral counterpart of the valley Hall effect [48]. Such an immediate correlation between them may have important experimental implications for measuring the VHC of sp 2 systems via the piezoelectric effect.
and a potential difference between two edges of the sample shown in Figure 4b, which in an analogous sense can be regarded as the "Hall voltage". To summarize, in a dynamic picture, the anomalous π-piezoelectric effect is a charge-non-neutral counterpart of the valley Hall effect [48]. Such an immediate correlation between them may have important experimental implications for measuring the VHC of sp 2 systems via the piezoelectric effect.

Application to Typical sp 2 Crystal and BNG Superlattice
Given the π-electrons' predominant contribution to the total piezoelectricity compared to that of the σ-electrons, the central concern in this section is the extent to which the piezoelectricity of real sp 2 materials can be determined by their π-electrons. Below, we first apply our π-band models to the typical sp 2 crystals, including h-SiC and BX (X = N, P, As, Sb), and then to the BNG superlattice with D3h symmetry to evaluate their piezoelectricity.

Intrinsic Piezoelectricity of Typical sp 2 Crystals
Let us calculate the π-piezoelectric coefficients of the typical sp 2 crystals including h-SiC and BX using our established models. The TB parameters for the π-bands of these unstrained crystals can be evaluated by fitting their DFT band structures (see Supplementary Materials Section S4). In Table 2, we first list the bond length a and the fitted t and Δ . Then, we calculated their bond polarity α and VHC VH e (corresponding to our models [5]) in the last three columns. To avoid drawing biased conclusions due to the underestimation of band gaps ( Δ ) in the DFT, the GW-corrected TB parameters cited from the extant literature, together with the data calculated from them, are also listed as comparisons in parentheses. Since the GW gaps were systematically larger than the DFT gaps, the values of ( ) 222 v e π calculated from the GW data were smaller than those from the DFT data.   [50], f [51], g [26], h [52], i [53], j [22], k [28], l [54].

Application to Typical sp 2 Crystal and BNG Superlattice
Given the π-electrons' predominant contribution to the total piezoelectricity compared to that of the σ-electrons, the central concern in this section is the extent to which the piezoelectricity of real sp 2 materials can be determined by their π-electrons. Below, we first apply our π-band models to the typical sp 2 crystals, including h-SiC and BX (X = N, P, As, Sb), and then to the BNG superlattice with D 3h symmetry to evaluate their piezoelectricity.

Intrinsic Piezoelectricity of Typical Sp 2 Crystals
Let us calculate the π-piezoelectric coefficients of the typical sp 2 crystals including h-SiC and BX using our established models. The TB parameters for the π-bands of these unstrained crystals can be evaluated by fitting their DFT band structures (see Supplementary Materials Section S4). In Table 2, we first list the bond length a and the fitted t and ∆. Then, we calculated their bond polarity α and VHC σ v VH . By adopting a unified constant β ppπ = 3.3 that was well fitted from graphene and h-BN [20,31], we further evaluated their π-piezoelectric coefficients e π(v) 222 using the TB (valley) model and compared the obtained e π(v) 222 with the clamped-ion DFT results e 222 (corresponding to our models [5]) in the last three columns. To avoid drawing biased conclusions due to the underestimation of band gaps (∆) in the DFT, the GW-corrected TB parameters cited from the extant literature, together with the data calculated from them, are also listed as comparisons in parentheses. Since the GW gaps were systematically larger than the DFT gaps, the values of e π(v) 222 calculated from the GW data were smaller than those from the DFT data. Despite this, we safely concluded from Table 2 that for all these crystals, the π-band contribution e π(v) 222 accounted for a major part of the total electronic piezoelectric coefficient e 222 , while the remaining proportion attributed to the σ-bands had a relatively minor contribution.
For instance, the ratio of e π(v) 222 /e 222 for the partially ionic BN was about 80% (70%), which confirmed the previous DFT calculation results [43]; for the nearly covalent BP and BAs, this ratio could be even higher (up to 95% (85%)), meaning that the piezoelectricity was almost entirely determined by the π-bands.
Moreover, we noticed that as the crystal changed from BN to BSb along the first column of Table 2, the bond polarity α decayed monotonically and was accompanied by an elongation in the bond length a. According to the bond polarity mechanism in Equation (16), the decline in α tended to enhance the π-piezoelectricity via increasing σ v H (α); meanwhile, the concomitant elongation of a restrained this enhancement via the bond-length mechanism. As illustrated in Figure 5, owing to the rapid decay in α, the increase in σ v VH (α) and e v 222 when going from BN to BP was dominated by the bond polarity mechanism. Thereafter, since σ v H (α) gradually saturated, the once-hidden bond-length mechanism came into play through the factor 1/a and partially canceled the feeble growth of σ v H (α), thus leading to the plateau or slight decline in e v 222 for the subsequent BAs and BSb. Therefore, the π-piezoelectricity in these sp 2 crystals was determined by the competition or balance between the above two mechanisms. This in turn would largely govern their overall piezoelectricity considering the subordinate role of σ-electrons. For example, the obvious difference in e 222 between the partially ionic BN and SiC (α > 1) and the nearly covalent BP, BAs, and BSb (α < 1) could be roughly attributed to the apparent bond polarity difference between these two groups.

Engineered Piezoelectricity of BNG Superlattice
Apart from the perfect sp 2 crystals, many engineered sp 2 systems such as graphene on substrates [34], F and H absorbed graphene [35], and BNG [36] can also be piezoelectric. Among these, the BNG superlattice seems to be the most promising platform for piezoelectric engineering due to its experimentally easily tailorable electronic structures [55]. To

Engineered Piezoelectricity of BNG Superlattice
Apart from the perfect sp 2 crystals, many engineered sp 2 systems such as graphene on substrates [34], F and H absorbed graphene [35], and BNG [36] can also be piezoelectric. Among these, the BNG superlattice seems to be the most promising platform for piezoelectric engineering due to its experimentally easily tailorable electronic structures [55]. To validate our models in guiding the engineering of piezoelectricity in sp 2 systems, we took as an example the D 3h -symmetric BNG superlattice shown in Figure 5a, the piezoelectricity of which was first studied using a DFT calculation in [36]. Its primitive vector was p(≥ 2) times as long as that of the pristine graphene. Accordingly, as shown in Figure 5b, the pristine BZ was folded into 1/p 2 of the original one to form a superlattice BZ. The periodically embedded (BN) 3 rings broke the inversion symmetry and opened a band gap of graphene by introducing sublattice-dependent perturbations into it (since the B and N atoms resided in different sublattices). This made the BNG superlattice piezoelectric.
According to the band-folding picture [56], the low-energy-band structures of BNG superlattices fall into two categories depending on whether p is a multiple of 3. In the p = 3n (n is an integer) scenarios in which (BN) 3 -mediated scattering between K and K valleys is suppressed, the BNG superlattice can be approximately treated as ISB graphene crystal within the virtual crystal approximation (VCA) [57][58][59]: both the on-site and hopping energies are averaged separately over each sublattice of the BNG to form a virtual sp 2 crystal; thus, the low-energy dispersion of the deformed BNG could still be described by the effective Hamiltonian in Equation (12) except that the parameters t(x) and ∆(x) depended on the BN concentration x = 3/p 2 . Predictably, their piezoelectricity could be sufficiently evaluated using our valley model. However, when p = 3n, the K and K points of the primitive BZ (red hexagon in Figure 6b) were folded onto the Γ points of the superlattice BZ [56], near which the two valleys were degenerate and the scattering between them was strongly enhanced, thus invalidating the VCA. In this situation, minimally describing the degenerate valleys would require a 4 × 4 valley-coupled Hamiltonian [60]; exploring the piezoelectricity in such systems entails a non-abelian model [7], which was beyond the scope of the current work and therefore was left for future study. In the following, we will focus only on the BNG superlattices with 3 p n ≠ configurations. We used Equation (16) at the mean-field VCA level to estimate their π-piezoelectric coefficients from their corresponding virtual sp 2 crystals. Since  In the following, we will focus only on the BNG superlattices with p = 3n configurations. We used Equation (16) at the mean-field VCA level to estimate their πpiezoelectric coefficients from their corresponding virtual sp 2 crystals. Since E π g = ∆ , the effective ISB strength of the virtual crystals could be derived from the DFT band gaps of the BNG superlattices: |∆(x)|= 4.11x(eV) [36]. If we assumed that their effective hopping varied linearly with x from −2.30 eV for h-BN [49] to −2.80 eV for graphene [39]; i.e., t(x) = −2.80 + 0.50x(eV), then the effective bond polarity was given by α(p) = 12.33/(2.80p 2 − 1.50). Inserting this into Equation (16), we obtained: In Figure 7, the calculated e v 222 (p) when adopting a = 1.42 Å and β ppπ = 3.3 for all of the considered BNG superlattice configurations (∆ > 0) were compared with the clamped-ion DFT results (e 222 (p)). It can be seen that our valley model, albeit only when including the π-band contributions, provided an accurate-enough estimation for the e 222 (p) of the BNG superlattices, especially in the large p (small x) cases. This was not surprising when considering that as the virtual sp 2 crystal's α(p) decayed with the enlarging p, its σ-piezoelectricity became gradually neglectable, so that e 222 (p) ≈ e v 222 (p). Since the bond length was fixed, e v 222 (p) increased with the decaying α(p) solely through the bond polarity mechanism. In the p → ∞ limit, it saturated at an ultrahigh value of e v 222 (∞) = eβ ppπ /2πa = 5.88×10 −10 C/m, which agreed well with the DFT result e 222 (∞) = 5.50 × 10 −10 C/m [36]. Such asymptotic piezoelectricity, in corresponding to a quantized VHC σ v VH (∞) = e 2 /h, (C v = 1), was protected by the marginal valley topology [46] (see also Section 3.2.1), and hence seemed robust against the change of ISB inducers. For example, the further DFT study in [36] showed that after replacing the (BN) 3 rings in BNG superlattices with D 3h hole defects, the calculated value for e 222 (∞) was found to be about 5.60 ± 0.4 × 10 −10 C/m, which was almost unchanged compared to its original value. Arguably, even when the ion relaxation was considered, this π-electron-governed and topologically robust piezoelectricity of the BNG superlattice was still expected. Within the rigid-ion picture, the ionic piezoelectricity of the relaxed BNG superlattice mainly arose from the strain-induced dipole change in the (BN)3 clusters (because the isocharged C 4+ ions contributed no net dipole), which, when averaged over the supercell, would decrease with the enlarging p just like the σ-piezoelectricity. Therefore, in the p → ∞ limit, the ion-relaxation correction to the piezoelectricity primarily affected the πelectrons' response, which was embodied as the renormalization of the electron-lattice coupling parameter ppπ β into rel ppπ β . When adopting an accurate value of 2.66 rel ppπ β = for relaxed-ion graphene [61], the asymptotic piezoelectric coefficient for the relaxed BNG superlattice was calculated as  Comparison between e v 222 (p) calculated using Equation (19) and the clamped-ion DFT results e 222 (p) [36] for the BNG superlattices with p = 3n configurations and ∆ > 0.
Arguably, even when the ion relaxation was considered, this π-electron-governed and topologically robust piezoelectricity of the BNG superlattice was still expected. Within the rigid-ion picture, the ionic piezoelectricity of the relaxed BNG superlattice mainly arose from the strain-induced dipole change in the (BN) 3 clusters (because the iso-charged C 4+ ions contributed no net dipole), which, when averaged over the supercell, would decrease with the enlarging p just like the σ-piezoelectricity. Therefore, in the p → ∞ limit, the ion-relaxation correction to the piezoelectricity primarily affected the π-electrons' response, which was embodied as the renormalization of the electron-lattice coupling parameter β ppπ into β rel ppπ . When adopting an accurate value of β rel ppπ = 2.66 for relaxed-ion graphene [61], the asymptotic piezoelectric coefficient for the relaxed BNG superlattice was calculated as e v−rel 222 (∞) = 4.74 × 10 −10 C/m, which again agreed well with the DFT result of e rel 222 (∞) = 4.50 × 10 −10 C/m [36]. This experimentally measurable piezoelectric coefficient was about one order larger than that of the F and H absorbed graphene C 2 HF (e rel 222 = 0.63 × 10 −10 C/m) [35] and three times larger than that of the h-BN monolayer (e rel 222 = 1.38 × 10 −10 C/m) [15]. Based on the above discussion, we concluded that, as in the cases of perfect sp 2 crystals, the piezoelectricity in the BNG superlattice was also controlled to a large extent by the π-electrons.
We will now discuss the feasibility of π-piezoelectric engineering in real sp 2 materials and propose some suggestions for this purpose based on our results. Compared with strong σ-bonds, low-energy π-bonds are more sensitive to structural, physical, or chemical modifications [34][35][36]55]. These external perturbations will inevitably reshape the ground-state geometry of π-electrons and thereby affect their piezoelectric response. Such a tunability of π-electrons combined with their decisive role in piezoelectricity makes the efficient engineering of piezoelectricity in sp 2 materials possible by predesigning their π-band structures. Because the anomalous piezoelectric response of π-bands stems from its non-trivial valley geometry (as reflected by the non-zero VHC), to maximally exploit the latent π-piezoelectricity, the first suggestion is that the engineering schemes should be designed to avoid destroying the valley structures. This in turn may explain why C 2 HF has a one-order-smaller piezoelectric coefficient than that of the BNG [35]: the absorbent F and H atoms fix all the non-local π-electrons of graphene into strong C-F and C-H σ-bonds and thereby push their energy far away from the Fermi-level, even beyond the C-C σ-bands [62], and then consequently destroy the low-energy valleys. Upon reserving the valley structure of π-bands, higher π-piezoelectricity would benefit from a shorter bond length or weaker bond polarity in the VCA sense. With the shortest bond lengths in the sp 2 family, graphene (a = 1.42 ) and h-BN (a = 1.44 ) seem to be the ideal starting crystals for π-piezoelectric engineering, as suggested by the bond-length mechanisms. Once the starting crystal is selected, external modifications generally cause hardly any dramatic change in its bond length, and thus manipulate its π-piezoelectricity mainly via the bond-polarity mechanisms. Therefore, another suggestion based on this consideration is that the π-piezoelectricity hidden in the starting crystals can be further released by engineering strategies that lower its effective VCA bond polarity (or equally narrow the band gap |∆| relative to t). For example, after hybridizing with carbon to form borocarbonitrides (B x C y N x ) [23], the band gap of the h-BN monolayer was dramatically narrowed, and as a result its relaxed-ion piezoelectric coefficient could be significantly enhanced from the original e rel 222 = 1.38 × 10 −10 C/m up to e rel 222 = 5.00 × 10 −10 C/m for the BC 2 N structures. Before closing this section, we will briefly discuss how the piezoelectricity in sp 2 monolayers can be observed. For crystals with a moderate band gap such as h-SiC, the piezoelectricity is usually observed via a transport measurement [16] in which the oscillating piezoelectric current is detected as a response to the cyclic deformation. However, this method is not suitable for sp 2 monolayers such as h-BN, the large band gap of which impedes the measurement of current. As an alternative, we highly recommend an innovative approach recently proposed in [63]; i.e., observing the piezoelectric effect by detecting the piezo charges induced by inhomogeneous in-plane deformations ε(r). In this scheme, the density of piezo charges is given by ρ piezo = −∇ r · P = e 222 [2∂ x ε xy + ∂ y (ε xx − ε yy )]. For a given strain configuration, e 222 can be inferred from the distribution of ρ piezo measured using electrostatic force microscopy (EFM) [19]. For example, in the simplest situation in which the lattice deformation has a constant strain gradient c along the y direction ε(r) = (0, cy + c 0 ), the piezoelectric coefficient can be easily determined from the uniform piezo charge density as e 222 = −ρ piezo /c.
In real EFM experiments, due to the in-plane polarizability of sp 2 monolayers [64], the piezoelectric field will be further screened, thus leading to the partial compensation of ρ piezo . In other words, the total charge density under inhomogeneous strains should contain the piezo-and screening-induced parts: ρ = ρ piezo + ρ ind . The latter is determined by ρ ind = χ 2D ∇ 2 r ϕ(r), where ϕ(r) is the piezo-charge-induced electric potential, and χ 2D is the 2D polarizability of the sp 2 monolayers [63]. Using the Poisson equation ∇ 2 r ϕ = −4πρ, the piezo charge density is related to the piezo potential via ρ piezo = −(χ 2D + 1/4π)∇ 2 r ϕ. In the simplest example considered above, one can further extract e 222 from the measured ϕ(y) as e 222 = (χ 2D + 1/4π)∂ 2 y ϕ(y)/c. Since the observable piezo potential ϕ(y) is reduced by the compensation charges ρ ind , ignoring this screening effect will lead to an incorrectly extracted piezoelectric coefficient e * 222 = ∂ 2 y ϕ(y)/4πc. Therefore, when observing the piezoelectricity of sp 2 monolayers via the EFM measurement, their in-plane polarizability must be carefully treated.

Conclusions
In summary, in this work, the bond-orbital-resolved electronic piezoelectricity in sp 2 -hybridized monolayer semiconductors was systematically investigated by combining the TB method and the geometric phase theory of polarization. We revealed in the TB calculation that their πand σ-piezoelectricity showed contrasting variations trends with the ISB strength ∆. Unlike the continuous decrease in the subordinate σ-piezoelectricity with a decaying |∆|, the predominant π-piezoelectricity increased piecewise as |∆|→ 0 and displayed critical behavior at the QPT point ∆ = 0 due to the non-analyticity of the π-electrons' ground-state geometry near this point. By focusing on the anomalous piezoelectric response of π-electrons, we further related the π-piezoelectricity to the valley Hall effect and reinterpreted it as a Hall-type response to the strain-induced pseudo-electric field in the low-energy valley model. With the help of this analytical model, we also identified the bond-length and bond-polarity mechanisms that underlie π-piezoelectricity and clarified its topological aspects. The validity of our theoretical predictions that the piezoelectricity of these materials is mainly dominated by the anomalous response of π-electrons was quantitatively justified by applying the π-band models to the typical sp 2 crystals h-SiC and BX, as well as the BNG superlattice. Our investigation not only deepens the understanding of piezoelectricity in real sp 2 materials, but also provides guidelines for tailoring their piezoelectricity, thus opening doors to π-electron piezoelectric engineering in these systems.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ma15217788/s1, Figure S1: Comparison between the TB band structure (red lines) and the DFT band structure (blue circles) for the h-SiC monolayer. Here, the Fermi energy was set at E F = 0. The supporting information contains the detailed steps needed to obtain Equation (9), Equation (13), and Equation (17), as well as the processes for fitting TB parameters for the considered sp 2 crystals.