Long range piezoelectricity effects in van der Waals heterobilayer systems beyond 1000 atoms

Twist angle is a relevant design and control component for the piezoelectric coefficients of van der Waals (vdW) heterostructures. This theoretical work assesses in high detail the impact of the twist angle on the piezoelectricity of two-dimensional (2D) heterobilayer systems. We expand the density-functional based tight-binding method to predict the piezoelectric coefficients of twisted and corrugated 2D heterobilayer structures with more than 1000 atoms. We showcase the method on hexagonal III–V/transition metal dichalcogenide vdW heterosystems. Our calculations yield a periodic relationship between the in-plane piezoelectric coefficients and the corresponding twist angles, indicating the tunability of the in-plane piezoelectricity. In contrast, the out-of-plane piezoelectricity is not twist angle dependent, but nonlinearly changes with the average interlayer distance.


Introduction
Piezoelectric two-dimensional (2D) materials attract growing interests in the fields of piezotronics and nanoelectromechanical systems, such as actuators, transducers, and energy harvesters [1][2][3][4][5][6][7][8][9][10].The piezoelectric coefficients of bilayer systems have been found to significantly exceed the sum of the respective monolayers [11,12].There is experimental evidence that indicates twist angles offer relevant design and control aspects for engineering the piezoelectricity of van der Waals (vdW) bilayer systems.For instance, some stacking Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.orders of hexagonal boron nitride (h-BN) bilayers induce outof-plane polarizations, and some twisted bilayer graphene was found to produce flexoelectricity [13,14].
This work assesses the impact of the twist angle on the piezoelectric coefficients of vdW heterobilayers on the examples of h-BN/h-BP (hexagonal boron phosphide) and h-BN/MoS2 (molybdenum disulfide) heterobilayer systems.It is found that the analytical description of the angle-dependence of the in-plane piezoelectric coefficient is correct for idealized vdW heterobilayer systems but deviates increasingly from the density-functional based tight-binding (DFTB) calculations the stronger the coupling strength between the layers is.
To obtain accurate and reliable predictions, all possible internuclear repulsion energies and long-distance effects caused by twist angle, corrugation, and strain, are carefully included in terms of well-converged observables vs. supercell sizes.However, the numerical load of density-functional theory (DFT) calculations prohibits expanding the supercells of heterobilayer systems with arbitrary twist angles to full convergence.Third order density-functional based tightbinding (DFTB3) theory is capable of modeling systems containing thousands of atoms and is therefore applied in this work [15,16].Our extended DFTB approach allows the calculations of the piezoelectricity in fully relaxed, twisted, and corrugated vdW heterobilayer systems.The DFTB parameters are optimized to accurately reproduce the DFT bandstructures and piezoelectric coefficients of the respective monolayers and untwisted heterobilayer using Heyd-Scuseria-Ernzerhof (HSE06) exchange-correlation functional as implemented in the Vienna ab initio simulation package (VASP) code [17][18][19][20].A periodic relationship between twist angle and the in-plane coefficients is observed in both h-BN/h-BP and h-BN/MoS2 heterobilayer simulations.Moreover, the calculated e 11 and e 22 of h-BN/h-BP systems deviate further from the symmetry predicted by the analytical model, which is an indication that stronger interlayer coupling exists in such systems.In contrast, the out-of-plane piezoelectric response appears to be almost independent of the twist angle.Instead, the out-ofplane piezoelectricity is exponentially dependent on the interlayer distance, which makes oscillations in the interlayer distance and corrugation strongly impact the out-of-plane piezoelectric coefficients.This out-of-plane sensitivity could be exploited for sensor and amplification applications.All these findings shed light on the design and optimization of the piezoelectricity in vdW heterostructures.
The rest of the paper is organized as follows: The next section details how the piezoelectric coefficients are extracted from DFT and DFTB calculations.This is followed by an explanation of the DFTB parametrization method.The result section lays out the transferability of DFTB parameters, the convergence of observables vs supercell sizes, and the results of in-plane and out-of-plane piezoelectric coefficients as a function of twist angles.

