ə@ ɛE ɪI ʊU ɑA ʌ2  1Electrical and Computer Engineering Department, University
of British Columbia, Canada
2College of Arts, Media and Design, Northeastern University, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# A comparative study of two-dimensional vocal tract acoustic modeling based
on Finite-Difference Time-Domain methods

###### Abstract

The two-dimensional (2D) numerical approaches for vocal tract (VT) modelling
can afford a better balance between the low computational cost and accurate
rendering of acoustic wave propagation. However, they require a high spatio-
temporal resolution in the numerical scheme for a precise estimation of
acoustic formants at the simulation run-time expense. We have recently
proposed a new VT acoustic modelling technique, known as the 2.5D Finite-
Difference Time-Domain (2.5D FDTD), which extends the existing 2D FDTD
approach by adding tube depth to its acoustic wave solver. In this work,
first, the simulated acoustic outputs of our new model are shown to be
comparable with the 2D FDTD and a realistic 3D FEM VT model at a low spatio-
temporal resolution. Next, a radiation model is developed by including a
circular baffle around the VT as head geometry. The transfer functions of the
radiation model are analyzed using five different vocal tract shapes for vowel
sounds /a/, /e/, /i/, /o/ and /u/.

Keywords: computational acoustics, vocal tract, FDTD, articulatory speech
synthesis

##  1 Introduction

Various articulatory models have been developed to approximate the
sophisticated upper vocal tract geometry and capture its acoustic features
using a physics-based acoustic wave solver. The 3D acoustic analysis [10, 11]
can precisely characterize the vocal tract’s spectral features while offering
better geometrical flexibility. However, such models are not pragmatic for
designing real-time speech synthesizers due to their high computational cost.
In contrast to 3D, the 1D vocal tract models are computationally lightweight.
Nevertheless, their over-simplified representation of a tube structure
substantially eliminates its geometrical details. Hence, the acoustic analysis
of 1D models is not accurate at higher formants.

As an alternative, the 2D acoustic wave solvers [2, 8] represent the complex
vocal tract geometry as a 2D contour (a mid-sagittal cut of a 3D VT geometry)
and capture the wave interaction only along the mid-sagittal plane. Due to the
dimensionality reduction, the 2D models are lightweight compare to 3D ones.
The most well-known 2D acoustic analysis methods for vocal tract modelling are
the finite element method (FEM) and FDTD. Since 2D models do not lump off-
plane wave interaction, direct use of cross-sectional areas from the vocal
tract MRI images to represent 2D contours inside a simulation grid yields
erroneous acoustic output. Alternatively, a nonlinear area function
transformation method needs to be implemented to tune the acoustic features of
2D models that can match the more complex 3D ones. However, this procedure
kills the performance boost of 2D models. Consequently, we have proposed a new
vocal tract acoustic model (2.5D FDTD) that improves the 2D model by capturing
the wave interaction across the mid-sagittal plane. The next section briefly
discusses the numerical implementation of a 3D tube using 2D and 2.5D FDTD
schemes.

##  2 Methodology

###  2.1 Acoustic Wave Solver

We implemented the FDTD numerical analysis technique to discretize the 2D
acoustic wave equation on a staggered Yee grid. The 2D Yee scheme consists of
a rectangular grid, define acoustic parameters (pressure ($p^{(n)}$) and
velocity ($v_{x}^{(n)}$, $v_{y}^{(n)}$)) at each grid point, a squared 2D cell
[1, 12], as shown in Figure 1. However, the employed technique for acoustic
simulation does not facilitate dynamic boundary conditions. Hence, a new
scalar field $0\leq\beta(x,y,t)\leq 1$ was introduced to the solver, which
could transit smoothly between $\beta=1$ (air) and $\beta=0$ (boundary). At
$\beta=0$, a prescribed velocity $v=\bm{v_{b}}$ was enforced to handle the
boundary condition. Specifically, per each time step ($n$), the wave solver
updates pressure component $p^{(n)}$ and velocity components ($v_{x}^{(n)}$,
$v_{y}^{(n)}$) at each grid point, by solving the below linear wave equations
in the time domain [14]:

$\displaystyle\frac{\partial p}{\partial t}+(1-\beta)p\ =\ -\rho c^{2}(\
\frac{\partial v_{x}}{\partial x}\ +\ \frac{\partial v_{y}}{\partial y}\ )$
(1) $\displaystyle\beta\frac{\partial\bm{v}}{\partial t}\ +\
\left(1-\beta\right)\bm{v}\ =\ -\beta^{2}\frac{\bm{\nabla}p}{\rho}\ +\
\left(1-\beta\right)\bm{v}_{b}$ (2)

The 2.5D FDTD follows the above rationale but improves upon it by
incorporating new impedance terms in its 2D acoustic wave solver, known as
tube depth. The tube’s depth $(D)$, i.e., its continuous extension along the
$z$ axis (with $x$ and $y$ axis being the dimension of starting 2D scheme), is
derived from the cross-sectional area of the tube and sampled at each point of
the scheme as shown in Figure 1. Then the resulting depths are mapped to their
respective acoustic parameters ($p$, $\bm{v_{x}}$, $\bm{v_{y}}$) as
demonstrated by [6]. The inclusion of tube depth allows the wave solver to
apprehend wave interactions across the mid-sagittal plane, in turn eliminating
the nonlinear cross-sectional area functions transformation requirement.
Therefore, it is expected that the time complexity of a 2.5D FDTD tube model
will be better than the existing 2D models. Basically, with the help of
following discrete update rules, we adopt a time marching algorithm (see
Algorithm 1) to sample the acoustic parameters at every grid point:

Figure 1: Representation of acoustic parameters ($p$, $\bm{v_{x}}$,
$\bm{v_{y}}$) and tube depth ($\bar{D}$, $D_{(x)}$ and $D_{(y)}$) inside the
2D Yee scheme

$\displaystyle p^{(n+1)}\ =\ \frac{\bar{D}p^{(n)}-\rho c^{2}\Delta t\
\bm{\widetilde{\nabla}}\cdot\bm{V}^{(n)}}{\bar{D}}$ (3)
$\displaystyle\bm{v}^{(n+1)}\ =\ \frac{\beta\bm{v}^{(n)}-\beta^{2}\Delta t\
\widetilde{\nabla}p^{(n+1)}/\rho\ +\ \Delta t(1-\beta)\bm{v_{b}}}{\beta+\Delta
t(1-\beta)}$ (4)

with:

$\displaystyle\bm{V}\ =\ (\ D_{(x)}v_{x},\ D_{(y)}v_{y}\ )$ (5)

Algorithm 1 FDTD Time-Marching algorithm

Input: VT area function $a(x)$, audio time $t$

1: Initialize physical constants: air density $\rho$, sound speed $c$,
boundary admittance coefficient $\mu$.

2: Initialize the simulation sampling rate $R$.

3: Set the temporal resolution ($\Delta t$) and grid resolution ($\Delta s$)
with $R$ and CFL condition.

4: Create a optimized grid having size $(M\times N)$.

5: Set the total number of time steps $T$ of the simulator using $t$ and
$\Delta t$.

6: Define boundary cells with $a(x)$ and normal acoustic impedance $Z$.

7: Set depth values ($\bar{D}$, $D_{(x)}$ and $D_{(y)}$) for each grid cell
from $a(x)$.

8: Define source excitation cells

9: Initialize source excitation velocity $v_{e}$.

10: Initialize acoustic components ($p$, $v_{x}$ and $v_{y}$) for each grid
cell.

11: for $n=1...T$ do

12: for $i=1...M$ do

13: for $j=1...N$ do

14: Update $p^{n+1}(i,j)$ with $v^{n}_{x}(i,j)$, $v^{n}_{y}(i,j)$ and $D(i,j)$
(Eq. 3 and Eq.5)

15: if $(i,j)=$ excitation cell then

16: $v^{n}_{x}\leftarrow v^{n}_{x}+v^{n}_{e}$ and $v^{n}_{y}\leftarrow
v^{n}_{y}+v^{n}_{e}$

17: end if

18: if $(i,j)=$ boundary cell then

19: $v^{n}_{b}\leftarrow p^{n}(i,j)/Z$

20: end if

21: Update $v^{n+1}_{x}(i,j)$ with $p^{n+1}(i,j)$ and $v^{n}_{b}$ (Eq. 4)

22: Update $v^{n+1}_{y}(i,j)$ with $p^{n+1}(i,j)$ and $v^{n}_{b}$ (Eq. 4)

23: end for

24: end for

25: end for

###  2.2 Vocal Tract Wall Losses

To simulate vocal tract wall reflection, we adopted the local reactive
boundary approach as proposed here [13, 10]. The method estimates wall losses
(locally reacting soft walls) and enforces boundary condition by using
$\bm{v_{w}}$, the normal particle velocity going into or coming from a wall.

$\bm{v_{w}}=\frac{p_{w}}{Z_{n}}\bm{\hat{n}}$ (6)

where $p_{w}$ is the current pressure value in an air cell, located in front
of the wall cell. And $Z_{n}$, the normal acoustic impedance can be computed
by using the normal sound absorption coefficient $\alpha_{n}$ as follows,

$Z_{n}=\rho c\frac{1+\sqrt{1-\alpha_{n}}}{1-\sqrt{1-\alpha_{n}}}$ (7)

The boundary absorption coefficient $\alpha_{n}$ can be derived from the
boundary admittance coefficient $\mu\in(0,1)$. During the simulation, the
solver combines equation (6) with equation (4) by setting
$\bm{v_{b}}=\bm{v_{w}}$ for all the velocity vectors which are positioned in
between an air cell and wall cell. However, this method does not comprehend
frequency dependent wall losses.

###  2.3 Excitation Model

For acoustic analysis and speech production, articulatory models require a
source excitation function just above the glottal-end to input acoustic energy
into the tube [5]. In this work, we are only characterizing the acoustic
behaviour of a vocal tract tube through its transfer function analysis. Hence,
a band-passed velocity pulse having frequency range 2Hz-20kHz was injected
near the glottis for various vocal tract shapes, and the corresponding
transfer functions (impulse response) were obtained at the mouth-end.

###  2.4 Radiation Model

Figure 2: (a) Simplified head geometry representation with circular mouth
aperture (b) Snapshot of acoustic wave propagation inside the vocal tract tube
and emanating from mouth aperture at time instant $t=6.1ms$ for vowel /a/.

The outward wave propagation at the mouth-end (lips) contributes to acoustic
energy losses in the vocal tract. It is one of the essential dissipation
mechanisms in voice production. Hence, a nonzero radiation model of a vocal
tract influences the synthesized audio output by lowering the formant
frequencies and increasing their bandwidths. It might be challenging to build
a model that can precisely attain the radiation effects, considering a
realistic human head geometry. However, it has already been demonstrated that
the inclusion of a spherical baffle (a simplified approximation of the actual
head geometry) around the 3D vocal tract model offers promising results for
vowel production [3]. Since the 2D model uses a mid-sagittal cut of a 3D vocal
tract, the spherical baffle can be replaced by a circular baffle for the 2.5D
radiation model as shown in Figure 2(a). Moreover, to illustrate acoustic
waves’ propagation emanating from the mouth toward infinity (Figure 2(b)), the
computational domain can be extended out of the vocal tract by using perfectly
matched layers (PMLs). PMLs help in eliminating spurious wave reflection at
the domain boundaries.

