††thanks: Brahmanandam Javvaji and Bo He contributed equally to this
work.††thanks: Brahmanandam Javvaji and Bo He contributed equally to this
work.

# High Flexoelectric Constants in Janus Transition-Metal Dichalcogenides

Brahmanandam Javvaji  Bo He  Institute of Continuum Mechanics, Leibniz
Universität Hannover, Appelstr. 11, 30167 Hannover, Germany  Xiaoying Zhuang
Corresponding author<EMAIL_ADDRESS>Institute of Continuum
Mechanics, Leibniz Universität Hannover, Appelstr. 11, 30167 Hannover, Germany
College of Civil Engineering, Tongji University, 1239 Siping Road, 200092
Shanghai, China  Harold S. Park  Corresponding author<EMAIL_ADDRESS>Department of Mechanical Engineering, Boston University, Boston, Massachusetts
02215, USA

###### Abstract

Due to their combination of mechanical stiffness and flexibility, two-
dimensional (2D) materials have received significant interest as potential
electromechanical materials. Flexoelectricity is an electromechanical coupling
between strain gradient and polarization. Unlike piezoelectricity, which
exists only in non-centrosymmetric materials, flexoelectricity theoretically
exists in all dielectric materials. However, most work on the
electromechanical energy conversion potential of 2D materials has focused on
their piezoelectric, and not flexoelectric behavior and properties. In the
present work, we demonstrate that the intrinsic structural asymmetry present
in monolayer Janus transition metal dichalcogenides (TMDCs) enables
significant flexoelectric properties. We report these flexoelectric properties
using a recently developed charge-dipole model that couples with classical
molecular dynamics simulations. By employing a prescribed bending deformation,
we directly calculate the flexoelectric constants while eliminating the
piezoelectric contribution to the polarization. We find that the flexoelectric
response of a Janus TMDC is positively correlated to its initial degree of
asymmetry, which contributes to stronger $\sigma-\sigma$ interactions as the
initial degree of asymmetry rises. In addition, the high transfer of charge
across atoms in Janus TMDCs leads to larger electric fields due to
$\pi-\sigma$ coupling. These enhanced $\sigma-\sigma$ and $\pi-\sigma$
interactions are found to cause the flexoelectric coefficients of the Janus
TMDCs to be several times higher than traditional TMDCs such as MoS2, whose
flexoelectric constant is already ten times larger than graphene.

##  I Introduction

Atomically thin two-dimensional (2D) materials have a variety of unique
physical properties that have made them attractive for many different
applications [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11]. An emerging area of interest
for 2D materials is electromechanical coupling, due to the desire to
miniaturize sensors and actuators to the micro and nanoscales. The most widely
studied electromechanical coupling mechanism is piezoelectricity, which has
also been investigated for 2D materials, including graphene, hexagonal boron
nitride (h-BN), transition metal dichalcogenides (TMDCs) and many others [12;
13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24]. We note that the majority of
this work has focused on the in-plane electromechanical properties.

TMDCs exhibit a unique three-layer atomic arrangement, where a metallic (M)
atom symmetrically bonds with two chalcogenide atoms (X and X) in the out-of-
plane direction. This structure makes them as good candidates for in-plane
stretching based piezoelectric materials [16], though a recent experimental
study measured how an out-of-plane indentation induced an electrical response
in MoS2 [25]. This deformation induced significant changes in the bond lengths
between M and X atoms, resulting in a strain gradient, which potentially
enabled flexoelectricity, which is a form of electromechanical coupling in
which electrical polarization is generated due to strain gradients [26; 27;
22; 28; 29; 30]. Recent work from the present authors [31] investigated the
bending flexoelectricity in various 2D materials, including graphene, graphene
allotropes, nitrides, graphene analogs of group-IV elements and TMDCs. That
study found that MoS2 (an MXX material) has a flexoelectric coefficient that
is ten times larger than graphene due to enhanced charge transfer resulting
from asymmetrical bending-induced changes in the M-X bond lengths.

Because of the impact of asymmetry in the bending-induced changes in the M-X
bond lengths in enhancing the flexoelectric properties of MoS2, we focus here
on potential flexoelectric effects in another class of TMDCs, the so-called
Janus TMDCs, which introduce an asymmetry in the MXX by replacing one of the X
layers of atoms with a different chalcogenide atom Y, resulting in an MXY
structure and intriguing physical properties [32; 33; 34; 35; 36], as recently
reviewed [37]. The changes in the atomic mass and electronic configurations of
X and Y atoms in MXY generates an out-of-plane as well as in-plane dipole
moment, which are absent in the MXX structure, due to the non-cancellation of
interactions between M-X and M-Y. Density functional theory (DFT) simulations
have been used to report high piezoelectric coefficients for various materials
in the Janus TMDCs family compared to conventional TMDCs [38], while a recent
experimental report reported the out-of-plane piezoelectric response for a
Janus MoSSe monolayer [39]. Furthermore, multi-layer Janus TMDCs have shown
very high out-of-plane piezoelectric coefficients due to the increase in
vertical dipole moments [38]. Overall, an electromechanical imbalance exists
due to the element changes between X and Y in MXY, which should facilitate
asymmetric deformation, strain gradients, and thus flexoelectricity, while
also impacting the induced dipole moments in the Janus TMDCs.

Various works have investigated the piezoelectric properties of Janus TMDCs.
For example, the tensile load induced in-plane piezoelectric coefficient
($d_{11}$) of the Janus TMDCs was found to be several times larger than their
out-of-plane shear piezoelectric coefficient ($d_{31}$) [38; 40]. where the
$d_{31}$ coefficient is absent in most other 2D materials, including MoS2, due
to reflection symmetry[16]. Additionally, the out-of-plane piezoelectric
coefficient ($d_{33}$) for monolayer MoSSe under compression was reported to
be several orders of magnitude smaller than $d_{11}$ [39]. Under compression,
large values for $d_{33}$ were reported for multi-layer Janus TMDCs, where the
monolayers are arranged such that the induced polarization does not cancel
out. However, for the multilayer structure, the out-of-plane elastic constants
increase with increasing number of layers, which requires an increasing amount
of force to induce deformation in the vertical direction [41], and thus is not
an optimal choice for electromechanical energy conversion. In contrast, Janus
TMDCs, as with most atomically thin 2D materials, are significantly easier to
bend rather than stretch[38; 42; 25]. The relative ease of bending 2D
materials, coupled with the fact that bending intrinsically generates a strain
gradient, indicates significant electromechanical energy conversion potential
for bending flexoelectricity as compared to in-plane piezoelectricity for
Janus TMDCs.