Geometry
The following method is applied on both DFT and the corresponding DFTB simulations.Monolayers (h-BN, h-BP, and MoS2) are constructed from the primitive hexagonal unit cell.The lattice constants and internal coordinates are allowed to fully relax until the maximum force exerted on each atom are smaller than 0.01 eV Å −1 .The initial commensurate supercell of the heterobilayer is created by stacking the separately optimized monolayers on top of each other with artificial boundary condition using the CellMatch code [21] with the dipole vectors of the two layers right above each other and pointing toward the +x direction.To make a twist angle of θ in the initial supercell setup, the upper layer is rotated by the angle of θ with respect to its original orientation as shown in figure 1.Then the bilayer structure is relaxed until the maximum force on each atom is less than 0.01 eV Å −1 .As detailed in later paragraphs, DFTB simulations allow increasing the size of the supercells until convergence of all observables is found.Typically, a supercell containing several hundreds of atoms is required for an accurate prediction on the heterobilayer piezoelectric coefficients.To ensure the accuracy of material predictions, the strain between different monolayers is relaxed with increasing supercell sizes until all variations of observables between different iterations are less than 5%.Detailed findings on the supercell sizes required for convergence are discussed in the result section.The atomic coordinates of a relaxed bilayer structure are expected to deviate from their original values due to intralayer and interlayer interactions.Corrugation is defined by the average deviation of the distance between nearest neighbors of different layers from the average interlayer distance (due to relaxation), which is derived by averaging the distances from each atom of a given layer to its three nearest atoms in the opposing layer.

Piezoelectric coefficient calculation
In Voigt notation the general piezoelectric tensor is represented as a 3 × 6 matrix e ij where 1 ⩽ i ⩽ 3 and 1 ⩽ j ⩽ 6 [22].Unperturbed h-BN, h-BP, and MoS2 monolayers possess a 6m2 point-group symmetry, leading to the existence of only one independent piezoelectric coefficient.According to the geometric phase approach [23], piezoelectric coefficients are solved by the differential change of 2D polarization with respect to uniaxial strain in a range from −0.015 to 0.015 in steps of 0.005 (see figure 3 in the result section).Atomic positions are relaxed at each strain step to generate the relaxed-ion piezoelectric coefficients.The clamped-ion piezoelectric coefficients are calculated with the same geometric phase approach procedure but without relaxing atomic positions at each strain configuration [17][18][19][20].The computational demands of DFT significantly limit the expansion of supercells for heterobilayer systems with arbitrary twist angle to achieve comprehensive convergence.To obtain a reliable piezoelectricity prediction, the geometry model of twisted bilayer must be sufficiently large until the error of each observable is smaller than 5%.