We create a fixed-sized circular baffle with a diameter of 0.20m around the
2.5D vocal tract contour inside the simulation domain. The circular baffle
center has been adjusted to fit in the entire vocal tract contour while
connecting its exit cross-section. Inside the 2.5D scheme, the region enclosed
by the vocal tract and baffle boundary is set to air cells. In order to absorb
outgoing waves and approximate infinite space, $6$ layers of PML are added at
all sides of the domain boundary.

##  3 Experimental Setup

Figure 3: 2.5D vocal tract tube representation for vowel /a/.

A comparative study was carried out for two different scenarios to
characterize the acoustic features of the 2.5D vocal tract model: (1)
numerical simulation with computational domain ended at the mouth aperture by
imposing Dirichlet boundary conditions (open-end boundary condition) and (2)
numerical simulation with free-field radiation, which includes head geometry
and PML layers (radiation condition). These conditions were examined using the
transfer functions analysis approach for five different vocal tract
shapes(vowel sounds: /a/, /e/, /i/, /o/ and /u/). The 2.5D vocal tract shapes
were derived using Story’s area function dataset [9], as it is the standard
one for the vocal tract acoustic analysis. We have only considered the vocal
tract geometry as a straight tube with circular cross-sections for this study
as shown in Figure 3.

For open-end boundary conditions, we use a precise 3D FEM vocal tract model as
the reference for the transfer function analysis. A virtual microphone is
placed 3mm inside the mouth opening to record pressure variations. In the case
of free-field radiation, we compare the transfer functions of the 2.5D model
with Story’s calculated formants, obtained from the frequency-domain analysis
of vocal tracts [7]. And the microphone is positioned 3mm outside of the mouth
to capture the radiating acoustic waves. Courant-Friedrichs-Lewy (CFL)
stability condition in two dimensions: $\Delta t\leq\Delta s/\sqrt{2}c$ is
imposed, where $c$ is the speed of sound.

The simulation generates $50$ms of synthesized audio for every vowel sounds.
During the simulation, we set the following physical parameters fixed: air
density $\rho=1.14$ kg/m3, boundary admittance $\mu=0.005$ and sound speed
$c=350$ m/s. We implement both the 2D and 2.5D vocal tract models in the
MATLAB environment. The simulation runs on a workstation equipped with an
Intel Core i7-8700K processor. The custom code for the 2.5D vocal tract model
is publicly available here111https://github.com/Debasishray19/Talking-Tube.

##  4 Model Validation

We follow the below steps to analyze the acoustic behaviour of the 2.5D model:

  1. (A)

Transfer function analysis with open-end boundary condition for the 2D FDTD
and 2.5D FDTD vocal tract models (vowel sound: /a/, /i/ and /u/).

  2. (B)

Model performance evaluation (simulation run-rime).

  3. (C)

Transfer function analysis for the 2.5D vocal tract model with free-field
radiation (vowel sound: /a/, /e/, /i/. /o/ and /u/).

For Step (A) and (B), we simulate each vowel sound under three different
spatial grid resolution for both 2D and 2.5D FDTD methods: low ($\Delta
s=0.74$mm), mid ($\Delta s=0.28$mm) and high ($\Delta s=0.18$mm). The high
spatial resolution should produce precise acoustic output as it yields vocal
tract geometry with greater details. However, it requires a larger simulation
domain size. The grid resolution values are empirically derived. Step (A) and
(B) examine the improvement in a 2.5D vocal tract model over its 2D simulation
due to the addition of tube depth. Step (C) investigates the change in
transfer function with the inclusion of a circular baffle for the radiation
model.

##  5 Result

We applied Fast Fourier Transformation (FFT) to the recorded pressure waves to
obtain their transfer function. Formants were extracted by considering local
maxima in transfer function curves. Here we are only analyzing the first three
formants since they are mainly responsible for distinguishing vowel sounds
[11]. The formants for the 3D FEM model can be found here [4].

###  5.1 Step A: Transfer Function Analysis

Transfer function analysis for the 2D and 2.5D simulation with open-end
boundary condition222Colormap: For 2D simulation, For 2.5D simulation.,333All
the formant frequencies are in Hz..

Resolution ($\Delta s$) | F1 | F2 | F3
---|---|---|---
low (% Error) |  660(-5.16) 700(0.57) |  1040(-2.61) 1040(-2.61) |  2960(-2.33) 3020(-0.36)
mid (% Error) |  660(-5.16) 680(-2.29) |  1040(-2.61) 1060(-0.74) |  3000(-1.02) 3060(0.95)
high(% Error) |  700(0.57) 680(-2.29) |  1060(-0.74) 1040(-2.62) |  3020(-0.36) 3040(0.29)
Table 1: Positional errors of the first three formants computed for vowel sound /a/, with respect to 3D FEM model. Resolution ($\Delta s$) | F1 | F2 | F3
---|---|---|---
low (% Error) |  240(-8.74) 260(-1.14) |  2120(0.42) 2160(2.32) |  3020(0.33) 3060(1.66)
mid (% Error) |  260(-1.14) 260(-1.14) |  2140(1.37) 2140(1.37) |  3000(-0.33) 3040(0.99)
high(% Error) |  260(-1.14) 260(-1.40) |  2140(1.37) 2160(2.32) |  3020(0.33) 3040(0.99)
Table 2: Positional errors of the first three formants computed for vowel sound /i/, with respect to 3D FEM model. Resolution ($\Delta s$) | F1 | F2 | F3
---|---|---|---
low (% Error) |  260(0.38) 260(0.38) |  700(-7.52) 720(-4.88) |  2260(-0.17) 2300(1.59)
mid (% Error) |  220(-15.0) 260(0.38) |  600(-20.7) 700(-7.52) |  2260(-0.17) 2300(1.59)
high(% Error) |  260(0.38) 260(0.38) |  700(-7.52) 720(-4.88) |  2280(0.70) 2300(1.59)
Table 3: Positional errors of the first three formants computed for vowel
sound /u/, with respect to 3D FEM model.

###  5.2 Step B: Model Performance

The computational cost of FDTD vocal tract model depends upon the spatial grid
resolution and simulation domain size. We implement an algorithm that enforces
an optimized domain size for FDTD models. The table below presents the
simulation domain size and time duration corresponding to each grid resolution
for vowel /u/.

Resolution ($\Delta s$) |  Domain Size (width$\times$height) |  Duration (Approx.) (in seconds)
---|---|---
low ($0.74mm$) |  $266\times 35$ |  280
mid ($0.28mm$) |  $702\times 85$ |  $4.98\text{e+}3$
high($0.18mm$) |  $1051\times 127$ |  $1.85\text{e+}4$
Table 4: Simulation run-time at different grid resolutions.

###  5.3 Step C: Transfer Functions Analysis

The table below presents transfer function analysis for the 2.5D vocal tract
models with the free-field radiation. We have only considered the low grid
resolution ($\Delta s=0.74$mm) for simulation of vocal tract radiation models.

Vowel |  F1 |  F2 |  F3
---|---|---|---
/a/ |  $640$ $-7.78\%$ |  $1020$ $8.28\%$ |  $3000$ $1.76\%$
/e/ |  $340$ $-29.75\%$ |  $1920$ $1.74\%$ |  $2320$ $-3.49\%$
/i/ |  $260$ $-20.00\%$ |  $2080$ $-2.75\%$ |  $2960$ $-0.53\%$
/o/ |  $480$ $-3.80\%$ |  $760$ $-1.04\%$ |  $2340$ $-3.34\%$
/u/ |  $260$ $-17.19\%$ |  $720$ $2.56\%$ |  $2180$ $-5.13\%$
Table 5: Formant frequencies and the percentage positional errors of the first
three formants for the vocal tract radiation model.

##  6 Discussion & Conclusion

The transfer function analysis is a standard practice to validate acoustic
features of area-based articulatory models. In the first step of our
evaluation, the 2.5D FDTD simulation produces positional errors below $2\%$ at
a low grid resolution. Moreover, this characteristic is consistent across all
three vowels. However, the positional error is high for the 2D FDTD simulation
at a low grid resolution in most cases. Table 1, 2 and 3 show that the
acoustic outputs of the 2D FDTD model can be tuned by increasing its
computational grid resolution whereas this approach doesn’t make much
difference for the 2.5D FDTD simulation. Table 4 shows for a real-time speech
synthesizer the low grid resolution is always desirable.

In the case of free-field of radiation, Table 5 presents the first three
formants and their positional errors. The positional error of the first
formant is high across all the vowels. And the second and third formants
produce relatively lesser error values. However, these positional errors are
determined relative to Story’s 1D vocal tract model, an oversimplified
representation of the actual vocal tract. Hence, as future work, the 2.5D
radiation model’s acoustic characteristics need to be compared with a 3D vocal
tract radiation model.

Additionally, our current 2.5D model should only be seen as simple
approximations of the voice radiation mechanism. Nevertheless, the influence
of few head details like the nose and cross-section of the mouth aperture is
worth studying. As well as, our current 2.5D FDTD model does not include
boundary losses for the tube depth. All these current limitations need to be
incorporated for the future development of the 2.5D vocal tract acoustic
model.

##  7 Acknowledgements

This work is supported by the Natural Sciences and Engineering Research
Council (NSERC).

## References

  * [1] Andrew Allen and Nikunj Raghuvanshi “Aerophones in flatland: Interactive wave simulation of wind instruments” In _ACM Transactions on Graphics (TOG)_ 34.4 ACM New York, NY, USA, 2015, pp. 1–11
  * [2] Marc Arnela and Oriol Guasch “Two-dimensional vocal tracts with three-dimensional behavior in the numerical generation of vowels” In _The Journal of the Acoustical Society of America_ 135.1 ASA, 2014, pp. 369–379
  * [3] Marc Arnela, Oriol Guasch and Francesc Alías “Effects of head geometry simplifications on acoustic radiation of vowel sounds based on time-domain finite-element simulations” In _The Journal of the Acoustical Society of America_ 134.4 Acoustical Society of America, 2013, pp. 2946–2954
  * [4] Marc Arnela Coll “Numerical production of vowels and diphthongs using finite element methods”, 2015
  * [5] Debasish Ray Mohapatra and Sidney Fels “Limitations Of Source-Filter Coupling In Phonation” In _arXiv preprint arXiv:1811.07435_ , 2018
  * [6] Debasish Ray Mohapatra, Victor Zappi and Sidney Fels “An extended two-dimensional vocal tract model for fast acoustic simulation of single-axis symmetric three-dimensional tubes” In _arXiv preprint arXiv:1909.09585_ , 2019
  * [7] Man Sondhi and Juergen Schroeter “A hybrid time-frequency domain articulatory speech synthesizer” In _IEEE Transactions on Acoustics, Speech, and Signal Processing_ 35.7 IEEE, 1987, pp. 955–967
  * [8] Matt Speed, Damian Murphy and David M Howard “Characteristics of two-dimensional finite difference techniques for vocal tract analysis and voice synthesis” In _Tenth Annual Conference of the International Speech Communication Association_ , 2009
  * [9] Brad H Story “Comparison of magnetic resonance imaging-based vocal tract area functions obtained from the same speaker in 1994 and 2002” In _The Journal of the Acoustical Society of America_ 123.1 ASA, 2008, pp. 327–335
  * [10] Hironori Takemoto, Parham Mokhtari and Tatsuya Kitamura “Acoustic analysis of the vocal tract during vowel production by finite-difference time-domain method” In _The Journal of the Acoustical Society of America_ 128.6 ASA, 2010, pp. 3724–3738
  * [11] Tomáš Vampola, Jaromír Horáček, Anne-Maria Laukkanen and Jan G Švec “Human vocal tract resonances and the corresponding mode shapes investigated by three-dimensional finite-element modelling based on CT measurement” In _Logopedics Phoniatrics Vocology_ 40.1 Taylor & Francis, 2015, pp. 14–23
  * [12] Kane Yee “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media” In _IEEE Transactions on antennas and propagation_ 14.3 IEEE, 1966, pp. 302–307
  * [13] Takatoshi Yokota, Shinichi Sakamoto and Hideki Tachibana “Visualization of sound propagation and scattering in rooms” In _Acoustical science and technology_ 23.1 ACOUSTICAL SOCIETY OF JAPAN, 2002, pp. 40–46
  * [14] Victor Zappi et al. “Towards real-time two-dimensional wave propagation for articulatory speech synthesis” In _Proceedings of Meetings on Acoustics 171ASA_ 26.1, 2016, pp. 045005 ASA