In order to compare the energy conversion potential between flexoelectricity
and piezoelectricity, it is essential to separate the contribution of the
flexoelectric response. However, the effect of flexoelectricity is measured in
terms of effective piezoelectricity in practice. For example, the experimental
work [25] determined the flexoelectric coefficient through a relationship with
the measured piezoelectric coefficient under assumptions of small length
scales and linear electric fields. Additional example [43] reported the
effective out-of-plane piezoelectricity from MoS2. The experimental study
reported the piezoelectric responses of a corrugated TMDC [44; 45], where
electrical polarization is mainly due to the local strain gradients that
govern flexoelectricity. A recent study addressed the effective
piezoelectricity from the effect of flexoelectricity in a non-piezoelectric
material [46]. Large-scale experimental approaches like axial stretching or
radial compression of a cylindrical rod wrapped with a non-piezoelectric
material could enable direct measurements of the flexoelectric coefficients
[47]. However, experimental methods to find the flexoelectric constants at the
nanoscale are unresolved. First-principle simulations can isolate the
flexoelectric effect by assuming unstable wrinkles in TMDCs [48]. The authors
previous work [31] provides a mechanical bending deformation which enables
calculation of the flexoelectric response by removing the out-of-plane
piezoelectric contribution to the total polarization.

In this work, we coupled classical molecular dynamics (MD) with a charge-
dipole model [49; 50] to investigate the bending flexoelectric response of the
Janus TMDC family. We first validate the simulation methodology with respect
to previous, DFT-calculated in-plane piezoelectric coefficients [38]. We then
propose a mechanical bending deformation to enable the direct measurement of
the flexoelectric response by eliminating the piezoelectric contribution to
the polarization. Our results show that the bending flexoelectric constants of
Janus TMDCs are significantly higher than that of traditional TMDCs such as
MoS2. The flexoelectric enhancement is found to emerge from the structural
asymmetry that is intrinsic to Janus TMDCs, which results in both stronger
$\sigma-\sigma$ and $\pi-\sigma$ interactions than is found in traditional
TMDCs.

##  II Simulation method

We introduce in this section the computational model we use to calculate the
electrical polarization due to mechanical deformation. Specifically, a
combination of short-range bonded interactions with long-range charge-dipole
(CD) interactions were considered for calculating the forces acting on a given
atomic system. According to the CD model [49; 50], each atom is assumed to
carry a charge $q$ and dipole moment $\mathbf{p}$. The short-range atomic
interactions are modeled using a Stillinger-Weber potential [51], which was
previously used to study the misfit strain induced buckling for lateral
heterostructures with different combinations of TMDCs (MoS2-WSe2; MoS2-WTe2;
MoS2-MoSe2 and MoS2-MoTe2). This potential was shown to capture the previously
reported spontaneous curling behavior of Janus TMDCs [52]. The CD model
requires a parameter $R$ (related to atomic polarizability $\alpha$) to
evaluate the charge and dipole for each atom, which is obtained by matching
the calculated polarizability $(\alpha_{\text{Cal.}})$ with DFT calculated
value $(\alpha_{\text{DFT}})$, where the complete details about the process
for calculating the CD parameters are given in Ref. [31]. Table 1 lists the
lattice parameters of each Janus TMDC, the parameter $R$,
$\alpha_{\text{DFT}}$ and $\alpha_{\text{Cal.}}$. All the simulations in this
study were conducted in the open-source molecular dynamics code LAMMPS [53].
The GAUSSIAN software [54] was employed to estimate $\alpha_{\text{DFT}}$. The
details of the atomic forces resulting from the CD model can be found in [22]
and references therein. blueThe current DFT calculations do not consider spin-
orbit coupling (SOC). We tabulated the spontaneous dipole moment for Janus
TMDCs from the present CD model $(p^{0}_{\text{Cal.}})$ and recent DFT reports
with and without SOC $(p^{0}_{\text{DFT}})$. While SOC does have some effect
on the results, Table. 1 shows that the comparison between
$p^{0}_{\text{Cal.}}$ and $p^{0}_{\text{DFT}}$ is quite satisfactory. This
shows that the present CD model is sufficiently accurate to predict the
flexoelectric properties.

Table 1: blueThe spontaneous dipole moment from the present CD model $(p^{0}_{\text{Cal.}}~{}\text{in e\AA})$ and recent DFT reports with and without SOC $(p^{0}_{\text{DFT}}~{}\text{in e\AA})$. Total polarizability estimated by DFT ($\alpha_{\text{DFT}}$ in $\text{\AA}^{3}$) and the present CD model ($\alpha_{\text{Cal.}}$ in $\text{\AA}^{3}$) [31]. Atomic polarizability ($R_{i}$ in Å) and lattice parameters (a,b are the lattice constants, $l_{1}$ and $l_{2}$ are the bond lengths for M-X and M-Y in Å, respectively) for the Janus TMDCs. For simplicity, in this study the Janus TMDCs are denoted as MXY with M = Mo and W; X,Y = S,Se and Te; X represents the chalcogenide atom with smaller atomic number while Y represents the chalcogenide atom with larger atomic number. Material |  blue$p^{0}_{\text{Cal.}}$ |  blue$p^{0}_{\text{DFT}}$ | $\alpha_{\text{DFT}}$ | $\alpha_{\text{Cal.}}$ | $R_{M}$ | $R_{X}$ | $R_{Y}$ | a,b | $l_{1}$ | $l_{2}$
---|---|---|---|---|---|---|---|---|---|---
MoS2 | 0.0 | 0.0111Reference with SOC [55] | 12.320 | 12.350 | 0.69 | 1.04 | 1.04 | 3.160222Reference [56] | 2.4202 | 2.4202
MoSSe | 0.032 | 0.0391333Reference with SOC [57],444Reference without SOC [34], 0.052555Reference with SOC [58],666Reference without SOC [59] | 13.456 | 13.450 | 0.84 | 1.14 | 1.06 | 3.2882 | 2.4162 | 2.5302
MoSTe | 0.0432 | 0.04123 | 15.590 | 15.593 | 1.0 | 1.16 | 1.02 | 3.3432 | 2.4322 | 2.7152
MoSeTe | 0.0649 | 0.0793 | 16.768 | 16.769 | 1.04 | 1.12 | 1.05 | 3.4122 | 2.5522 | 2.7172
WSSe | 0.0479 | 0.03623, 0.055 | 16.824 | 16.822 | 0.94 | 1.18 | 1.08 | 3.2322 | 2.4212 | 2.5382
WSTe | 0.0574 | 0.03893 | 19.489 | 19.489 | 1.08 | 1.18 | 1.08 | 3.3442 | 2.4382 | 2.7202
WSeTe | 0.0587 | 0.07573 | 21.225 | 21.224 | 1.08 | 1.24 | 1.08 | 3.4132 | 2.5592 | 2.7222