DFTB3 method
We summarize known aspects of the DFTB3 method for later reference in the parameterization section.The DFTB method solves Kohn-Sham equations with low enough computational load to allow for simulations with more than a 1000 of atoms [16].The third-order Taylor expansion of the DFT total energy functionals around a chosen reference density ρ 0 (r), which is the superposition of electron densities of neutral atoms, can be expressed as [16,24,28,29]: .
Note the first order expansion term is cancelled by elements of the second order.The first order term E DFTB1 [ρ 0 ] can be further divided into two parts [16]: where n i represents the occupation number of a molecular orbital ψ i and V rep ab denotes the distance-dependent repulsive potential between atoms a and b, which is determined as explained in the DFTB parameterization section.ρ a 0 and ρ b 0 are reference densities of the respective atoms, and R ab is their interatomic distance.The molecular orbital ψ i (r) can be written as a linear combination of pseudoatomic orbitals: where ϕ µ is the basis function of orbital µ centered at atom a located at R a and c µi is the basis set coefficient.Pseudoatomic orbitals ϕ µ are determined by solving the Kohn-Sham equations with a confining potential Vconf,µ (r): where T is the kinetic energy operator and V is composed of external potential, Coulomb repulsion, and exchange correlation energy [16,[24][25][26].The confining potential operator Vconf,µ reads on its diagonal: where the confinement radius r 0,µ and the exponent σ µ are fitted parameters as discussed in the next section [16,[24][25][26].We use as additional fitting parameters the orbital energies of the free spherical atoms ϵ free µ , which is defined in [24].Note that the off-diagonal elements of the Hamiltonian are computed by [16,[24][25][26]: The second-order E γ and third-order E Γ terms are higher order corrections to the total electronic energy in selfconsistent-charge DFTB (termed as 'DFTB2' and 'DFTB3' in [16]).They are associated with chemical hardness and its derivative respectively, which are computed directly from DFT (see [16,[24][25][26] for more details).

DFTB parameterization
All fitting parameters previously mentioned are considered to be material-and orbital-dependent.They depend on atom types, in this work Mo, S, B, and N atoms, as well as the angular momentum of the respective atomic orbital (s, p, and d).We fitted all parameters so that DFTB results for the bandstructures and piezoelectric coefficients agree with the respective results of DFT calculations with HSE06 functionals solved by VASP [17][18][19][20].Please note that an accurate prediction of the piezoelectric coefficients requires even deep lying valence bands of DFTB agree very well with DFT results (see figure 2).
Repulsive potentials are dependent on atom types only and typically assumed to be transferable between different materials [27].In this study, a two-step fitting process is applied to obtain all the repulsive potentials.First of all, the preliminary repulsive potentials are extracted from the energy differences between DFT HSE06 data and DFTB calculations.Then we fit all the relevant repulsive potentials simultaneously with the genetic algorithm until the DFT HSE06 results are reproduced.To cover different types of interactions, including intralayer and interlayer interactions, repulsive potentials are fitted to a range longer than the distance between nearest neighbors.To ensure transferability, the reference systems considered in the repulsive potential parameterizations include all respective monolayers, heterobilayers, and all relevant atom pairs (dimers).The fitting target is to minimize the total energy differences between DFT and DFTB calculations of all reference systems simultaneously.The quality of the repulsive potential is sensitive to the choice of division points and the smoothness of the interpolation between them.Both have to be chosen carefully to avoid deviations of the DFTB-predicted piezoelectricity from DFT results.To avoid unphysical forces exerted on atoms, repulsive potentials are described by a fourth-order spline function, which is continuously differentiable up to the second-order derivative in each interval [27].In addition to the total energy, we impose the constraint that the fitted repulsive potentials reproduce the ionic contribution to the piezoelectric tensor of the monolayers.

Parameter transferability
The DFTB parameters are chosen to accurately reproduce DFT monolayer and untwisted heterobilayer bandstructures of VASP HSE06 calculations.For an accurate and reliable piezoelectric coefficient prediction with DFTB, deep lying valence bands of DFT calculations have to be faithfully reproduced (illustrated in figure 2), since the polarization mostly depends on the ionic cores and the occupied valence bands.Our results show that both the DFTB-calculated ionic and the electronic contributions to the piezoelectric coefficient e 11 agree with HSE06 DFT data and are transferable to all relevant strain constellations, as shown in figure 3.

Convergence vs. supercell size
For twisted heterobilayer structures, a well-converged simulation requires sufficiently large initial supercells.Figure 4 indicates that the interlayer distance gradually converges with respect to the supercell size.Before convergence, the average interlayer spacing can vary significantly and rapidly with the supercell size.Typically, a supercell containing several hundreds of atoms is necessary for well-converged piezoelectricity calculations (see figure 4).A geometry model with more than 1000 atoms can be required for the systems with long commensurate unit cell lengths.This convergence affects all observables including piezoelectric coefficients and elastic constants.Therefore, for all the simulations in this work the supercells have been extended until convergence was achieved, which resulted in systems of at least 600 atoms.To illustrate the underlying lack of periodicity within ranges of 1 nm, figure 5 shows a zoom-in of a converged h-BN/h-BP system with a twist angle of 10 degrees.The in-plane interlayer distance and corresponding corrugation landscape shows no short-range order.Average interlayer distance and its standard deviation of h-BN/h-BP heterobilayer structures with various twist angles as a function of supercell size.Typically, a supercell containing at least 600 atoms is necessary for well-converged results., with e hBN 22 (θ = 30 • ) the maximum amplitude of e hBN 22 of hBN in direction 2 which equals the amplitude of e hBN 11 .The calculated results of the twisted bilayer system deviate from the idealized results for all twist angles depending on, for instance, the charge transfer between the layers, the break of inversion symmetry, corrugations, and nonlinear, twist angle dependent interference effects of the same: The local corrugation that is specific to a given twist angle determines the local strain fields and the distance between the layers.Local strain directly enters the piezoelectric coefficient while the distance between the layers determines the strength of interlayer interaction and the charge transfer between the layers and thus the local polarization.The twisted h-BN/h-BP system shows a larger deviation from this idealized result than h-BN/MoS2 which indicates a stronger interlayer coupling in h-BN/h-BP.If the distance between the monolayers is artificially increased, the analytical forms of the in-plane piezoelectric coefficients are recovered.Furthermore, it is important to note that the 120 degrees periodicity of the idealized system is still maintained  when realistic interlayer coupling is considered.Nevertheless, figures 6 and 7 show the tunability of the piezoelectricity of heterobilayers as a function of twist angle.

Out-of-plane piezoelectric coefficients
Out-of-plane piezoelectric response arises from asymmetric charge distribution and broken inversion symmetry  in heterobilayers along the z-direction.This is illustrated in figures 8 and 9 for the h-BN/h-BP and h-BN/MoS2 heterostructures, respectively.In contrast to the in-plane coefficients, e 33 is mostly independent of the twist angle: Twisting the layers is a rotation around direction 3 of the laboratory system.Therefore, the projection of polarization in direction 3 does not immediately change with twist angles.However, the average interlayer distance is twist angle dependent (see figure 4).Since the charge transfer between layers depends on the interlayer distance, the out-of-plane dielectric response is still twist angle dependent.With the irregular corrugation for given twist angles (figure 5) and the rapid oscillations of supercell-size averaged interlayer distances (see figure 4) the out-of-plane piezoelectric coefficient shows a rather irregular twist angle dependence in figures 8 and 9.
Additionally, we performed a series of numerical experiments to investigate the correlation between the out-of-plane piezoelectric coefficient and interlayer distance.The interlayer spacing is shifted upward by 10% and downward up to 30% after structural relaxation.We find that e 33 changes nonlinearly with respect to this interlayer shift, and its absolute value increases by a factor of 10x compared to the equilibrium distance (see figure 10).

Conclusion
We introduce an efficient DFTB-based approach that allows converged piezoelectric coefficient predictions of relaxed, twisted heterobilayer systems beyond 1000 atoms on a single compute node.The prediction of twist angle dependent piezoelectric coefficients of heterobilayers converges with supercell sizes of around 600 atoms only.Our results unveil controllable in-plane piezoelectricity in both h-BN/h-BP and h-BN/MoS2 heterobilayer structures.The corresponding outof-plane piezoelectric response, on the other hand, strongly depends on the interlayer distance.That distance fluctuates with pronounced monolayer corrugations, which do not show systematic twist angle dependence.Therefore, in first order, e 33 .is independent of the twist angle and mostly constant.However, the out-of-plane piezoelectric response nonlinearly depends on the interlayer distance and can change by a factor of 10x with a reduction of interlayer distance of 30%.

Figure 1 .
Figure 1.Schematic of h-BN/h-BP heterobilayer structure.The angle θ between the dipole vectors of upper h-BN layer and lower h-BP layer (gradient black arrows) is the twist angle.

Figure 2 .
Figure 2. A comparison of the bandstructures of the reference heterobilayer system.The squares are calculated by DFTB with parameterized DFTB parameters.The gray lines are solved by DFT using HSE06 functional, which is implemented in the VASP software.Our DFTB calculations accurately reproduce the deep lying valence bands as well as the first few conduction bands.

Figure 3 .
Figure 3. Applied uniaxial strain along the x-axis results in a change in polarization along the same axis for monolayer h-BN.The DFTB-calculated data (black squares and triangles) agrees with the results solved by HSE06 DFT (gray solid and dotted lines), including both ionic contribution (dotted) and electronic contribution (solid).

Figure 4 .
Figure 4. Average interlayer distance and its standard deviation of h-BN/h-BP heterobilayer structures with various twist angles as a function of supercell size.Typically, a supercell containing at least 600 atoms is necessary for well-converged results.

Figures 6
Figures 6 and 7 show the DFTB-calculated in-plane piezoelectric coefficients e 11 (black symbols) and e 22 (grey symbols) of twisted h-BN/h-BP and h-BN/MoS2 heterobilayer structures as a function of twist angle, respectively.In both figures, results of an analytical formula for

Figure 5 .
Figure 5. Zoom into the interlayer spacing of the relaxed h-BN/h-BP heterobilayer system twisted at θ = 10 degrees and its corresponding corrugation landscape projection on the xy-plane.

Figure 6 .
Figure 6.In-plane piezoelectric coefficients e 11 (black) and e 22 (gray) of twisted h-BN/h-BP heterostructures vs. twist angle.The symbols and the dotted lines are the DFTB results and the analytical formulas derived from isolated monolayers, respectively.

Figure 7 .
Figure 7. In-plane piezoelectric coefficients e 11 (blacke) and e 22 (gray) of twisted h-BN/MoS 2 heterostructures vs. twist angle.The symbols and the dotted lines are the DFTB results and the analytical formulas derived from isolated monolayers, respectively.

Figure 8 .
Figure 8.The solid line represents out-of-plane piezoelectric coefficient e 33 of h-BN/h-BP heterobilayer structures vs. twist angle.The dotted line shows the average interlayer spacing of the same system as a function of twist angle.Similar fluctuating behavior illustrates that the piezoelectric coefficient e 33 strongly depends on the interlayer distance.

Figure 9 .
Figure 9.The solid line represents out-of-plane piezoelectric coefficient e 33 h-BN/MoS 2 heterobilayer structures vs. twist angle.The dotted line shows the average interlayer spacing as a function of twist angle.Similar fluctuating behavior illustrates that the piezoelectric coefficient e 33 strongly depends on the interlayer distance.

Figure 10 .
Figure 10.Out-of-plane piezoelectric coefficient e 33 of MoS2/h-BN as a function of interlayer distance and the twist angle.

Figure A1 .
Figure A1.Conventional unit cell of a hexagonal 2D monolayer oriented with anion-cation bonds along laboratory direction 1.(a) The ideal, unstrained unit cell does not host any net polarization P 1 in direction 1.(b) When the structure of (a) is strained in direction 1, a finite net polarization P 1 results from the sum of all local polarizations (e.g.p upper 1 , p middle1

Figure A2 .
Figure A2.The structure of figure A1(a) rotated by 30 degrees.(a) The net polarization of the unstrained structure vanishes.(b) When the structure of (a) is strained in direction 1, the sum of all local polarizations has no net projection in direction 1, i.e.P 1 = 0. To illustrate the mutual canceling of projections of pairs of local polarizations in direction 1, some representative local polarization pairs (e.g.p LR , p RL , p BR , and p BL ) are depicted with red arrows.