# Creating a three dimensional intrinsic electric dipole on rotated CrI3
bilayers

Shiva P. Poudel,1 Juan M. Marmolejo-Tejada,2 Joseph E. Roll,1 Martín A.
Mosquera,2 and Salvador Barraza-Lopez1 1\. Department of Physics, University
of Arkansas, Fayetteville, AR 72701, USA and MonArk NSF Quantum Foundry,
University of Arkansas, Fayetteville, AR 72701, USA
2\. Department of Chemistry and Biochemistry, Montana State University,
Bozeman, MT 59717 USA and MonArk NSF Quantum Foundry, Montana State
University, Bozeman, MT 59717 USA

###### Abstract

Two-dimensional (2D) materials are being explored as a novel multiferroic
platform. One of the most studied magnetoelectric multiferroic 2D materials
are antiferromagnetically-coupled (AFM) CrI3 bilayers. Neglecting magnetism,
those bilayers possess a crystalline point of inversion, which is only removed
by the antiparallel spin configuration among its two constituent monolayers.
The resultant intrinsic electric dipole on those bilayers has a magnitude no
larger than 0.04 pC/m, it points out-of-plane, and it reverts direction when
the–Ising-like–cromium spins are flipped (toward opposite layers versus away
from opposite layers). The combined presence of antiferromagnetism and a weak
intrinsic electric dipole makes this material a two-dimensional
magnetoelectric multiferroic. Here, we remove the crystalline center of
inversion of the bilayer by a relative $60^{\circ}$ rotation of its
constituent monolayers. This process enhances the out-of-plane intrinsic
electric dipole tenfold with respect to its magnitude in the non-rotated AFM
bilayer and also creates an even stronger and switchable in-plane intrinsic
electric dipole. The ability to create a three-dimensional electric dipole is
important, because it enhances the magnetoelectric coupling on this
experimentally accessible 2D material, which is explicitly calculated here as
well.

##  I Introduction

Ferroicity is a symmetry-breaking quantum phenomenon [1] that gives rise to
domain physics. Magnetoelectric multiferroics [2, 3] lack spatial and time
inversion symmetries [4, 5, 6, 7, 8, 9, 10, 11, 12], and the experimental
exfoliation of an atomically-thin material (graphene) [13, 14, 15] is now
evolving into a field aiming to deliver advanced functional multiferroic
materials at the atomic-thickness limit.

Effects from coupled magnetic and electric fields have been studied for over a
century: Röntgen recorded magnetic fields when dielectric objects move within
a constant electric field [16], Pierre Curie concluded that magnetoelectric
effects are due to broken symmetries [1], and Debye coined the term
“magnetoelectric directional effect” [17]. The first multiferroic material was
grown in 1959 [18], and the demonstration of intrinsic magnetoelectric
coupling in antiferromagnetically-coupled Cr2O3 was rationalized by
Dzyaloshinskii [19], and experimentally demonstrated by Astrov [5]. The main
crystallographic families for bulk multiferroics [6, 20, 9, 10, 21, 7, 22, 11,
23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36] are perovskites [37,
2], compounds with hexagonal structure [38], boracites [39], and BaMF4
compounds [40, 41].

Initial forays into multiferroics within the 2D materials community are
studies of quantum paraelastic/ferromagnetic SnO monolayers [42, 43, 44], and
of ferroelastic/ferroelectric (experimentally available) group-IV
monochalcogenide monolayers [45, 46, 47, 48, 49]. Magnetoelectric couplings
induced by the proximal placement of ferroelectric and ferromagnetic 2D
materials have been discussed, too [50, 51, 52, 53], along with theoretical
and experimental reports on the change of magnetic ordering by charge doping
[54, 55], and proposals for atomistic defects as sources of multiferroic
behavior [56]. Analytical models exploring magnetoelectric couplings in 2D
materials exist as well [57, 58]. Notably, layered multiferroic NiI2 was
experimentally demonstrated recently [59, 60].

A recent and comprehensive study of ferroelectricity on bilayers based on an
analysis of non-magnetic groups appeared recently [61]. Nevertheless, being
one of the most studied magnetoelectric bilayer systems, the electric control
of magnetization on AFM-coupled CrI3 bilayers [62, 63, 64, 65, 66] (this is,
their magnetoelectric coupling) was experimentally demonstrated back in 2019
through second harmonic generation [67]. The coupling relies on the broken
inversion symmetry induced by the alternating (Ising-like) out-of-plane
magnetization at opposite monolayers [68, 67], which is furnished when
considering the full magnetic group. The magnetic ordering of monolayers [69]
and bilayers [70] can also be tuned by uniaxial strain.

Here, we optimize the intrinsic electric dipole moment $\mathbf{P}$ of a CrI3
bilayer by a relative $60^{\circ}$ rotation among layers (a concept recently
applied to MnSe bilayers [71]). This procedure results on a ten-fold increase
of the out of plane intrinsic electric dipole and it gives rise to an even
larger in-plane intrinsic electric dipole [61]. This effect is unlike what is
reported on MnSe bilayers, which only develop an out-of-plane intrinsic
electric dipole [71]. The three-dimensional intrinsic polarization enhances
the magnetoelectric coupling of CrI3 bilayers, which is calculated here
explicitly.

The manuscript is organized as follows: The numerical methods are described in
Sec. II. The results and discussion are provided in Sec. III. Conclusions can
be found in Sec. IV.

##  II Methods

This study is based on density functional theory with the PBE approximation
for exchange-correlation [72] as implemented in the VASP package [73, 74, 75].
The implementation of spin-orbit-coupling is inconsistent with the use of ab
initio van der Waals corrections on this code, and empirical van der Waals
corrections due to Grimme [76] were utilized for that reason (see Fig. S1 and
Table SI in the Supplementary Material [77] for an extended discussion of
structure versus exchange-correlation functionals). Structural optimizations
were performed using a $15\times 15\times 1$ $k-$point mesh. We used an energy
cutoff of 500 eV for the expansion of electronic states in terms of plane
waves, an energy convergence criteria of $10^{-6}$ eV, and a force convergence
limit of $10^{-4}$ eV/Å. Structural optimizations were performed in
configurations whereby the angle among in-plane lattice vectors
$\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ was kept fixed at 120∘. All calculations
included dipole corrections, and the lattice constant parallel to the two-
dimensional material’s normal was set to $\mathbf{a}_{3}=50$ Å to further
reduce out-of-plane interactions among periodic copies (see Fig. S2(a) in Ref.
[77] for an energy convergence test versus $a_{3}$). In order to treat the
strong on-site Coulomb interaction of Cr d orbitals, we used the DFT+U method
introduced by Dudarev et al. [78], with an effective value of $U$ equal to 3
eV ($J=0$) [79]. Structural optimizations were originally carried without
spin-orbit coupling (SOC); SOC was turned on for structures at local minima,
and to compute one-dimensional cuts of the energy landscape.

Having a well-converged electronic density is crucial to determine the
intrinsic electric dipole moment $\mathbf{P}$, which was calculated using the
Berry phase theory of polarization [80]. $\mathbf{P}$ remained unchanged on a
denser, $71\times 71\times 1$ $k-$point mesh. In addition, the convergence of
$\mathbf{P}$ against the magnitude of $\mathbf{a}_{3}$ was guaranteed as well
(see Fig. S2(b) for a test of $P_{z}$ versus $a_{3}$ [77]).

Phonon dispersion calculations were performed with the PHONOPY tool [81] on a
$3\times 3\times 1$ supercell. Atomic displacements were set at 0.005 $\AA{}$.
In those calculations, the $k-$point grid was set to $5\times 5\times 1$.
Cutoff energy and energy convergence criteria remained at 500 eV and $10^{-4}$
eV/$\AA{}$, respectively. The force constant (Hessian) matrix at the
$\Gamma-$point was obtained as well; eigenvalues and eigenvectors of this
Hessian were utilized to calculate the magnetoelectric coupling tensor
$\alpha$ using the method developed by Íñiguez [82].

##  III Results and discussion

###  III.1 AB and AB’ CrI3 bilayers

The CrI3 monolayers [83, 62] depicted on Fig. 1(a) are one of many available
two-dimensional ferromagnets and antiferromagnets [84, 85, 86, 87, 88, 89].
Crucial for an eventual creation of an in-plane electric dipole moment are the
voids within star-of-David patterns on this structure. No such voids are
present on hexagonal boron nitride nor transition metal dichalcogenide
monolayers, for which no in-plane intrinsic electric dipoles occur on their
(non-moiré, rotated by $60^{\circ}$) ferroelectric bilayer configurations [90,
91, 92, 92, 93, 94, 95, 96]. This monolayer belongs to the crystalline space
group $P\overline{3}1m$ (space group 162) [97]. (A listing of all its
potential magnetic groups can be found in Table SII [77].)

Figure 1: (a) CrI3 monolayer. (b) CrI3 bilayer on its AB FM ($uuuu$)
configuration prior to the relative rotation: FM coupling prevents the removal
of a center of inversion and the development of a net intrinsic electric
dipole; see Table 1. (c) CrI3 bilayer with AB’-stacking. The crystalline point
of inversion–not considering the magnetic configuration–is shown.

When analyzing CrI3 bilayers, previous studies explore two stable structures
(labeled AB and AB’) that are formed by relative translations between layers
[62, 63, 64, 70, 97, 98]: Starting with an AA structure in which the two
monolayers sit atop one another, translating the top layer by
$2\mathbf{a}_{1}/3+\mathbf{a}_{2}/3$ yields the AB structure displayed in Fig.
1(b). Not considering magnetism, it belongs to space group $P\overline{3}$
(space group 147). (Possible magnetic groups are listed in Table SII [77].) If
one translates the lower monolayer originally set on the AB conformation by
$\mathbf{a}_{1}/3$, the resultant structure is AB’ [Fig. 1(c)], which belongs
to the space group $C2/m$ (12) [97] without the inclusion of magnetism (see
Table SII for potential magnetic groups [77]).