We first validated the CD model by comparing to previously reported in-plane
piezo coefficients [38]. In this study, the in-plane piezoelectric
coefficients of the Janus TMDCs (MXY, where M = Mo, W; X, Y = S, Se, and Te,
where the atomic mass of X is smaller than the atomic mass of Y) are obtained
by subjecting a square Janus TMDC sample with dimensions 80 Å$\times$ 80 Å to
tensile loading in the y-direction as shown in Fig. 1(a). The initial
configuration is flat, and the relaxed atomic configurations of the Janus TMDC
are obtained through energy minimization, after which the atomic charge
$q_{i}$ and dipole moments $\mathbf{p}_{i}$ are derived from the CD model for
atom $i$. The relaxed configuration is bent due to spontaneous curling [52;
32], which arises from the structural asymmetry between the M-X and M-Y
layers, resulting in stretching in the Y layer and compression in the X layer.

The total polarization $\mathbf{P}$ of the system is calculated by
$\mathbf{P}=\frac{1}{V}\left(\sum_{i=1}^{n}\mathbf{p}_{i}\right)$ where $V$ is
the volume and $n$ is the total number of atoms in the system. The in-plane
strain is calculated by $\epsilon^{yy}=\frac{ly-
ly_{\text{ini}}}{ly_{\text{ini}}}$, where $ly$ and $ly_{\text{ini}}$ are the
deformed length and initial length in the y-direction, respectively. The
loading scheme for calculating the in-plane piezoelectric coefficient and
polarization ($P^{y}$) - strain ($\epsilon^{yy}$) diagrams are presented in
Figs. 1(a) and (b). The polarization ($P^{y}$) - strain ($\epsilon^{yy}$)
diagrams are shifted to have zero initial polarization by subtracting the
polarization caused by the initial spontaneous curling [52]. A linear relation
is observed between the polarization $P^{y}$ and given strain $\epsilon^{yy}$
in Fig. 1(b). The slope of the linear relation yields the in-plane
piezoelectric coefficient of the Janus TMDCs. The calculated piezoelectric
coefficients $(d_{yyy}~{}\text{or}~{}d_{11})$ are in good agreement with the
reported DFT values [38] (see Table 2), which validates the effectiveness of
the CD model.

Figure 1: (a) Loading scheme for calculating the in-plane piezoelectric coefficient $d_{11}$ and (b) polarization $P^{y}$ (C/m) vs strain $\epsilon^{yy}$ for the Janus TMDC systems. (c) Schematic of applied bending deformation to calculate the flexoelectric constants. (d) Polarization $P^{z}$ (C/m2) versus strain gradient $K_{eff}~{}(\text{nm}^{-1})$ for MoS2, MoSSe, MoSTe and WSTe. Markers indicate the data from CD model and soild lines indicate the linear fitting. Table 2: In-plane piezoelectric coefficients $d_{yyy}$ ( $\times 10^{-10}$ C/m) for Janus TMDCs using proposed CD model in comparison to previous DFT results ($d_{yyy}^{\text{DFT}}$). $\mu_{zyzy}$ (nC/m) are the bending flexoelectric constants while $l_{2}-l_{1}\>(\text{\AA})$ is the initial asymmetry for Janus TMDCs. Material | $d_{yyy}$ | $d_{yyy}^{\text{DFT}}$ | $l_{2}-l_{1}$ | $\mu_{zyzy}$
---|---|---|---|---
MoS2 | 3.95 | 3.56777Reference [38] | 0.0 | 0.032
MoSSe | 4.099 | 3.747 | 0.114 | 0.117
MoSeTe | 4.676 | 4.357 | 0.165 | 0.120
MoSTe | 4.733 | 4.537 | 0.283 | 0.125
WSSe | 3.144 | 2.577 | 0.117 | 0.089
WSeTe | 3.209 | 3.347 | 0.163 | 0.092
WSTe | 3.327 | 3.487 | 0.282 | 0.114

##  III Results and discussion

In this section, we describe the bending scheme used to study the
flexoelectric properties of Janus TMDCs. To determine the bending
flexoelectric coefficients for the Janus TMDCs, the loading scheme illustrated
in Fig. 1(c) is applied. Specifically, the following displacement field is
applied to the atomic system

$u^{z}=K\frac{y^{2}}{2},$ (1)

where $y$ represents the atom coordinate in the $y$ direction, and
$\frac{1}{2}K$ represents the given strain gradient of the bending plane. We
note that the imposed displacement field in Eq. (1) is imposed starting from
the relaxed, or spontaneously curved configuration in Fig. 1(c). Once the
deformation is prescribed, the boundary region atoms are held fixed while the
interior atoms are allowed to relax to energy minimizing positions using the
conjugate-gradient algorithm, after which the point charges $q_{i}$ and dipole
moments $\mathbf{p}_{i}$ are found for each atom by CD model. The bending
flexoelectric constant $\mu_{zyzy}$ can then be obtained as the slope of the
resulting plot between polarization and strain gradient.

In general, an imposed deformation would induce polarization, which has
contributions from both piezoelectric and flexoelectric effects. The current
bending deformation in Eq. (1) induces only the strain component
$\epsilon^{yz}$ and the strain gradient term
$(\frac{\partial\epsilon^{yz}}{\partial y})$, while the remaining components
of the strain and strain gradient tensors are zero. Therefore, the total
induced polarization in the z-direction is

$P^{z}=d_{zyz}\epsilon^{yz}+\mu_{zyzy}\frac{\partial\epsilon^{yz}}{\partial
y},$ (2)

where $d_{zyz}$ and $\mu_{zyzy}$ are the bending piezoelectric coefficient and
flexoelectric coefficient, respectively. The inset of Fig. 2 represents the
atomic configuration of MoSSe obtained using the OVITO software [60] when the
applied curvature $K=0.05~{}\text{nm}^{-1}$. The atoms are colored according
to the $\epsilon^{yz}$ component of strain, which is calculated from Eq. (3)
as [61; 62]

$\epsilon_{i}^{\alpha\beta}=\frac{1}{2}\left[F_{i}^{\beta\alpha}F_{i}^{\alpha\beta}-\delta^{\alpha\beta}\right],$
(3)