The AB bilayer is observed in low-temperature samples ($T<210$ K) [99] and has
the magnetic spins of Cr atoms ferromagnetically (FM) aligned in an Ising-like
(out-of-plane) configuration. To set ideas, each monolayer has two Cr atoms on
its unit cell (navy blue dots on Fig. 1), and each Cr atom contributes a spin
magnetic dipole moment $\mu$ of about 3 Bohr magnetons ($\mu_{B}=9.274\times
10^{-24}$ J/T) along the $z-$axis. Without loss of generality, we set each of
the four Cr atoms on the unit cell (u.c.) to
$\boldsymbol{\mu}=(0,0,+3\mu_{B}$), and label such FM setting $uuuu$ (for up-
up-up-up), where the first two labels refer to the Cr atoms in the lower
monolayer and the remaining two labels indicate the direction of spin on the
two atoms at the upper monolayer. A displacement by either $\mathbf{a}_{1}/3$
or $\mathbf{a}_{2}/3$ takes the structure shown into the widely studied AB’
configuration [100, 70, 97, 98, 79, 67, 65, 83], with a preferred AFM
configuration which is explicitly labeled $uudd$ and in which Cr atoms at each
monolayer have out-of-plane magnetic dipoles pointing toward the opposite
monolayer ($d$ stands for spin “down” in this context). Without considering
their magnetic moment, antiferromagnetically coupled (AB’) CrI3 bilayers stem
from a bulk monoclinic phase and belong to space group $C2/m$ (space group
12). All possible magnetic groups with an AB’ configuration are listed in
Table SII [77].

The lack of a center of inversion on AFM AB’ bilayers was experimentally
verified through second harmonic generation and it is driven by their magnetic
configuration [67, 65, 83]. In other words, its non-zero magnetoelectric
coupling [67] originates from the removal of a center of inversion by an
antiparallel spin. Theory indicates a FM coupling among CrI3 monolayers in the
AB (ground state) configuration [70, 97, 98, 79], which is inconsistent with
experiment, where an AFM coupling is observed. That theoretical result is
independent of the exchange-correlation functional, the van der Waals
correction scheme, and the Hubbard (U) term on the Cr atoms with and without
spin-orbit coupling (SOC) [79]. The ground state structure at bilayer
thickness is still a matter of debate–it has been suggested that the AFM AB’
structure is created at room temperature and “trapped” by the extra layers
used to create experimental devices [70, 97, 98]. Nevertheless, switching
among AFM and FM magnetic configurations can be experimentally performed at
will [62, 63, 64].

###  III.2 Rotated CrI3 bilayers ($s_{1}$, $s_{2,1}$, and $s_{2,2}$)

While the electric control of magnetism in CrI3 bilayers reported in the
recent past [62, 63, 64] relies on an external electric field, one precept of
this work is to engineer an intrinsic electric dipole by a relative rotation
[61] to enhance the multiferroic properties of those bilayers; such rotations
are a new handle for magnetoelectric multiferroic behavior that is not
available on non-layered multiferroics [6, 20, 9, 10, 21, 7, 22, 11, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 2, 38, 39, 40, 41].

Following a concept introduced by Wu and collaborators [101, 102] (and
experimentally verified for hexagonal boron nitride and transition metal
dichalcogenide bilayers [90, 91, 92, 92, 93, 94, 95]), we rotated the upper
monolayer on the AB structure by $60^{\circ}$. This removes the bilayer’s
crystalline center of inversion, and induces an intrinsic electric dipole. The
rotated bilayer structure–which we labeled $s_{1}$ (for structure 1)–is
displayed on Fig. 1(c) and it belongs to the crystalline space group $Cm$
(space group 8) (see Table SII for its possible magnetic groups). The group
$Cm$ does not forbid the creation of an intrinsic in-plane electric dipole
that is parallel to the mirror plane depicted on Fig. 1(c) [61].

Figure 2: (a) CrI3 bilayer at a local minimum after a $60^{\circ}$ rotation
with respect to the axis shown. This bilayer structure is dubbed $s_{1}$. The
removal of a crystalline inversion center is made evident in the lowermost
side views. $P_{z}$ can be reverted by sliding the lower monolayer by (c)
about $2\mathbf{a}_{1}/3$ on a bilayer structure labelled $s_{2,1}$, or (d) by
about $2\mathbf{a}_{2}/3$, on a structure labeled $s_{2,2}$. The different
orientations of the mirror plane for structures $s_{2,1}$ and $s_{2,2}$ are
highlighted.

An additional displacement of the lower monolayer in bilayer $s_{1}$ [Fig.
2(a)] by about $2\mathbf{a}_{1}/3$ [as explicitly displayed on Fig. 2(b)] or
by about $2\mathbf{a}_{2}/3$ [Fig. 2(c)] still results on a structure
belonging to crystalline space group $Cm$ (8). We label $s_{2,1}$ the
structure seen on Fig. 2(b)–in which the lower monolayer was originally
displaced by $2\mathbf{a}_{1}/3$–to tell it apart from the bilayer $s_{2,2}$
shown in Fig. 2(c), in which the lower monolayer was initially displaced by
about $2\mathbf{a}_{2}/3$ instead.

A comparison among Figs. 2(a) and 2(b) shows that the mirror plane of bilayer
$s_{2,1}$ remains parallel to that seen on bilayer $s_{1}$. Thus, besides the
natural creation of an out-of-plane intrinsic electric polarization by the
$60^{\circ}$ relative rotation, this mirror plane allows for a non-zero in-
plane intrinsic electric dipole parallel to $\mathbf{a}_{1}$ on bilayer
$s_{2,1}$ as well. Nevertheless, the mirror plane turns out to be parallel to
the $x-$axis (this is, parallel to $\mathbf{a}_{1}+\mathbf{a}_{2}$) for the
$s_{2,2}$ bilayer [Fig. 2(c)], allowing for a possible in-plane intrinsic
electric dipole along the $x-$axis on the $s_{2,2}$ bilayer. This way, the
rotated CrI3 bilayers may develop an in-plane component of the intrinsic
electric dipole that can be controlled by sliding. The creation of an in-plane
intrinsic electric dipole through a relative $60^{\circ}$ rotation among
monolayers lends additional handles for the engineering of magnetoelectric
couplings by geometrical/mechanical means.

Fig. 3 shows the total energy and the global minima for the AB and (rotated)
$s_{1}$ FM bilayers, and for three additional AFM configurations ($uudd$,
$udud$, and $uddu$) for both AB and $s_{1}$ structures. The energetics of the
AB bilayer are consistent with previous results [97, 70]. The numerical
results reported on Fig. 3 were obtained without SOC. (Electronic band
structures for AB and $s_{1}$ bilayers show similar band gaps; those can be
found in Fig. S3 [77].) Note that–though the $s_{1}$ structure sits higher in
energy than the AB bilayer on Fig. 3(c)–“falling back” into the AB structure
requires a macroscopic rotation of the entire upper layer in unison, an event
that is statistically rare. In other words, structure $s_{1}$ sits at a stable
local energy minimum.

Figure 3: Side views illustrating four possible (Ising-like) magnetic
configurations for (a) the AB CrI3 bilayer and for (b) the rotated bilayer on
the $s_{1}$ configuration. The FM configuration is labeled $uuuu$, and three
AFM configurations ($uudd$, $udud$, and $uddu$) were considered as well. (c)
Total structural energy around global (AB) and local ($s_{1}$) minima for CrI3
bilayers as a function of the lattice parameter $a_{0}$ and magnetic
configuration, without SOC. The $s_{1}$ local minima on the FM configuration
sits 139 K/u.c. above the global minima, while the AFM $s_{1}$ structure with
a $uudd$ magnetic configuration sits 202 K/u.c. above that global minima (see
Table SIII for numerical details [77]). Figure 4: (a) Energy landscape as a
function of sliding among CrI3 monolayers for FM ($uuuu$) and AFM ($uudd$)
magnetic configurations (without SOC), and corresponding optimal vertical
separation among monolayers $D$ (see Fig. 1). Structure $s_{1}$ sits at the
point from which lattice vectors $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ join,
and it is periodic upon lattice translations. Structure $s_{2,1}$ ($s_{2,2}$)
sits $0.6389\mathbf{a}_{1}$ ($0.6806\mathbf{a}_{2}$) away from $s_{1}$
locations. (b) One-dimensional cuts of the energy landscape along (i)
$\mathbf{a}_{1}$ and (ii) $\mathbf{a}_{2}$ with and without SOC turned on. (c)
Structural stability of $s_{1}$ bilayer, and vibrational spectra at saddle
point $p$ displaying a few imaginary modes. (d) Induced intrinsic three-
dimensional electric dipole along the straight path joining structures $s_{1}$
and $s_{2,1}$. The inset shows the change in magnetization at a nearly
degenerate crossover among AFM and FM configurations at point $p$.

Fig. 4(a) displays the energy landscape and the optimal relative height among
monolayers upon relative sliding, on the bilayer configuration in which the
upper monolayer was rotated by $60^{\circ}$. This was calculated for FM
($uuuu$) and AFM ($uudd$) magnetic configurations without SOC. The point at
which lattice vectors $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ join in Fig. 4(a)
provides the energy and separation among monolayers for structure $s_{1}$ on
Figs. 1(c) and 3(c).

Working with a primitive u.c., the energy and structural sampling was carried
out by displacing the upper layer by
$n\mathbf{a}_{1}/12+m\mathbf{a}_{2}/12$–$n$ and $m$ ran from 0 to 11–with
respect to structure $s_{1}$. So that bilayers do not fall back into local
minima while computing the landscape, the atomic in-plane coordinates were
kept fixed, and only the out-of-plane components were allowed to relax. The
landscape and height profiles depicted upon multiple u.c.s on Fig. 4(a) rely
on periodicity, and the energies were expressed with respect to that obtained
for the AB FM configuration (see Fig. 3 and Table SIII [77]).

A more detailed calculation was performed along one-dimensional cuts along
$\mathbf{a}_{1}$ or $\mathbf{a}_{2}$ (Fig. 4(b)). In those calculations, all
atoms were allowed to move around the local minima, and SOC was turned on.
Unlike the energy landscape of rotated hBN or transition metal dichalcogenide
(TMDC) bilayers, which have identical maximum energies as the layers move
either along the $\mathbf{a}_{1}$ or $\mathbf{a}_{2}$ axes (see, e.g., Refs.
[96] and [103]), the rotated CrI3 bilayer displays barriers of dissimilar
height, with the largest height assigned to a displacement along
$\mathbf{a}_{1}$. This anisotropy was not discussed in Ref. [61]. The
identical barriers observed on rotated hBN and TMDC bilayers [96, 103]
guarantees that the local minima be exactly at a distance $2\mathbf{a}_{1}/3$
or $2\mathbf{a}_{2}/3$. To the contrary, and as seen on Fig. 4(b), the
dissimilar heights of the barriers makes the $s_{2,1}$ structure to be located
at 0.6389$\mathbf{a}_{1}$ instead, while the $s_{2,2}$ structure sits at
0.6806$\mathbf{a}_{2}$, and it justifies our choice of words (“about 2/3”)
made earlier on.

Calculations in which SOC is turned on require copious additional
computational resources, and those were only pursued for structures along the
one-dimensional cuts of the landscape in Fig. 4(b) for that reason.
Surprisingly, SOC did not break the energy degeneracy of the minima $s_{1}$
and $s_{2,1}$ upon inversion of the magnetic moment significantly, as the
total energy difference upon a magnetic swap from FM configurations $uuuu$ to
$dddd$ was only 0.15 K/u.c. Calculations of magnetic anisotropy barriers can
be found on Fig. S4 [77].

The structural stability of bilayer $s_{1}$ alluded to earlier is established
with the vibrational spectra shown in Fig. 4(c), which lacks imaginary modes.
In addition, and in preparation for the calculation of the magnetoelectric
coupling tensor, we also considered the bilayer structure situated at the
(unstable) saddle point $p$ in the landscape (see imaginary modes), for which
both FM and AFM magnetic configurations turn degenerate. Point $p$ is shown on
Fig. 4(b). Fig. 4(d) shows $\mathbf{P}$ as the structure is turned from
$s_{2,1}$ to $s_{1}$ along $\mathbf{a}_{1}$. For comparison, the reported
magnitudes of the out-of-plane intrinsic electric dipole are 1.0 pC/m for a
MoS2 bilayer [101], 0.5 pC/m for a WSe2 bilayer [95], and 0.4 pC/m for a WTe2
bilayer [102, 104].

###  III.3 Three-dimensional intrinsic electric polarization of rotated CrI3
bilayers

While $P_{z}$ swaps sign as the bilayer evolves from structure $s_{1}$ into
$s_{2,1}$ or $s_{2,2}$ as expected [101, 102], there is an in-plane intrinsic
electric dipole moment in Fig. 4(d). In the remainder of this work we will (i)
justify this observation, (ii) argue that it can be switched, and (iii)
calculate the magnetoelectric coupling of the rotated CrI3 bilayer.

We start by reporting the magnitude of the intrinsic electric dipoles for the
AB and AB’ bilayers in Table 1. For context, electric dipole moments as large
as 5 pC/m have been recently reported experimentally [105].

Table 1: Intrinsic electric dipole moments per u.c. area for CrI3 bilayers at AB and AB’ configurations. FM configurations have all spins pointing up ($uuuu$) or down ($dddd$), while the studed AFM configurations were $uudd$ and $dduu$. Digits within parenthesis indicate the precision of calculations. Configuration |  $\mathbf{P}$ (pC/m)
---|---
AB FM ($uuuu$) | $(-0.0(1),0.0(0),-0.0(1))$
AB FM ($dddd$) | $(-0.0(1),0.0(0),-0.0(1))$
AB’1 FM ($uuuu$) | $(-0.0(1),0.0(0),-0.0(1))$
AB’1 FM ($dddd$) | $(-0.0(1),0.0(0),-0.0(1))$
AB’2 FM ($uuuu$) | $(-0.0(1),0.0(0),-0.0(1))$
AB’2 FM ($dddd$) | $(-0.0(1),0.0(0),-0.0(1))$
Configuration |  $\mathbf{P}$ (pC/m)
AB AFM ($uudd$) | $(-0.0(1),0.0(0),-0.0(3))$
AB AFM ($dduu$) | $(-0.0(1),0.0(0),+0.0(2))$
AB’1 AFM ($uudd$) | $(-0.0(1),0.0(0),-0.0(4))$
AB’1 AFM ($dduu$) | $(-0.0(1),0.0(0),+0.0(2))$
AB’2 AFM ($uudd$) | $(-0.0(1),0.0(0),-0.0(4))$
AB’2 AFM ($dduu$) | $(-0.0(1),0.0(1),+0.0(2))$

Table 1 conveys the following information: (i) FM-coupled AB and AB’ bilayers
display negligible intrinsic electric dipole moments. Since there is a center
of inversion in those structures, the reported $\mathbf{P}$ on for the FM-
coupled structures conveys the numerical precision achieved here, which is of
the order of $\pm 0.01$ pC/m.

On the other hand, the lack of a center of inversion facilitated by AFM
magnetic configurations gives rise to an expected net $P_{z}$ [67, 65, 83] in
Table 1. In AFM structures, $P_{z}$ does not change sign in going from
structure AB onto any of the two energy-degenerate structures AB’, which are
arrived at by a displacement of the upper monolayer by $\mathbf{a}_{1}/3$ (for
a structure labelled AB’1) or $\mathbf{a}_{2}/3$ (AB’2). The two structures
AB’1 and AB’2 were optimized independently, hence providing a cross-check of
numerical results. In short, FM AB and AB’ configurations have a center of
inversion and do not develop an intrinsic electric dipole; AFM AB and AB’
configurations, on the other hand, break inversion symmetry due to their spin
configuration and give rise to a non-ionic out-of-plane net intrinsic electric
dipole moment ($P_{z}$) that swaps sign upon spin reversal and regardless of
the (AB or AB’) atomistic configuration. We are unaware of any previous
computational study discussing the magnitude of the intrinsic electric dipole
on experimentally available CrI3 bilayer configurations as the one provided in
Table 1.

The effectiveness of the relative $60^{\circ}$ rotation procedure in enhancing
$P_{z}$ is demonstrated in Table 2, where a ten-fold increase with respect to
the experimentally available AB’ AFM bilayer is registered. Rounding up the
reported values to 0.3 pm/C (a value close to the 0.2 pm/C reported in Ref.
[106]), and dividing by $\simeq$12 Å–which would correspond to the thickness
of a u.c. in a bulk 3D stack–one gets 0.025 $\mu$C/cm2, which sits halfway
among the values reported in Refs. [61] and [107] and may be mostly due to the
different exchange-correlation functionals employed in our works. Furthermore,
switching of $P_{z}$ is achieved by sliding from structure $s_{1}$ into either
$s_{2,1}$ or $s_{2,2}$ [101, 102]. A more remarkable result is the creation of
an even larger in-plane intrinsic electric dipole like the one present on
group-IV monochalcogenide monolayers [47, 108, 49].

Table 2: Intrinsic electric dipole moments per u.c. area for CrI3 bilayers at the $s_{1}$, $s_{2,1}$, and $s_{2,2}$ configurations. Digits within parenthesis indicate the precision of calculations. The removal of a center of inversion upon rotation makes $P_{z}$ to increase an order of magnitude with respect to its value on AFM AB and AB’ configurations (Table 1). In addition, $P_{z}$ does not change sign upon reversal of antiferromagnetic polarization ($uudd$ versus $dduu$), but only upon relative displacements among monolayers. A larger in-plane intrinsic electric dipole swapping direction at $s_{2,1}$ and $s_{2,2}$ is now created. Configuration |  $\mathbf{P}$ (pC/m)
---|---
$s_{1}$ FM ($uuuu$) | $(+0.9(3),-1.6(0),+0.2(5))$
$s_{1}$ FM ($dddd$) | $(+0.9(3),-1.6(0),+0.3(0))$
$s_{2,1}$ FM ($uuuu$) | $(+0.9(4),-1.6(0),-0.2(5))$
$s_{2,1}$ FM ($dddd$) | $(+0.9(4),-1.6(1),-0.3(1))$
$s_{2,2}$ FM ($uuuu$) | $(-1.8(4),+0.0(1),-0.2(5))$
$s_{2,2}$ FM ($dddd$) | $(-1.8(7),+0.0(1),-0.3(0))$
Configuration |  $\mathbf{P}$ (pC/m)
$s_{1}$ AFM ($uudd$) | $(+0.8(0),-1.3(7),+0.2(5))$
$s_{1}$ AFM ($dduu$) | $(+0.8(0),-1.3(7),+0.2(5))$
$s_{2,1}$ AFM ($uudd$) | $(+0.8(0),-1.3(7),-0.2(5))$
$s_{2,1}$ AFM ($dduu$) | $(+0.8(0),-1.3(7),-0.2(5))$
$s_{2,2}$ AFM ($uudd$) | $(-1.6(1),+0.0(1),-0.2(5))$
$s_{2,2}$ AFM ($dduu$) | $(-1.6(1),+0.0(1),-0.2(5))$

We explain the existence of an in-plane intrinsic electric dipole (which
provides new vistas for magnetoelectric couplings on this novel multiferroic)
by recourse to symmetry: $s_{1}$, $s_{2,1}$, and $s_{2,2}$ CrI3 bilayers all
belong to crystalline space group $Cm$. Nevertheless, while the $s_{1}$ and
$s_{2,1}$ bilayers have a mirror plane parallel to
$\mathbf{a}_{1}=a_{0}\left(1/2,-\sqrt{3}/2,0\right)$ (Figures 1(c) and 1(d)),
bilayer $s_{2,2}$ has a mirror plane parallel to
$\mathbf{a}_{1}+\mathbf{a}_{2}$ ($x-$axis on Fig. 2). Any in-plane dipole for
structures $s_{1}$ through $s_{2,1}$, if it is to exist, must be parallel to
$\mathbf{a}_{1}$ (and this is the case for $P_{x}$ and $P_{y}$ on Fig. 4(d)).
In other words, its magnitude must satisfy the ratio $-\sqrt{3}\simeq-1.73$.
This is precisely documented in Table 3. As for structure $s_{2,2}$, the in-
plane electric dipole can only be parallel to the $x-$axis, which is what one
sees in Table 2. Furthermore, the magnitude of the in-plane dipole is nearly
identical at $s_{1}$, $s_{2,1}$, and $s_{2,2}$, and it must switch in going
from structure $s_{2,1}$ to structure $s_{2,2}$ through $s_{1}$.

The second check to justify the creation of in-plane components of intrinsic
electric polarization is numerical: adding up all vectors that can be built
from all four Cr atoms to the nearest I atoms on the opposite monolayer within
the u.c. ($\Delta_{x},\Delta_{y},\Delta_{z}$), one gets the results listed in
Table 3. A top view of the atoms involved is displayed in Fig. 5 for bilayers
AB, AB’1, $s_{1}$, and $s_{2,1}$. The observed relative displacements indeed
respect the crystalline symmetry of group $Cm$: they are parallel to the
mirror planes.

Table 3: Sum of vectors pointing from Cr atoms to the first and second nearest I neighbors located in the opposite monolayer. This sum is carried out over all Cr atoms in the u.c. $P_{y}/P_{x}$ as reported in Table 2 correlates with $\Delta_{y}/\Delta_{x}$, and hence justifies the existence of a non-zero in-plane intrinsic dipole moments on rotated CrI3 bilayers. (${a_{1,y}}/{a_{1,x}}=-\sqrt{3}\simeq-1.73$.) Configuration | ($\Delta_{x},\Delta_{y},\Delta_{z}$) ($\AA{}$/u.c.) | $P_{y}/P_{x}$ | $\Delta_{y}/\Delta_{x}$
---|---|---|---
AB FM ($uuuu$) | $(0.000,0.000,0.000)$ | - | -
AB AFM ($uudd$) | $(0.000,0.000,0.000)$ | - | -
$s_{1}$ FM ($uuuu$) | $(+0.216,-0.388,+10.372)$ | $-1.79$ | $-1.72$
$s_{1}$ AFM ($uudd$) | $(+0.216,-0.382,+10.374)$ | $-1.77$ | $-1.72$
$s_{21}$ FM ($uuuu$) | $(+0.221,-0.410,-10.369)$ | $-1.70$ | $-1.85$
$s_{21}$ AFM ($uudd$) | $(+0.247,-0.423,-10.306)$ | $-1.71$ | $-1.71$
$s_{22}$ FM ($uuuu$) | $(-0.453,-0.021,-10.373)$ | $0.00$ | $+0.05$
$s_{22}$ AFM ($uudd$) | $(-0.453,-0.021,-10.373)$ | $0.00$ | $+0.05$
Figure 5: Schematic top views of Cr atoms and nearest I atoms on the opposite
monolayer; vectors joining those atoms are utilized to determine
$(\Delta_{x},\Delta_{y},\Delta_{z})$ in Table 3.