where $\epsilon_{i}$ is the atomic strain for atom $i$, $F_{i}^{\alpha\beta}$
is the $i$th atomic component of deformation gradient, $\delta$ is the
Kronecker delta, $\alpha$ and $\beta$ are the coordinate components.

As shown in Fig. 2, the strain $\epsilon^{yz}$ varies from -0.018 to 0.018
along the $y-$direction in the MoSSe sheet. These strain values represent
averaged values of the atomic strain, which were found by averaging over 24
equal-width bins along the y-direction. The linear variation of atomic strain
indicates that total strain is zero and maintains the point group symmetry
posed by MXY [38]. Overall, the observed symmetry removes the piezoelectric
part of polarization $d_{zyz}\epsilon^{yz}=0$, which supports the assumption
that the imposed bending deformation removes the piezoelectric contribution to
the polarization.

Figure 2: Strain profile $\epsilon_{yz}$ versus y coordinates $Y$. Inset shows
a contour plot of the strain $\epsilon_{yz}$.

The slope of the plot in Fig. 2 between $\epsilon^{yz}$ and the $y-$coordinate
of each bin leads to the effective strain gradient
$K_{eff}=\partial\epsilon^{yz}/\partial y$. The current value of $K_{eff}$
differs from $\frac{1}{2}K$ (from Eq. (1)) by about 15%. Theoretically, the
numerical value of $K_{eff}$ should be equal to $\frac{1}{2}K$ under the
imposed bending deformation. For example, previous works found that for
conventional TMDCs (MoS2, WS2 and CrS2), the effective strain gradient is
equal to half of the given value of $K$ [31] under the same bending
deformation. However, the observed difference in Janus TMDCs is due to the
spontaneous curling effect [52; 32], which arises from the structural
asymmetry between the M-X and M-Y layers, resulting in stretching in the Y
layer and compression in the X layer. In order to account for the spontaneous
deformation, the effective strain gradient $K_{eff}$ is used in our
calculations of the flexoelectric coefficient.

We note that a similar type of curling, i.e. ripplocations, have been observed
in MoS2 by introducing line defects (the removal of sulfur atoms along a line)
[63], which generate local strain gradients and thus flexoelectricity.
However, in the present work we restrict ourselves to calculations of
intrinsic flexoelectricity for spontaneously curled Janus TMDCs.

###  III.1 Flexoelectric effect in MoS2 and MoSSe

In this section, we examine MoS2 and MoSSe to illustrate the effects of
structural asymmetry on the resulting flexoelectric properties. MoS2 is chosen
as a TMDC that does not have structural asymmetry (MXX), while MoSSe is chosen
as a representative Janus TMDC which does have structural asymmetry (MXY). To
aid in the analysis of the resulting flexoelectric constants, we note that the
dipole moment of an atom $\mathbf{p}$ depends primarily on three factors: the
effective atomic polarizability ($R$), the charge induced electric field
$(E^{q-z})$, and the dipole induced electric field $(E^{p-z})$. We focus, as
shown previously in Eq. (2), on the the out-of-plane (z-direction) dipole
moment $p^{z}$ and associated polarization $P^{z}$.

Figure 3: Atomic configurations of (a) undeformed MoS2, (b) deformed MoS2 at
$K_{eff}$ = 0.05 nm-1, (c) undeformed MoSSe and (d) deformed MoSSe at
$K_{eff}$ = 0.05 nm-1.

The bond length between atoms A-B and A-C in Fig. 3(a) is equal to
$2.42~{}\text{\AA}$ $(l_{Mo-S})$ for unbent MoS2. Because the A-B and A-C bond
lengths are the same, there is no initial structural asymmetry for MoS2, and
the bonds induce equal and opposite electric fields, which cause the total
dipole moment, and thus total polarization of each MoS2 unit cell to be zero.
However, significant changes in the bond lengths are observed in the deformed
state (Fig. 3(b)) for a curvature of $K_{eff}$ = 0.05 nm-1. Bond A-B is
stretched from 2.42 to 2.56 Å while bond A-C is compressed from 2.42 to 2.23
Å. This difference in bond lengths breaks the electric field symmetry and
increases the total dipole moment. At $K_{eff}$ = 0.05 nm-1, the total
electric field difference $\Delta(E^{p-z}+E^{q-z})$ with respect to the
initial (undeformed) configuration is 53.06 V/Å, which increases the total
polarization from $0$ to 0.0016 C/m2. The changes to $E^{q-z}$ and $E^{p-z}$
at $K_{eff}$ = 0.05 nm-1 are 45.6 and 7.4 V/Å, respectively. The significant
contribution of $E^{q-z}$ implies the increasing importance of $\pi-\sigma$
interactions in generating the dipole moment $p^{z}$. The $\pi-\sigma$
interactions originate from the coupling between valence electrons and bonding
electrons [64], which are also interpreted as pyramidalization [65; 29; 66].
Furthermore, a recent DFT study reported [67] electron transfer from the $p$
orbitals of S atoms to the $d_{z^{2}}$ orbital of Mo atoms. This electron
transfer modifies the charges on the Mo and S atoms, which generates local
charge induced electric fields, which are captured through $E^{q-z}$ within
the current CD model. Further changes in the bond length between Mo and S
atoms enhance the $\pi-\sigma$ coupling [36], resulting in large $E^{q-z}$ and
$P^{z}$. The variation of total polarization with the effective strain
gradient for MoS2 is given in Fig. 1(d), where the slope represents the
flexoelectric constant for MoS2 $(\mu_{\text{MoS}_{2}})$, which is found to be
$0.032$ nC/m.

In contrast to MoS2, the Janus TMDC MoSSe has a structural asymmetry in the
undeformed configuration between the Mo-S and Mo-Se atomic layers.
Specifically, the A-B $(l_{Mo-Se})$ bond length in Fig. 3(c) is equal to
$2.53~{}\text{\AA}$, while the A-C bond length $(l_{Mo-S})$ is
$2.416~{}\text{\AA}$. The bending deformation of MoSSe to a curvature of
$K_{eff}$ = 0.05 nm-1 stretches $l_{Mo-Se}(l_{2})$ from 2.53 to
$2.72~{}\text{\AA}$ and shrinks $l_{Mo-S}(l_{1})$ from 2.416 to
$2.20~{}\text{\AA}$. The initial bond length asymmetry in MoSSe is further
increased due to the given deformation and helps to produce larger dipole
moments $p^{z}$, compared to MoS2. The dipole moment $p^{z}$ is also related
to the polarizability of the atomic system, where the polarizability of MoSSe
$(\alpha_{\text{MoSSe}})$ is 1.09 times larger than that of MoS2
$(\alpha_{\text{MoS}_{2}})$ as shown in Table. 1, which shows the important
effect of structural asymmetry for the Janus TMDCs.

In addition, the initial bond length asymmetry in unbent MoSSe induces an
intrinsic electric field, which is not present in MoS2, as shown in recent DFT
simulations [34]. That work also reported the non-overlapping of out-of-plane
wavefunctions for electrons and holes due to this electric field, which
implies a weak bonding between the electron-hole pair, and which reduces the
bandgap by pushing the $d$ orbitals of the metal atom closer to the Fermi
level [68]. The $d$ orbital shifting may enhance the charge transfer process
through $\pi-\sigma$ coupling, which represents an easy transfer of charges
from the S or Se atom to the Mo atom. This phenomenon is reflected in our CD
model as the calculated charge on Mo (atom A in Fig. 3(d)) in MoSSe is
$0.733e$, which is significantly larger than for Mo in MoS2 (atom A in Fig.
3(b)), which is $0.282e$, and as such the charge acquired by the Mo atom in
MoSSe is $2.6$ times higher than in MoS2.

For a curvature of $K_{eff}$=0.05 nm-1, the charge transfer-induced change of
$E^{q-z}$ in MoSSe is 117.6 V/Å, which is exactly $2.6$ times higher than the
field induced in MoS2, and which reflects stronger $\pi-\sigma$ coupling in
MoSSe. Furthermore, the value of $E^{p-z}$ is higher in MoSSe (38.4 V/Å) than
in MoS2 (7.4 V/Å), which represents a stronger dipole interaction
($\sigma-\sigma$ coupling) in MoSSe than MoS2. Recent studies on the
electronic properties of strained Janus TMDCs [36] suggest an increased
coupling between the $p$ orbitals of S/Se atoms with the in-plane bonding
orbitals $d_{x^{2}-y^{2}}$ and $d_{xy}$ of Mo atom as a function of changes in
bond angle. In MoSSe, the angle B-A-C in Fig. 3 varies from $81.1^{\circ}$ to
$74.2^{\circ}$ between the initial and deformed states ($K_{eff}$ = 0.05
nm-1), while for MoS2, angle variations are $81.93^{\circ}$ to $78.66^{\circ}$
for the initial and deformed states. The reduction in bond angles between
MoSSe and MoS2 may increase the contribution of $E^{p-z}$ in MoSSe.

Overall, $E^{q-z}$ is higher than $E^{p-z}$ in MoSSe, which represents that
$\pi-\sigma$ coupling is dominant over $\sigma-\sigma$ coupling. DFT
simulations have shown that coupling between $d_{z^{2}}$ and $p$ orbitals
($\pi-\sigma$ coupling) is stronger than the coupling between $p$ and
$d_{x^{2}-y^{2}}$ and $d_{xy}$ ($\sigma-\sigma$ coupling) orbitals [36; 68] in
Janus TMDCs. The total electric field increment $\Delta(E^{p-z}+E^{q-z})$ in
MoSSe is 156.07 V/Å at a curvature of $K_{eff}$ = 0.05 nm-1, which is $2.94$
times higher than in MoS2 at the same curvature. The increased polarizability
and electric field increases the total dipole moment of MoSSe to $3.24$ times
higher than in MoS2. The flexoelectric constant for MoSSe
$(\mu_{\text{MoSSe}})$ is found to be $0.117$ nC/m, which is $3.6$ times
higher than flexoelectric coefficient of MoS2, which is consistent with the
larger polarizability of MoSSe. The polarizability is directly related to the
dielectric constant ($\varepsilon$) of the material [69]. For instance,
$\alpha$ of MoSSe and MoS2 are 13.45 and 12.35 $\text{\AA}^{3}$, respectively
as shown in Table 1. The DFT calculated values for $\varepsilon$ for these
materials are 8.67 and 8.05 [34]. Thus, the increased flexoelectric constant
for MoSSe over MoS2 is in agreement with the fact that felxeoelectric effect
scales with the material dielectric constant [70; 71].

As previously shown in Table 2, the in-plane piezoelectric coefficient for
MoSSe is comparable with that of MoS2. However, the out-of-plane flexoelectric
coefficient is $3.6$ times higher than MoS2 (Table 2), which implies that
there are relative benefits to bending flexoelectricity as compared to in-
plane piezoelectricity when comparing Janus TMDCs to standard TMDCs. In order
to understand this further, we compare the change in total electric field and
asymmetry in bond length between in-plane tensile deformation and out-of-plane
bending deformation for MoSSe for the same strain energy density, which is
defined as the sum of atomic stress times the atomic strain over the volume of
the deformed system. The change in $E^{q-z}$ and $E^{p-z}$ for bending
deformation were previously noted as 117.6 and 38.4 V/Å for a curvature of
$K_{eff}$ = 0.05 nm-1, while for tensile deformation the values are 11.9 and
18.7 V/Å, respectively. This shows that the induced electric fields are higher
in bending than in tension, which generates high dipole moments. It is also
observed that the difference in bond length ($l_{2}-l_{1}$) for bending is
higher than in tension, for the same unit cell compared in tension and
bending. This reflects the larger bond length asymmetry that is induced in
bending, which supports the charge transfer and enhanced $\pi-\sigma$ coupling
based electric fields.

To further compare the resulting electromechanical coupling between standard
TMDCs and Janus TMDCs, we have computed the electrical energy density for
atomic configurations at the same strain energy density for MoSSe and MoS2.
The electrical energy density is defined as the sum of the dot product between
induced polarization and electric field over the volume. The atomic
configuration of MoSSe at $K_{eff}$=0.05 nm-1 gives an electrical energy
density of $2.6~{}\times 10^{9}$ J/m3 and a strain energy density of
$1.13~{}\times 10^{9}$ J/m3. By selecting an atomic configuration for MoS2
under bending deformation with the same strain energy density results in an
electrical energy density that is about 39% of the MoSSe electrical energy
density. The electrical energy density under tensile deformation of MoSSe is
about 84% of the electrical energy density for MoSSe under bending. These
results show that the higher values of electric fields and large dipole
moments under bending deformation point to the advantage of bending as
compared to stretching in generating strong electromechanical coupling in
Janus TMDCs.

###  III.2 Flexoelectric effect among Janus TMDCs