###  III.4 Calculation of the magnetoelectric coupling at point $p$ within the
energy landscape

We end this work with a calculation of the magnetoelectric tensor for the
rotated CrI3 bilayer. Following Íñiguez [82], we seek a point in the landscape
for which both intrinsic electric and magnetic dipoles are zero. As seen on
Fig. 4, the in-plane components of the intrinsic polarization remain constant,
and so we used the change of $\mathbf{P}$ ($\Delta\mathbf{P}$) away from the
polarization at point $p$ here. The eigenvalues $C_{n}$ of the bilayer’s
Hessian at $p$ listed in Table 4 range in between $-$0.011 and 13.556 eV/Å,
which compare well with Ref. [82]. We calculated the intrinsic electric dipole
for eigenvalues greater than zero (IR modes get revealed by giving non-zero
net electric polarizations), and set all associated eigenvectors to be $u=0.1$
Å in magnitude.

Table 4:  Hessian eigenvalues for rotated CrI3 bilayer at point $p$ (eV/Å). $-$0.011 | 0.000 | 0.000 | 0.000 | 0.007 | 0.263
---|---|---|---|---|---
0.917 | 1.022 | 1.088 | 1.110 | 1.456 | 1.482
2.413 | 2.457 | 2.457 | 2.483 | 2.497 | 2.524
2.714 | 2.731 | 3.552 | 3.614 | 3.629 | 3.672
4.002 | 4.089 | 4.138 | 4.183 | 4.831 | 4.900
4.935 | 4.963 | 6.137 | 6.177 | 6.203 | 6.267
9.319 | 9.356 | 9.830 | 9.917 | 9.935 | 10.047
10.836 | 10.944 | 10.971 | 11.098 | 13.556 | 13.556

Since $\boldsymbol{\mu}$ changes stepwise from (0,0,0) to (0,0,12$\mu_{B}$)
away from $p$, we used a value for the magnetic polarity equal to
$\mathbf{p}^{m}=(0,0,\frac{12\mu_{B}\mu_{0}}{u})$, where $\mu_{0}=1.2566\times
10^{-6}$ N/A2 is the permeability of vacuum. This way, our expression for the
magnetoelectric tensor takes the form (SI units):

$\displaystyle\alpha_{3,j}=\alpha_{j,3}=\frac{12\mu_{B}\mu_{0}}{0.01\text{\AA}^{2}}\sum_{n=5}^{48}\frac{\Delta
P_{n,j}}{C_{n}}$ (1) $\displaystyle=(-31.7,6.7,-4.3)\times 10^{-20}\text{ s},$

where the lack of 1/m reflects the fact that $\mathbf{P}$ is reported in units
of charge per length here. Dividing by the bilayer thickness ($\simeq 10$ Å)
and multiplying by the speed of light [109], one gets

$\alpha=\begin{pmatrix}0.00&0.00&-9.30\\\ 0.00&0.00&2.01\\\
-9.30&2.01&-1.29\end{pmatrix}\times 10^{-2}$ (2)

in gaussian units. This way, $\alpha$ here is one order of magnitude smaller
than the values reported for Cr2O3 by Íñiguez.

CrI3 bilayers are experimentally known to become paramagnetic at a few Kelvin.
$AB^{\prime}$ bilayers will lose their ferroelectric properties at that point,
while rotated $s_{1}$, $s_{2,1}$ or $s_{2,2}$ bilayers may still maintain
their ferroelectric properties at temperatures well above those for which
their intrinsic magnetism is lost.

The magnetoelectric tensor thus found implies that the out-of-plane magnetic
moment can be most easily switched by applying an electric field along the
$x-$direction, or that the in-plane $x-$component of the intrinsic electric
dipole will switch direction most easily as the Ising (out-of-plane) spin
swaps sign. The explicit calculation of $\alpha$ brings the engineering of
two-dimensional magnetoelectric multiferroics on rotated magnetic bilayers
onto a more quantitative footing.

##  IV Conclusion

Seeking to enhance the out-of-plane intrinsic electric dipole of CrI3 bilayers
by a relative rotation, we discovered that those bilayers develop an
unexpected and sizeable in-plane intrinsic electric dipole moment as well.
This is due to the unique atomistic structure of those magnetic materials that
does not forbid the creation of in-plane intrinsic electric dipoles. The
magnetoelectric coupling tensor $\alpha$ was estimated explicitly as well.
Those results enhance our understanding, and provide new vistas into
magnetoelectric multiferroics on two-dimensional material platforms.

###### Acknowledgements.

Calculations from the Arkansas team were performed at Cori at NERSC, a DOE
facility funded under Contract DE-AC02-05CH11231, and at the University of
Arkansas’ Pinnacle supercomputer, funded by the U.S. National Science
Foundation, the Arkansas Economic Development Commission, and the Office of
the Vice Provost for Research and Innovation. J.M.M.T. and M.A.M. thank
Montana State University, Bozeman, for startup support and computational
resources within the Tempest Research Cluster. All authors acknowledge
financial support from the MonArk NSF Quantum Foundry, supported by the
National Science Foundation Q-AMASE-i program under NSF Award DMR-1906383.