From the previous section, it is clear that the asymmetry in bond lengths
between layers of MXY induce large dipole moments through the increase in
induced electric fields. Table 2 lists the out-of-plane bending flexoelectric
coefficients of the Janus TMDCs along with the initial bond length difference
$l_{2}-l_{1}$ in Fig. 3, from which a positive correlation is identified for
both the MoXY and WXY Janus TMDCs. We note that the out-of-plane piezoelectric
coefficients also show a similar dependence on bond length asymmetry [38]. We
mechanistically examine this correlation further using the electric fields due
to charge-dipole ($E^{q-z}$) and dipole-dipole ($E^{p-z}$) interactions from
the CD model.

Under the prescribed bending scheme, the increase in the total electric field
$E^{q-z}+E^{p-z}$ increases with strain gradient $K_{eff}$ for every Janus
TMDC. However, the relative contribution from $E^{q-z}$ or $E^{p-z}$ to
$E^{q-z}$+$E^{p-z}$ varies between elements of the Janus TMDCs group. In the
MoXY group for a strain gradient of $K_{eff}$ = 0.1 nm-1, the contribution
from the charge-dipole interaction induced electric field $E^{q-z}$ in MoS2,
MoSSe, MoSeTe, and MoSTe to the increase in the total electric field
($E^{q-z}$+$E^{p-z}$) is 84.76%, 75.48%, 70.33% and 56.29%, respectively. The
contribution from the dipole-dipole interaction induced electric field
$E^{p-z}$ to the increase in the total electric field ($E^{q-z}$+$E^{p-z}$) is
then 15.24%, 24.52%, 29.67% and 43.71% in MoS2, MoSSe, MoSeTe, and MoSTe,
respectively.

The bending deformation further develops a coupling among the induced dipole
moments via the $\sigma-\sigma$ interactions, which raises the contribution of
$E^{p-z}$. For example in MoSTe, at $K_{eff}$ = 0.01 nm-1, $E^{q-z}$ and
$E^{p-z}$ are 56.82 (V/Å) and 33.99 (V/Å), respectively. When $K_{eff}$ is
increased to 0.05 nm-1, these values increase to 280.5 (V/Å) and 216.99 (V/Å).
The enhanced dipole-dipole interaction $E^{p-z}$ is due to the reduction in
bond angle via the increased bond length asymmetry due to bending. The
decrease in bond angle (B-A-C in Fig. 3) between initial and deformed states
($K_{eff}$ = 0.05 nm-1) is $6.9^{\circ}$ for MoSSe and $7.4^{\circ}$ for
MoSTe. The difference in $E^{p-z}$ for atom A in MoSSe is $0.38$ V/Å and for
atom A in MoSTe is $0.53$ V/Å (not shown in Fig. 3). This further confirms
that the increased reduction in the bond angle helps in increasing the dipolar
interactions ($\sigma-\sigma$ coupling), which are measured in the form of
$E^{p-z}$. Note that atoms A, B, and C are selected at exactly the same unit
cell locations in MoSSe and MoSTe. The cumulative effect of bond angle
reduction within the unit cell and across the unit cell in the given atomic
system makes the contribution of $E^{p-z}$ significant. The total effect of
increased electric fields and polarizability helps in inducing high
polarization and thus high flexoelectric coefficient for MoSTe over other
Janus TMDCs. A similar trend is found in the WXY group in Table 2, where high
charge transfer and higher bond angle reduction in WSTe compared to other
elements in WXY group. This makes the flexoelectric coefficient of WSTe higher
than WSSe and WSeTe.

Table. 2 also shows that the MoXY group exhibits higher out-of-plane bending
flexoelectric coefficients than the WXY group. For example, MoSTe and WSTe
have the highest flexoelectric coefficients in their groups, where the
flexoelectric coefficient of MoSTe is about 1.09 times larger than WSTe. At
$K_{eff}$ = 0.05 nm-1, the total electric field change in MoSTe is found to be
1.4 times larger than that of WSTe while the polarizability of MoSTe is 0.79
times that of WSTe. Though the polarizability of MoSTe is lower than WSTe, the
larger induced electric fields results in the flexoelectric constant of MoSTe
being slightly larger than that of WSTe.

###  III.3 Bending against the initial spontaneous curvature

Up to this point, we have studied the bending-induced flexoelectric response
of MXY materials while applying the out-of-plane bending toward to their
initial curling direction (towards S in the case of MoSTe). In this section,
we apply bending deformation against the initial curling direction (towards Te
in the case of MoSTe) to study the effect of bending direction on the induced
polarization. We followed the same simulation procedure described in Sec II
except for inserting a minus sign in Eq. (1) to indicate that the applied
deformation is opposite to the initial curling direction.

Fig. 4 shows the variation of $P^{z}$ with $K_{eff}$. The dashed line at
$K_{eff}=K_{flat}$ represents the flattening of the MoSTe sheet from the
initial spontaneously curved state. The notable feature of Fig. 4 is that the
total polarization $P^{z}$ increases up to $K_{flat}$ and decreases
afterwards. This is due to an increase in bond length between Mo and S atoms
which reduces the charge transfer when $K_{eff}$ is larger than $K_{flat}$,
which is seen in the form of decreasing $E^{q-z}$. Specifically, the charge on
an Mo atom at $K_{flat}$ and $1.4\times K_{flat}$ are $0.55e$ and $0.43e$,
respectively, which leads to a reduction in $E^{q-z}$. A simultaneous decrease
in $E^{p-z}$ is also observed due to an increase in the angle B-A-C, which
increases from about $3.96^{\circ}$ from $K_{flat}$ to $1.4\times K_{flat}$.
While the slope of the polarization $P^{z}$ variation with $K_{eff}$ is
similar to the flexoelectric coefficient obtained in the case of bending
towards the spontaneous curling direction, more energy is required to deform
towards the Te layer in MoSTe as compared to the S layer in MoSTe. As a
result, bending towards the S layer in MoSTe is a better choice for energy
conversion. A similar observation is found in the case of WSTe.

Figure 4: The variation of polarization $P^{z}$ with the effective strain
gradient $K_{eff}$ for MoSTe and WSTe when bending against their initial
spontaneous curvature.

##  IV Conclusion

We have utilized classical molecular dynamics coupled with a charge-dipole
model to study the bending flexoelectric response of Janus TMDCs. We further
employed a prescribed bending deformation that enabled us to directly
calculate the flexoelectric constants while eliminating the piezoelectric
contribution to the polarization. In doing so, we found that Janus TMDCs have
flexoelectric constants that are several times larger than traditional TMDCs
such as MoS2. The mechanism underlying this was found to be bond length
asymmetry for the Janus TMDCs between the M-X and M-Y atoms. This bond length
asymmetry was found to lead to stronger $\sigma-\sigma$ interactions with
increasing initial asymmetry, along with stronger $\pi-\sigma$ interactions
due to increased charge transfer, which combine to result in increased
polarization for Janus TMDCs.

The present results also demonstrate the enhanced electromechanical coupling
that results from out-of-plane bending as compared to in-plane stretching for
Janus TMDCs. Specifically, by comparing equivalent states of strain energy
density, we found that the electrical energy density is higher for Janus TMDCs
as compared to MoS2, while the electrical energy density under tension is
smaller than that of bending for Janus TMDCs. These facts imply that there may
be benefits in electromechanical energy conversion for Janus TMDCs by
utilizing out-of-plane bending flexoelectricity as compared to in-plane
piezoelectricity.

## Acknowledgement

The authors gratefully acknowledge the sponsorship from the ERC Starting Grant
COTOFLEXI (No. 802205).

## References

  * Pop  _et al._ [2012] E. Pop, V. Varshney, and A. K. Roy, MRS bulletin 37, 1273 (2012).
  * Ghosh  _et al._ [2010] S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau, and A. A. Balandin, Nature materials 9, 555 (2010).
  * Fei and Yang [2014] R. Fei and L. Yang, Nano letters 14, 2884 (2014).
  * Yu  _et al._ [2017] S. Yu, X. Wu, Y. Wang, X. Guo, and L. Tong, Advanced Materials 29, 1606128 (2017).
  * Wang  _et al._ [2012] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Nature nanotechnology 7, 699 (2012).
  * Hanakata  _et al._ [2016] P. Z. Hanakata, A. Carvalho, D. K. Campbell, and H. S. Park, Physical Review B 94, 035304 (2016).
  * Andrew  _et al._ [2012] R. C. Andrew, R. E. Mapasha, A. M. Ukpong, and N. Chetty, Physical review B 85, 125428 (2012).
  * Javvaji  _et al._ [2016] B. Javvaji, P. Budarapu, V. Sutrakar, D. R. Mahapatra, M. Paggi, G. Zi, and T. Rabczuk, Computational Materials Science 125, 319 (2016).
  * Kim  _et al._ [2015] S. J. Kim, K. Choi, B. Lee, Y. Kim, and B. H. Hong, Annual Review of Materials Research 45, 63 (2015).
  * Sun  _et al._ [2016] Z. Sun, A. Martinez, and F. Wang, Nature Photonics 10, 227 (2016).
  * Mas-Ballesté  _et al._ [2011] R. Mas-Ballesté, C. Gómez-Navarro, J. Gómez-Herrero, and F. Zamora, Nanoscale 3, 20 (2011).
  * Wang  _et al._ [2015] X. Wang, H. Tian, W. Xie, Y. Shu, W.-T. Mi, M. A. Mohammad, Q.-Y. Xie, Y. Wang, J.-B. Xu, and T.-L. Ren, NPG Asia Materials 7, e154 (2015).
  * Song  _et al._ [2017] X. Song, F. Hui, T. Knobloch, B. Wang, Z. Fan, T. Grasser, X. Jing, Y. Shi, and M. Lanza, Appl. Phys. Lett 111, 083107 (2017).
  * Zheng  _et al._ [2017] X.-Q. Zheng, J. Lee, and P. X.-L. Feng, Microsystems and Nanoengineering 3, 17038 (2017).
  * Zelisko  _et al._ [2014] M. Zelisko, Y. Hanlumyuang, S. Yang, Y. Liu, C. Lei, J. Li, P. M. Ajayan, and P. Sharma, Nature Communications 5, 4284 (2014).
  * Blonsky  _et al._ [2015] M. N. Blonsky, H. L. Zhuang, A. K. Singh, and R. G. Hennig, ACS nano 9, 9885 (2015).
  * Kundalwal  _et al._ [2017a] S. I. Kundalwal, S. A. Meguid, and G. J. Weng, Carbon 117, 462 (2017a).
  * Chandratre and Sharma [2012] S. Chandratre and P. Sharma, Applied Physics Letters 100, 023114 (2012).
  * Duerloo and Reed [2013] K.-A. N. Duerloo and E. J. Reed, Nano letters 13, 1681 (2013).
  * Zhang [2017] J. Zhang, Nano Energy 41, 460 (2017).
  * Zhou  _et al._ [2016] Y. Zhou, W. Liu, X. Huang, A. Zhang, Y. Zhang, and Z. L. Wang, Nano Research 9, 800 (2016).
  * Javvaji  _et al._ [2018] B. Javvaji, B. He, and X. Zhuang, Nanotechnology 29, 225702 (2018).
  * Wu  _et al._ [2014] W. Wu, L. Wang, Y. Li, F. Zhang, L. Lin, S. Niu, D. Chenet, X. Zhang, Y. Hao, T. F. Heinz, _et al._ , Nature 514, 470 (2014).
  * Duerloo  _et al._ [2012] K.-A. N. Duerloo, M. T. Ong, and E. J. Reed, J. Phys. Chem. Lett. 3, 2871 (2012).
  * Brennan  _et al._ [2017] C. J. Brennan, R. Ghosh, K. Koul, S. K. Banerjee, N. Lu, and E. T. Yu, Nano Lett, 17, 5464 (2017).
  * Majdoub  _et al._ [2008] M. Majdoub, P. Sharma, and T. Cagin, Physical Review B 77, 125424 (2008).
  * He  _et al._ [2018] B. He, B. Javvaji, and X. Zhuang, Physica B: Condensed Matter 545, 527 (2018).
  * Kalinin and Meunier [2008] S. V. Kalinin and V. Meunier, Physical Review B 77, 033403 (2008).
  * Dumitrică  _et al._ [2002] T. Dumitrică, C. M. Landis, and B. I. Yakobson, Chemical physics letters 360, 182 (2002).
  * Kundalwal  _et al._ [2017b] S. I. Kundalwal, S. A. Meguid, and G. J. Weng, Carbon 117, 462 (2017b).
  * Zhuang  _et al._ [2019] X. Zhuang, B. He, B. Javvaji, and H. S. Park, Physical Review B 99, 054105 (2019).
  * Wang  _et al._ [2018] J. Wang, H. Shu, T. Zhao, P. Liang, N. Wang, D. Cao, and X. Chen, Physical Chemistry Chemical Physics 20, 18571 (2018).
  * Guo [2018] S. D. Guo, Physical Chemistry Chemical Physics 20, 7236 (2018).
  * Jin  _et al._ [2018a] H. Jin, T. Wang, Z. R. Gong, C. Long, and Y. Dai, Nanoscale 10, 19310 (2018a).
  * Shi and Wang [2018] W. Shi and Z. Wang, Journal of Physics Condensed Matter 30, 215301 (2018).
  * Huang  _et al._ [2018] Z. Huang, J. Zhong, H. Liu, Y. Wu, C. Tang, C. He, X. Qi, and L. Xue, Journal of Applied Physics 125, 082516 (2018).
  * Li  _et al._ [2018] R. Li, Y. Cheng, and W. Huang, Small 14, 1 (2018).
  * Dong  _et al._ [2017] L. Dong, J. Lou, and V. B. Shenoy, ACS nano 11, 8242 (2017).
  * Lu  _et al._ [2017] A. Y. Lu, H. Zhu, J. Xiao, C. P. Chuu, Y. Han, M. H. Chiu, C. C. Cheng, C. W. Yang, K. H. Wei, Y. Yang, Y. Wang, D. Sokaras, D. Nordlund, P. Yang, D. A. Muller, M. Y. Chou, X. Zhang, and L. J. Li, Nature Nanotechnology 12, 744 (2017).
  * Zhang  _et al._ [2019] C. Zhang, Y. Nie, S. Sanvito, and A. Du, Nano letters 19, 1366 (2019).
  * Hu and Huang [2017] L. Hu and X. Huang, RSC Advances 7, 55034 (2017).
  * Jiang  _et al._ [2013] J. W. Jiang, Z. Qi, H. S. Park, and T. Rabczuk, Nanotechnology 24, 435705 (2013).
  * Jin  _et al._ [2018b] H. J. Jin, W. Y. Yoon, and W. Jo, ACS Applied Materials and Interfaces 10, 1334 (2018b).
  * Kang  _et al._ [2018] S. Kang, S. Jeon, S. Kim, D. Seol, H. Yang, J. Lee, and Y. Kim, ACS Appl. Mater. Interfaces 10, 27424 (2018).
  * Kang  _et al._ [2019] S. Kang, S. Kim, S. Jeon, W. S. Jang, D. Seol, Y. M. Kim, J. Lee, H. Yang, and Y. Kim, Nano Energy 58, 57 (2019).
  * Abdollahi  _et al._ [2019] A. Abdollahi, N. Domingo, I. Arias, and G. Catalan, Nature Communications 10, 1266 (2019).
  * Zhang  _et al._ [2017] S. Zhang, K. Liu, T. Wu, M. Xu, and S. Shen, Journal of Applied Physics 122, 144103 (2017).
  * Shi  _et al._ [2018] W. Shi, Y. Guo, Z. Zhang, and W. Guo, Journal of Physical Chemistry Letters 9, 6841 (2018).
  * Mayer  _et al._ [2006] A. Mayer, P. Lambin, and R. Langlet, Applied physics letters 89, 063117 (2006).
  * Mayer [2007] A. Mayer, Physical Review B 75, 045407 (2007).
  * Jiang [2019] J.-W. Jiang, Acta Mechanica Solida Sinica 32, 17 (2019).
  * Xiong  _et al._ [2018] Q.-l. Xiong, J. Zhou, J. Zhang, T. Kitamura, and Z.-h. Li, Physical Chemistry Chemical Physics 20, 20988 (2018).
  * Plimpton [1995] S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Frisch  _et al._ [2016] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, “Gaussian 16 Revision B.01,” (2016), gaussian Inc. Wallingford CT.
  * Defo  _et al._ [2016] R. K. Defo, S. Fang, S. N. Shirodkar, G. A. Tritsaris, A. Dimoulas, and E. Kaxiras, Physical Review B 94, 155310 (2016).
  * Hu  _et al._ [2018] T. Hu, F. Jia, G. Zhao, J. Wu, A. Stroppa, and W. Ren, Physical Review B 97, 235404 (2018).
  * Xia  _et al._ [2018] C. Xia, W. Xiong, J. Du, T. Wang, Y. Peng, and J. Li, Physical Review B 98, 165424 (2018).
  * Li  _et al._ [2017] F. Li, W. Wei, P. Zhao, B. Huang, and Y. Dai, Journal of Physical Chemistry Letters 8, 5959 (2017).
  * Ji  _et al._ [2018] Y. Ji, M. Yang, H. Lin, T. Hou, L. Wang, Y. Li, and S. T. Lee, Journal of Physical Chemistry C 122, 3123 (2018).
  * Stukowski [2009] A. Stukowski, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2009).
  * Falk and Langer [1997] M. L. Falk and J. S. Langer, Physical Review E 57, 7192 (1997).
  * Shimizu  _et al._ [2007] F. Shimizu, S. Ogata, and J. Li, Materials Transactions 48, 2923 (2007).
  * Fabbri  _et al._ [2016] F. Fabbri, E. Rotunno, E. Cinquanta, D. Campi, E. Bonnini, D. Kaplan, L. Lazzarini, M. Bernasconi, C. Ferrari, M. Longo, G. Nicotra, A. Molle, V. Swaminathan, and G. Salviati, Nature Communications 7, 1 (2016).
  * Gleiter [1987] R. Gleiter, Pure and Applied Chemistry 59, 1585 (1987).
  * Surya  _et al._ [2012] V. J. Surya, K. Iyakutti, H. Mizuseki, and Y. Kawazoe, Computational Materials Science 65, 144 (2012).
  * Nikiforov  _et al._ [2014] I. Nikiforov, E. Dontsova, R. D. James, and T. Dumitrică, Physical Review B 89, 155437 (2014).
  * Pike  _et al._ [2017] N. A. Pike, B. Van Troeye, A. Dewandre, G. Petretto, X. Gonze, G. M. Rignanese, and M. J. Verstraete, Physical Review B 95, 201106(R) (2017).
  * Er  _et al._ [2018] D. Er, H. Ye, N. C. Frey, H. Kumar, J. Lou, and V. B. Shenoy, Nano Letters 18, 3943 (2018).
  * Pan  _et al._ [2016] L. Pan, B. Zou, and L. J. Shi, Physics Letters, Section A: General, Atomic and Solid State Physics 380, 2227 (2016).
  * Yudin and Tagantsev [2013] P. V. Yudin and A. K. Tagantsev, Nanotechnology 24, 432001 (2013).
  * Zubko  _et al._ [2013] P. Zubko, G. Catalan, and A. K. Tagantsev, Annual Review of Materials Research 43, 387 (2013).