## References

  * Curie [1894] M. P. Curie, Sur la symétrie dans les phénomènes physiques, symétrie d’un champ électrique et d’un champ magnétique, J. Phys. Theor. Appl. 3, 393 (1894).
  * Smolenskii and Chupis [1982] G. A. Smolenskii and I. E. Chupis, Ferroelectromagnets, Sov. phys., Usp. 25, 475 (1982).
  * Fiebig  _et al._ [2013] M. Fiebig, V. Eremenko, and I. Chupis,  _Magnetoelectric Interaction Phenomena in Crystals_ , NATO Science Series II: Mathematics, Physics and Chemistry (Springer Netherlands, 2013).
  * Landau and Lifshitz [1984] L. D. Landau and E. M. Lifshitz,  _Electrodynamics of Continuous Media_ (Pergamon, New York, 1984).
  * Astrov [1960] D. N. Astrov, The magnetoelectric effect in antiferromagnets, J. Exptl. Theoret. Phys. (U.S.S.R.) 38, 984 (1960).
  * Kimura  _et al._ [2003] T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y. Tokura, Magnetic control of ferroelectric polarization, Nature 426, 55 (2003).
  * Kimura  _et al._ [2005] T. Kimura, G. Lawes, T. Goto, Y. Tokura, and A. P. Ramirez, Magnetoelectric phase diagrams of orthorhombic RMnO3 (R=Gd, Tb, and Dy), Phys. Rev. B 71, 224425 (2005).
  * Eerenstein  _et al._ [2006] W. Eerenstein, N. D. Mathur, and J. F. Scott, Multiferroic and magnetoelectric materials, Nature 442, 759 (2006).
  * Fiebig [2005] M. Fiebig, Revival of the magnetoelectric effect, J. Phys. D: Appl. Phys. 38, R123 (2005).
  * Hill [2000] N. A. Hill, Why are there so few magnetic ferroelectrics?, J. Phys. Chem. B 104, 6694 (2000).
  * Van Aken  _et al._ [2004] B. B. Van Aken, T. T. Palstra, A. Filippetti, and N. A. Spaldin, The origin of ferroelectricity in magnetoelectric YMnO3, Nat. Mater. 3, 164 (2004).
  * Heron  _et al._ [2011] J. T. Heron, M. Trassin, K. Ashraf, M. Gajek, Q. He, S. Y. Yang, D. E. Nikonov, Y.-H. Chu, S. Salahuddin, and R. Ramesh, Electric-field-induced magnetization reversal in a ferromagnet-multiferroic heterostructure, Phys. Rev. Lett. 107, 217202 (2011).
  * Novoselov  _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Two-dimensional gas of massless Dirac fermions in graphene, Nature 438, 197 (2005).
  * Zhang  _et al._ [2005] Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Experimental observation of the quantum Hall effect and Berry’s phase in graphene, Nature 438, 201 (2005).
  * Katsnelson [2012] M. I. Katsnelson,  _Graphene: Carbon in Two Dimensions_ (Cambridge University Press, 2012).
  * Röntgen [1888] W. C. Röntgen, Ueber die durch bewegung eines im homogenen electrischen felde befindlichen dielectricums hervorgerufene electrodynamische kraft, Ann. Phys. 271, 264 (1888).
  * Debye [1926] P. Debye, Bemerkung zu einigen neuen versuchen über einen magneto-elektrischen richteffekt, Z. Physik 36, 300 (1926).
  * Smolenskii and Agranovskaya [1960] G. A. Smolenskii and A. I. Agranovskaya, Dielectric polarization of a number of complex compounds, Sov. phys., Solid state 1, 1429 (1960).
  * Dzyaloshinskii [1959] I. E. Dzyaloshinskii, On the magneto-electrical effects in antiferromagnets, J. Exptl. Theoret. Phys. (U.S.S.R.) 10, 628 (1959).
  * Wang  _et al._ [2003] J. Wang, J. B. Neaton, H. Zheng, V. Nagarajan, S. B. Ogale, B. Liu, D. Viehland, V. Vaithyanathan, D. G. Schlom, U. V. Waghmare, N. A. Spaldin, K. M. Rabe, M. Wuttig, and R. Ramesh, Epitaxial BiFeO3 multiferroic thin film heterostructures, Science 299, 1719 (2003).
  * Neaton  _et al._ [2005] J. B. Neaton, C. Ederer, U. V. Waghmare, N. A. Spaldin, and K. M. Rabe, First-principles study of spontaneous polarization in multiferroic $\mathrm{Bi}\mathrm{Fe}{\mathrm{O}}_{3}$, Phys. Rev. B 71, 014113 (2005).
  * Spaldin and Fiebig [2005] N. A. Spaldin and M. Fiebig, The renaissance of magnetoelectric multiferroics, Science 309, 391 (2005).
  * Duan  _et al._ [2006] C.-G. Duan, S. S. Jaswal, and E. Y. Tsymbal, Predicted magnetoelectric effect in Fe/BaTiO3 multilayers: Ferroelectric control of magnetism, Phys. Rev. Lett. 97, 047201 (2006).
  * Fennie and Rabe [2006] C. J. Fennie and K. M. Rabe, Magnetic and electric phase control in epitaxial ${\mathrm{EuTiO}}_{3}$ from first principles, Phys. Rev. Lett. 97, 267602 (2006).
  * Sergienko and Dagotto [2006] I. A. Sergienko and E. Dagotto, Role of the Dzyaloshinskii-Moriya interaction in multiferroic perovskites, Phys. Rev. B 73, 094434 (2006).
  * Picozzi  _et al._ [2007] S. Picozzi, K. Yamauchi, B. Sanyal, I. A. Sergienko, and E. Dagotto, Dual nature of improper ferroelectricity in a magnetoelectric multiferroic, Phys. Rev. Lett. 99, 227201 (2007).
  * Tsymbal  _et al._ [2012a] E. Y. Tsymbal, E. R. A. Dagotto, C.-B. Eom, and R. Ramesh, eds.,  _Multifunctional Oxide Heterostructures_ (Oxford University Press, 2012).
  * Lee and Rabe [2011] J. H. Lee and K. M. Rabe, Coupled magnetic-ferroelectric metal-insulator transition in epitaxially strained ${\mathrm{srcoo}}_{3}$ from first principles, Phys. Rev. Lett. 107, 067601 (2011).
  * Valencia  _et al._ [2011] S. Valencia, A. Crassous, L. Bocher, V. Garcia, X. Moya, R. O. Cherifi, C. Deranlot, K. Bouzehouane, S. Fusil, A. Zobelli, A. Gloter, N. D. Mathur, A. Gaupp, R. Abrudan, F. Radu, A. Barthélémy, and M. Bibes, Interface-induced room-temperature multiferroicity in BaTiO3, Nat. Mater. 10, 753 (2011).
  * Tsymbal  _et al._ [2012b] E. Y. Tsymbal, A. Gruverman, V. Garcia, M. Bibes, , and A. Barthélémy, Ferroelectric and multiferroic tunnel junctions, MRS Bulletin 37, 138 (2012b).
  * Zhao  _et al._ [2014] H. J. Zhao, W. Ren, Y. Yang, J. Iñiguez, X. M. Chen, and L. Bellaiche, Near room-temperature multiferroic materials with tunable ferromagnetic and electrical properties, Nat. Commun. 5, 4021 (2014).
  * Spaldin and Ramesh [2019] N. Spaldin and R. Ramesh, Advances in magnetoelectric multiferroics, Nat. Mater. 18, 203 (2019).
  * Huang  _et al._ [2020] Y.-L. Huang, D. Nikonov, C. Addiego, R. V. Chopdekar, B. Prasad, L. Zhang, J. Chatterjee, H.-J. Liu, A. Farhan, Y.-H. Chu, M. Yang, M. Ramesh, Z. Q. Qiu, B. D. Huey, C.-C. Lin, T. Gosavi, J. Iñiguez, J. Bokor, X. Pan, I. Young, L. W. Martin, and R. Ramesh, Manipulating magnetoelectric energy landscape in multiferroics, Nat. Commun. 11, 2836 (2020).
  * Ponet  _et al._ [2022] L. Ponet, S. Artyukhin, T. Kain, J. Wettstein, A. Pimenov, A. Shuvaev, X. Wang, S.-W. Cheong, M. Mostovoy, and A. Pimenov, Topologically protected magnetoelectric switching in a multiferroic, Nature 607, 81 (2022).
  * Li  _et al._ [2022] P. Li, X. Zhou, and Z. Guo, Intriguing magnetoelectric effect in two-dimensional ferromagnetic/perovskite oxide ferroelectric heterostructure, npj Comput. Mater. 8, 20 (2022).
  * Shen  _et al._ [2023] X. Shen, F. Wang, X. Lu, and J. Zhang, Two-dimensional multiferroics with intrinsic magnetoelectric coupling in A-Site ordered perovskite monolayers, Nano Lett. 23, 735 (2023).
  * Schmid [1973] H. Schmid, On a magnetoelectric classification of materials, Int. J. Magnetism 4, 337 (1973).
  * Sugie  _et al._ [2002] H. Sugie, N. Iwata, and K. Kohn, Magnetic ordering of rare earth ions and magnetic-electric interaction of hexagonal RMnO3 (R=Ho, Er, Yb or Lu), J. Phys. Soc. Japan 71, 1558 (2002).
  * Schmid  _et al._ [1965] H. Schmid, H. Rieder, and E. Ascher, Magnetic susceptibilities of some 3d transition metal boracites, Solid State Commun. 3, 327 (1965).
  * Scott [1980] J. F. Scott, Spectroscopy of incommensurate ferroelectrics, Ferroelectrics 24, 127 (1980).
  * Aizu [1970] K. Aizu, Possible species of ferromagnetic, ferroelectric, and ferroelastic crystals, Phys. Rev. B 2, 754 (1970).
  * Seixas  _et al._ [2016] L. Seixas, A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Multiferroic two-dimensional materials, Phys. Rev. Lett. 116, 206803 (2016).
  * Bishop  _et al._ [2019] T. B. Bishop, E. E. Farmer, A. Sharmin, A. Pacheco-Sanjuan, P. Darancet, and S. Barraza-Lopez, Quantum paraelastic two-dimensional materials, Phys. Rev. Lett. 122, 015703 (2019).
  * Pacheco-Sanjuan  _et al._ [2019] A. Pacheco-Sanjuan, T. B. Bishop, E. E. Farmer, P. Kumar, and S. Barraza-Lopez, Evolution of elastic moduli through a two-dimensional structural transformation, Phys. Rev. B 99, 104108 (2019).
  * Mehboudi  _et al._ [2016] M. Mehboudi, A. M. Dorio, W. Zhu, A. van der Zande, H. O. H. Churchill, A. A. Pacheco-Sanjuan, E. O. Harriss, P. Kumar, and S. Barraza-Lopez, Two-dimensional disorder in black phosphorus and monochalcogenide monolayers, Nano Lett. 16, 1704 (2016).
  * Wang and Qian [2017] H. Wang and X. Qian, Two-dimensional multiferroics in monolayer group IV monochalcogenides, 2D Mater. 4, 015042 (2017).
  * Chang  _et al._ [2016] K. Chang, J. Liu, H. Lin, N. Wang, K. Zhao, A. Zhang, F. Jin, Y. Zhong, X. Hu, W. Duan, Q. Zhang, L. Fu, Q.-K. Xue, X. Chen, and S.-H. Ji, Discovery of robust in-plane ferroelectricity in atomic-thick SnTe, Science 353, 274 (2016).
  * Barraza-Lopez  _et al._ [2018] S. Barraza-Lopez, T. P. Kaloni, S. P. Poudel, and P. Kumar, Tuning the ferroelectric-to-paraelectric transition temperature and dipole orientation of group-IV monochalcogenide monolayers, Phys. Rev. B 97, 024110 (2018).
  * Barraza-Lopez  _et al._ [2021] S. Barraza-Lopez, B. M. Fregoso, J. W. Villanova, S. S. P. Parkin, and K. Chang, Colloquium: Physical properties of group-IV monochalcogenide monolayers, Rev. Mod. Phys. 93, 011001 (2021).
  * Cheng  _et al._ [2021] H.-X. Cheng, J. Zhou, C. Wang, W. Ji, and Y.-N. Zhang, Nonvolatile electric field control of magnetism in bilayer CrI3 on monolayer In2Se3, Phys. Rev. B 104, 064443 (2021).
  * Hu  _et al._ [2022] C. Hu, J. Chen, E. Du, W. Ju, Y. An, and S.-J. Gong, Ferroelectric control of band alignments and magnetic properties in the two-dimensional multiferroic VSe2/In2Se3, J. Phys.: Condens. Matter 34, 425801 (2022).
  * Zhu  _et al._ [2022] H. Zhu, Y. Gao, Y. Hou, Z. Gui, and L. Huang, Tunable magnetic anisotropy in two-dimensional CrX3/AlN (X=I, Br, Cl) heterostructures, Phys. Rev. B 106, 134412 (2022).
  * Gong  _et al._ [2019] C. Gong, E. M. Kim, Y. Wang, G. Lee, and X. Zhang, Multiferroicity in atomic van der Waals heterostructures, Nat. Commun. 10, 2657 (2019).
  * Lu  _et al._ [2022] Y. Lu, H. Wang, L. Wang, and L. Yang, Mechanism of carrier doping induced magnetic phase transitions in two-dimensional materials, Phys. Rev. B 106, 205403 (2022).
  * Verzhbitskiy  _et al._ [2020] I. A. Verzhbitskiy, H. Kurebayashi, H. Cheng, J. Zhou, S. Khan, Y. P. Feng, and G. Eda, Controlling the magnetic anisotropy in Cr2Ge2Te6 by electrostatic gating, Nat. Electron. 3, 460 (2020).
  * Gao  _et al._ [2021] Y. Gao, M. Gao, and Y. Lu, Two-dimensional multiferroics, Nanoscale 13, 19324 (2021).
  * Yīng and Zülicke [2022] Y. Yīng and U. Zülicke, Magnetoelectricity in two-dimensional materials, Adv. Phys.: X 7, 2032343 (2022).
  * Fumega and Lado [2022] A. O. Fumega and J. L. Lado, Moiré-driven multiferroic order in twisted CrCl3, CrBr3 and CrI3 bilayers (2022), arXiv:2207.01416.
  * Ju  _et al._ [2021] H. Ju, Y. Lee, K.-T. Kim, I. H. Choi, C. J. Roh, S. Son, P. Park, J. H. Kim, T. S. Jung, J. H. Kim, K. H. Kim, J.-G. Park, and J. S. Lee, Possible persistence of multiferroic order down to bilayer limit of van der waals material NiI2, Nano Lett. 21, 5126 (2021).
  * Song  _et al._ [2022] Q. Song, C. A. Occhialini, E. Ergeçen, B. Ilyas, D. Amoroso, P. Barone, J. Kapeghian, K. Watanabe, T. Taniguchi, A. S. Botana, S. Picozzi, N. Gedik, and R. Comin, Evidence for a single-layer van der Waals multiferroic, Nature 602, 601 (2022).
  * Ji  _et al._ [2023] J. Ji, G. Yu, C. Xu, and H. J. Xiang, General theory for bilayer stacking ferroelectricity, Phys. Rev. Lett. 130, 146801 (2023).
  * Huang  _et al._ [2018] B. Huang, G. Clark, D. R. Klein, D. MacNeill, E. Navarro-Moratalla, K. L. Seyler, N. Wilson, M. A. McGuire, D. H. Cobden, D. Xiao, W. Yao, P. Jarillo-Herrero, and X. Xu, Electrical control of 2D magnetism in bilayer CrI3, Nat. Nanotech. 13, 544 (2018).
  * Jiang  _et al._ [2018a] S. Jiang, L. Li, Z. Wang, K. F. Mak, and J. Shan, Controlling magnetism in 2D CrI3 by electrostatic doping, Nat. Nanotech. 13, 549 (2018a).
  * Jiang  _et al._ [2018b] S. Jiang, J. Shan, and K. F. Mak, Electric-field switching of two-dimensional van der Waals magnets, Nat. Mater. 17, 406 (2018b).
  * Mak  _et al._ [2019] K. F. Mak, J. Shan, and D. Ralph, Probing and controlling magnetic states in 2D layered magnetic materials, Nat. Rev. Phys. 1, 646 (2019).
  * Suárez Morell  _et al._ [2019] E. Suárez Morell, A. León, R. Hiroki Miwa, and P. Vargas, Control of magnetism in bilayer CrI3 by an external electric field, 2D Mater. 6, 025020 (2019).
  * Sun  _et al._ [2019] Z. Sun, Y. Yi, T. Song, G. Clark, B. Huang, Y. Shan, S. Wu, D. Huang, C. Gao, Z. Chen, M. McGuire, T. Cao, D. Xiao, W.-T. Liu, W. Yao, X. Xu, and S. Wu, Giant nonreciprocal second-harmonic generation from antiferromagnetic bilayer CrI3, Nature 572, 497 (2019).
  * Sivadas  _et al._ [2016] N. Sivadas, S. Okamoto, and D. Xiao, Gate-controllable magneto-optic kerr effect in layered collinear antiferromagnets, Phys. Rev. Lett. 117, 267203 (2016).
  * Webster and Yan [2018] L. Webster and J.-A. Yan, Strain-tunable magnetic anisotropy in monolayer CrCl3, CrBr3, and CrI3, Phys. Rev. B 98, 144411 (2018).
  * León  _et al._ [2020] A. M. León, J. W. González, J. Mejía-López, F. Crasto de Lima, and E. Suárez Morell, Strain-induced phase transition in CrI3 bilayers, 2D Mater. 7, 035008 (2020).
  * Liu  _et al._ [2023] K. Liu, X. Ma, S. Xu, Y. Li, and M. Zhao, Tunable sliding ferroelectricity and magnetoelectric coupling in two-dimensional multiferroic MnSe materials, npj Comput. Mater. 9, 16 (2023).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132, 154104 (2010).
  * [77] Supplementary materials contain: (i) an analysis of structural parameters against exchange correlation functionals and SOC, (ii) a table of possible magnetic groups, (iii) tests of energy and intrinsic electric dipole convergence against the vacuum in the unit cell and $k-$point meshes, (iv) values of lattice parameters and energies at global and local minima, (v) electronic band structures for AB and $s_{1}$ bilayers, and (vi) a calculation of magnetic anisotropy barriers.
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B 57, 1505 (1998).
  * Jiang  _et al._ [2019] P. Jiang, C. Wang, D. Chen, Z. Zhong, Z. Yuan, Z.-Y. Lu, and W. Ji, Stacking tunable interlayer magnetism in bilayer CrI3, Phys. Rev. B 99, 144401 (2019).
  * Vanderbilt [2018] D. Vanderbilt,  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, 2018).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Íñiguez [2008] J. Íñiguez, First-principles approach to lattice-mediated magnetoelectric effects, Phys. Rev. Lett. 101, 117201 (2008).
  * Huang  _et al._ [2017] B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit, Nature 546, 270 (2017).
  * Gong  _et al._ [2017] C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia, and X. Zhang, Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals, Nature 546, 265 (2017).
  * Kim  _et al._ [2019] M. Kim, P. Kumaravadivel, J. Birkbeck, W. Kuang, S. G. Xu, D. G. Hopkinson, J. Knolle, P. A. McClarty, A. I. Berdyugin, M. Ben Shalom, R. V. Gorbachev, S. J. Haigh, S. Liu, J. H. Edgar, K. S. Novoselov, I. V. Grigorieva, and A. K. Geim, Micromagnetometry of two-dimensional ferromagnets, Nat. Electron. 2, 457 (2019).
  * Cai  _et al._ [2019] X. Cai, T. Song, N. P. Wilson, G. Clark, M. He, X. Zhang, T. Taniguchi, K. Watanabe, W. Yao, D. Xiao, M. A. McGuire, D. H. Cobden, and X. Xu, Atomically thin CrCl3: An in-plane layered antiferromagnetic insulator, Nano Lett. 19, 3993 (2019).
  * Klein  _et al._ [2019] D. R. Klein, D. MacNeill, Q. Song, D. T. Larson, S. Fang, M. Xu, R. A. Ribeiro, P. C. Canfield, E. Kaxiras, R. Comin, and P. Jarillo-Herrero, Enhancement of interlayer exchange in an ultrathin two-dimensional magnet, Nat. Phys. 15, 1255 (2019).
  * Bedoya-Pinto  _et al._ [2021] A. Bedoya-Pinto, J.-R. Ji, A. K. Pandeya, P. Gargiani, M. Valvidares, P. Sessi, J. M. Taylor, F. Radu, K. Chang, and S. S. P. Parkin, Intrinsic 2D-XY ferromagnetism in a van der Waals monolayer, Science 374, 616 (2021).
  * Ghazaryan  _et al._ [2018] D. Ghazaryan, M. T. Greenaway, Z. Wang, V. H. Guarochico-Moreira, I. J. Vera-Marun, J. Yin, Y. Liao, S. V. Morozov, O. Kristanovski, A. I. Lichtenstein, M. I. Katsnelson, F. Withers, A. Mishchenko, L. Eaves, A. K. Geim, K. S. Novoselov, and A. Misra, Magnon-assisted tunnelling in van der Waals heterostructures based on CrBr3, Nat. Electron. 1, 344 (2018).
  * Vizner Stern  _et al._ [2021] M. Vizner Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. Ben Shalom, Interfacial ferroelectricity by van der Waals sliding, Science 372, 1462 (2021).
  * Woods  _et al._ [2021] C. R. Woods, P. Ares, H. Nevison-Andrews, M. J. Holwill, R. Fabregas, F. Guinea, A. K. Geim, K. S. Novoselov, N. R. Walet, and L. Fumagalli, Charge-polarized interfacial superlattices in marginally twisted hexagonal boron nitride, Nat. Commun. 12, 347 (2021).
  * Moore  _et al._ [2021] S. L. Moore, C. J. Ciccarino, D. Halbertal, L. J. McGilly, N. R. Finney, K. Yao, Y. Shao, G. Ni, A. Sternbach, E. J. Telford, B. S. Kim, S. E. Rossi, K. Watanabe, T. Taniguchi, A. N. Pasupathy, C. R. Dean, J. Hone, P. J. Schuck, P. Narang, and D. N. Basov, Nanoscale lattice dynamics in hexagonal boron nitride moiré superlattices, Nat. Commun. 12, 5741 (2021).
  * Weston  _et al._ [2022] A. Weston, E. G. Castanon, V. Enaldiev, F. Ferreira, S. Bhattacharjee, S. Xu, H. Corte-León, Z. Wu, N. Clark, A. Summerfield, T. Hashimoto, Y. Gao, W. Wang, M. Hamer, H. Read, L. Fumagalli, A. V. Kretinin, S. J. Haigh, O. Kazakova, A. K. Geim, V. I. Fal’ko, and R. Gorbachev, Interfacial ferroelectricity in marginally twisted 2D semiconductors, Nat. Nanotech. 17, 390 (2022).
  * Wang  _et al._ [2022] X. Wang, K. Yasuda, Y. Zhang, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, L. Fu, and P. Jarillo-Herrero, Interfacial ferroelectricity in rhombohedral-stacked bilayer transition metal dichalcogenides, Nat. Nanotech. 17, 367 (2022).
  * Liu  _et al._ [2022] Y. Liu, S. Liu, B. Li, W. J. Yoo, and J. Hone, Identifying the transition order in an artificial ferroelectric van der Waals heterostructure, Nano Lett. 22, 1265 (2022).
  * Marmolejo-Tejada  _et al._ [2022] J. M. Marmolejo-Tejada, J. E. Roll, S. P. Poudel, S. Barraza-Lopez, and M. A. Mosquera, Slippery paraelectric transition-metal dichalcogenide bilayers, Nano Lett. 22, 7984 (2022).
  * Sivadas  _et al._ [2018] N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, and D. Xiao, Stacking-dependent magnetism in bilayer CrI3, Nano Lett. 18, 7658 (2018).
  * Wang and Sanyal [2021] D. Wang and B. Sanyal, Systematic study of monolayer to trilayer CrI3: Stacking sequence dependence of electronic structure and magnetism, J. Phys. Chem. C 125, 18467 (2021).
  * Qiu  _et al._ [2021] Z. Qiu, M. Holwill, T. Olsen, P. Lyu, J. Li, H. Fang, H. Yang, K. S. Kashchenko, Mikhailand Novoselov, and J. Lu, Visualizing atomic structure and magnetism of 2D magnetic insulators via tunneling through graphene, Nat. Commun. 12, 70 (2021).
  * McGuire  _et al._ [2015] M. A. McGuire, H. Dixit, V. R. Cooper, and B. C. Sales, Coupling of crystal structure and magnetism in the layered, ferromagnetic insulator CrI3, Chem. Mater. 27, 612 (2015).
  * Li and Wu [2017] L. Li and M. Wu, Binary compound bilayer and multilayer with vertical polarizations: Two-dimensional ferroelectrics, multiferroics, and nanogenerators, ACS Nano 11, 6382 (2017).
  * Yang  _et al._ [2018] Q. Yang, M. Wu, and J. Li, Origin of two-dimensional vertical ferroelectricity in WTe2 bilayer and multilayer, J. Phys. Chem. Lett. 9, 7160 (2018).
  * Zhou  _et al._ [2015] S. Zhou, J. Han, S. Dai, J. Sun, and D. J. Srolovitz, van der Waals bilayer energetics: Generalized stacking-fault energy of graphene, boron nitride, and graphene/boron nitride bilayers, Phys. Rev. B 92, 155438 (2015).
  * Wang  _et al._ [2023] C. Wang, L. You, D. Cobden, and J. Wang, Towards two-dimensional van der Waals ferroelectrics, Nat. Mater. 10.1038/s41563-022-01422-y (2023).
  * Niu  _et al._ [2022] R. Niu, Z. Li, X. Han, Z. Qu, D. Ding, Z. Wang, Q. Liu, T. Liu, C. Han, K. Watanabe, T. Taniguchi, M. Wu, Q. Ren, X. Wang, J. Hong, J. Mao, Z. Han, K. Liu, Z. Gan, and J. Lu, Giant ferroelectric polarization in a bilayer graphene heterostructure, Nat. Commun. 13, 6241 (2022).
  * Zhong  _et al._ [2023] T. Zhong, L. Cheng, Y. Ren, and M. Wu, Theoretical studies of sliding ferroelectricity, magnetoelectric couplings, and piezo-multiferroicity in two-dimensional magnetic materials, Chem. Phys. Lett. 818, 140430 (2023).
  * Wu and Li [2021] M. Wu and J. Li, Sliding ferroelectricity in 2D van der Waals materials: Related physics and future opportunities, Proc. Natl. Acad. Sci. (USA) 118, e2115703118 (2021).
  * Chang  _et al._ [2020] K. Chang, F. Küster, B. J. Miller, J.-R. Ji, J.-L. Zhang, P. Sessi, S. Barraza-Lopez, and S. S. P. Parkin, Microscopic manipulation of ferroelectric domains in SnSe monolayers at room temperature, Nano Lett. 20, 6590 (2020).
  * Rivera [1994] J.-P. Rivera, On definitions, units, measurements, tensor forms of the linear magnetoelectric effect and on a new dynamic method applied to Cr-Cl boracite, Ferroelectrics 161, 165 (1994).