# Fully Differentiable and Interpretable Model for VIO
with 4 Trainable Parameters

Zexi Chen†, Haozhe Du†, Yiyi Liao, Yue Wang, Rong Xiong †: Both authors
contribute equally.All authors are with the State Key Laboratory of Industrial
Control Technology and Institute of Cyber-Systems and Control, Zhejiang
University, Zhejiang, P.R.China. Yue Wang is the corresponding author
<EMAIL_ADDRESS>
###### Abstract

Monocular visual-inertial odometry (VIO) is a critical problem in robotics and
autonomous driving. Traditional methods solve this problem based on filtering
or optimization. While being fully interpretable, they rely on manual
interference and empirical parameter tuning. On the other hand, learning-based
approaches allow for end-to-end training but require a large number of
training data to learn millions of parameters. However, the non-interpretable
and heavy models hinder the generalization ability. In this paper, we propose
a fully differentiable, interpretable, and lightweight monocular VIO model
that contains only 4 trainable parameters. Specifically, we first adopt
Unscented Kalman Filter as a differentiable layer to predict the pitch and
roll, where the covariance matrices of noise are learned to filter out the
noise of the IMU raw data. Second, the refined pitch and roll are adopted to
retrieve a gravity-aligned BEV image of each frame using differentiable camera
projection. Finally, a differentiable pose estimator is utilized to estimate
the remaining 4 DoF poses between the BEV frames. Our method allows for
learning the covariance matrices end-to-end supervised by the pose estimation
loss, demonstrating superior performance to empirical baselines. Experimental
results on synthetic and real-world datasets demonstrate that our simple
approach is competitive with state-of-the-art methods and generalizes well on
unseen scenes.

##  I Introduction

Visual-inertial odometry (VIO) is a fundamental component of the visual-
inertial simultaneous localization and mapping system with many applications
in robotics and autonomous driving. As a cheap and efficient solution,
monocular VIO has attracted growing interest recently. However, monocular VIO
is a challenging task considering the limited sensory information.

Traditional monocular VIO approaches have demonstrated promising results based
on classical feature detection and feature matching [20, 10, 22, 3], see Fig.
1 (top). These methods are interpretable and generalize well, but require
manual interference and empirical parameter tuning. Moreover, these methods
only leverage sparse features, resulting in the degraded performance in
textureless regions [9].

More recently, learning-based approaches address these limitations by learning
the ego-motion in an end-to-end fashion from the raw, dense monocular images,
as illustrated in Fig. 1 (middle). However, they lack interpretability and
thus pose a new challenge to generalization. While many works improve the
interpretability by predicting an intermediate optical flow [19, 18, 5, 2] or
depth map [25, 31, 28, 15], the following pose regression network remains
uninterpretable. Furthermore, existing learning-based approaches rely on a
large number of training images to learn millions of parameters while the ego-
motion has essentially 6 DoF only. We thus ask the following question: Is it
possible to obtain a fully interpretable model with a minimal set of
parameters being learned end-to-end?

Figure 1: Comparison of VIO methods. Top: the standard visual-inertial
odometry framework that requires manual tuning parameters, middle: most
learning based VIO approaches learns depth or optical flow prediction before a
pose regression, and are loaded with tons of heavy parameters for training,
and bottom: our differential interpretable model based learning approach that
solve the problem in an end to end manner, with only 4 trainable parameters.

In this paper, we aim to combine the advantages of the traditional methods
with those of the learning-based approaches. As illustrated in Fig. 1
(bottom), our key idea is to design a fully differentiable, interpretable
model that allows for replacing empirical parameter tuning with learning from
the data. Additionally, our method exploits dense appearance information
instead of sparse features. As it is difficult to differentiate the 6 DoF pose
estimation that involves robust and non-linear optimization, we focus on the 4
DoF pose estimation conditioned on the estimated pitch and roll angles. This
approximation usually works in driving scenarios, e.g., urban roads and
highways, and has gained traction recently [21].

More formally, we consider the monocular VIO as a function that outputs a
relative pose taking as input i) two consecutive frames, and ii) IMU data.
First, we utilize a Differentiable Unscented Kalman Filter (DUKF) to filter
out noise from the IMU raw measurements and retrieve pitch and roll. As the
covariance matrices of noise in both the motion model and the measurement
model are not known precisely in practice, classical UKF relies on empirical
parameter tuning given a specific scene. In contrast, we exploit the fact that
UKF is differentiable and optimize the covariance matrices in an end-to-end
manner leveraging the training data. This requires the estimator of the
remaining 4 DoF pose to be differentiable as well. Therefore, in the second
step, we estimate the remaining 4 DoF poses by converting monocular images to
bird-eye view conditioned on the denoised pitch and roll [21]. The bird-eye
view representation allows us to obtain a global optimal solution using
Differentiable Phase Correlation (DPC) [4] in a single forward pass. Note that
DPC leverages dense appearance information and does not introduce any
trainable parameters. More importantly, it allows for back-propagating the
gradient to the DUKF module to update the covariance matrices. The combination
of DUKF and DPC provides us a fully differentiable, interpretable, and
lightweight model for VIO, named BEVO, with 4 parameters of the covariance
matrices as the only trainable parameters in total. We summarize our
contributions as follows:

  * •

We present a fully differentiable, interpretable, and lightweight monocular
VIO framework that combines the advantages of both traditional and learning-
based approaches.

  * •

We propose to use the Unscented Kalman Filter as a differentiable layer and
leverage the Differentiable Phase Correlation on BEV images for pose
estimation. This novel combination allows for learning the covariance matrices
from data in an end-to-end manner.

  * •

Extensive experiments on KITTI, CARLA, and AeroGround demonstrate that,
surprisingly, our simple approach shows competitive results compared to the
state-of-the-art.

We will make our code publicly available here.

##  II Related Works

In this section, we will introduce the related works of VIO, mainly divided
into two parts, traditional methods and learning-based methods.

###  II-A Traditional Methods

Visual-inertial odometry aims to fuse data from the camera and inertial
measurement unit to estimate the ego-motion. Traditional VIO methods are
mainly based on filtering and optimization. Mourikis et al. [17] propose a
Multi-State Constraint Kalman Filter (MSCKF) method that utilizes the EKF to
estimate poses. Moreover, Li et al. [12] improve the MSCKF approach by
ensuring the correct observability properties and performing online estimation
of calibration parameters. Sun et al. [23] present a stereo version MSCKF
which is robust and efficient. OKVIS [11] optimize through keyframe while
VINS-Mono [20] is a state estimator based on nonlinear optimization, which
contains a tightly coupled visual-inertial odometry and performs global pose
graph optimization. These robust methods can generalize well but require
empirical parameter tuning which is labor intensive.

###  II-B Learning-based Methods

VINet [6] is the first end-to-end learning-based method for visual-inertial
odometry which eliminates the need for manual synchronization and calibration.
DeepVO [26] uses Recurrent Convolutional Neural Networks to learn feature
representation in visual odometry problems. Wang et al. [27] present TartanVO,
which can generalize to multiple datasets and real-world scenarios. DeepVIO
[9] merges 2D optical flow features and IMU data to provide absolute
trajectory estimation, during which the depth and dense point cloud are
estimated. More recent works, e.g., SelfVIO [1], CodeVIO [34], UnDeepVO [13],
Li et al. [14], also take advantage of depth estimation to achieve high pose
estimation accuracy. However, all methods above train a large network with
millions of parameters, resulting in heavy models and are merely interpretable
with weak generalization ability. Therefore, we set to solve this problem by
introducing a fully interpretable model with only 4 trainable parameters.

Figure 2: The overall pipeline of the proposed BEVO. As a demonstration, the
frequency of the front image is 10Hz while the frequency of IMU is 100Hz. For
each frame of IMU, the pitch rate, roll rate, and 3D acceleration are fed into
the differentiable UKF to learn noise filtering (elaborated in the left part
of the figure). However, due to the frequency difference, only frames with
front images have back-propagated gradients (red dashed lines) for the DUKF
training.

##  III Fully Differentiable and Interpretable Model for VIO

We aim to estimate the odometry with the input of one monocular camera and IMU
data by a differentiable and interpretable model. Following the pipeline in
Fig. 2, we employ a trainable UKF to integrate IMU data for pitch and roll
estimation. Then the denoised pitch and roll, together with the corresponding
front image are fed into the second differentiable module, for differentiable
pose estimation of the remaining 4 DoF pose. To this goal, we transform the
image to the BEV images conditioned on the denoised pitch and roll. Following
the BEV projection, the $\mathbb{SIM}(2)$ pose of two consecutive BEV images
is estimated by the DPC. Note that the UKF is a continuous forward function
and is differentiable in essence. In this work, we emphasize its backward path
by referring to it as DUKF, to distinguish it from the traditional UKF
utilization which focuses only on the forward path. As a result, the only
trainable parameters are the 4 numbers that make up the covariance matrices of
motion and measurement in the DUKF and are trained end-to-end with the pose
loss altogether.

###  III-A Differentiable Unscented Kalman Filter

For the DUKF, there are two basic models in the pipeline, namely, the motion
model and measurement model.

Motion Model: Given a consecutive series of pitch rate and roll rate from the
IMU raw data from time $t-1$ to $t$, and the filtered pitch and roll in time
$t-1$, we have the state in time $t-1$ as:

$x_{t-1}=\begin{bmatrix}\alpha_{t-1},&\beta_{t-1}\end{bmatrix}^{T},$ (1)

where $\alpha$, and $\beta$ are the pitch and roll angle respectively.
Following the sampling theory of UKF [24], the sampled state in $t-1$ becomes:

$X_{t-1}=\Gamma(x_{t-1}),$ (2)

where $\Gamma$ is the sampling rule. Given the motion model $G$, we have the
initial guess of pitch and roll in time $t$ as:

$\bar{X}_{t}^{*}=G(X_{t-1})=\begin{bmatrix}\alpha_{t-1}+\omega_{\alpha_{t-1}},&\beta_{t-1}+\omega_{\beta_{t-1}}\end{bmatrix}^{T}$
(3) $\displaystyle\bar{\mu}_{t}=\sum{w_{m}\bar{X}_{t}^{*}}$ (4)
$\displaystyle\bar{\sigma}_{t}=\sum{w_{c}(\bar{X}_{t}^{*}-\bar{\mu}_{t})(\bar{X}_{t}^{*}-\bar{\mu}_{t})^{T}}+\boldsymbol{O_{t}},$
(5)

where $\omega_{\alpha}$, and $\omega_{\beta}$ are the pitch and roll rate from
the gyroscope, and $\bar{\mu}_{t}$ and $\bar{\sigma}_{t}$ are the mean and
variance of the initial guess of $\bar{X}_{t}^{*}$, respectively. Note that
$w_{m}$ and $w_{c}$ are the weights for each sampled state in
$\bar{X}_{t}^{*}$ based on the theory in [24], and $\boldsymbol{O_{t}}$ is the
covariance matrix for the noise of the motion model to be trained. Since
$\boldsymbol{O_{t}}$ is a $2\times 2$ diagonal matrix, we parameterize it with
the two parameters on the diagonal.

Measurement Model: After updating the covariance matrix $\boldsymbol{O_{t}}$,
we retrieve a new sampled state $\bar{X}_{t}$ in time $t$:

$\displaystyle\bar{x}_{t}\sim
N(\bar{\mu}_{t},\bar{\sigma}_{t})\triangleq\begin{bmatrix}\bar{\alpha}_{t},&\bar{\beta}_{t}\end{bmatrix}^{T}$
(6) $\displaystyle\bar{X}_{t}=\Gamma(\bar{x}_{t}),$ (7)

where $x_{t}$ satisfies the normal distribution. Given the measurement model
$H$, the initial guess of measurement is

$\displaystyle\bar{Z}_{t}=H(\bar{X}_{t})=\begin{bmatrix}\bar{\alpha}_{t},&\bar{\beta}_{t}\end{bmatrix}^{T}$
(8) $\displaystyle\bar{M}_{t}=\sum{w_{m}\bar{Z}_{t}}$ (9)
$\displaystyle\bar{\Sigma}_{t}=\sum{w_{c}(\bar{Z}_{t}-\bar{M}_{t})(\bar{Z}_{t}-\bar{M}_{t})^{T}}+\boldsymbol{Q_{t}},$
(10)

where $\bar{M}_{t}$ and $\bar{\Sigma}_{t}$ are the mean and variance of
$\bar{X}_{t}$ respectively, and $\boldsymbol{Q_{t}}$ is the covariance matrix
for the noise of the measurement model to be trained and is parametrized
following the same rule as $\boldsymbol{O_{t}}$.

Now, considering a real-world measurement is obtained:

$\displaystyle\hat{\alpha_{t}}=arctan(acc^{x}_{t}/{acc^{z}_{t}})$ (11)
$\displaystyle\hat{\beta_{t}}=-arctan(acc^{y}_{t}/\sqrt{{acc^{x}_{t}}^{2}+{acc^{z}_{t}}^{2}})$
(12) $\displaystyle
Z_{t}=\begin{bmatrix}\hat{\alpha_{t}},\hat{\beta_{t}}\end{bmatrix}^{T},$ (13)

where $\hat{\alpha_{t}}$ and $\hat{\beta_{t}}$ are the measurement of pitch
and roll in time $t$ with respect to the raw acceleration data in each axes
$acc^{x}_{t}$, $acc^{y}_{t}$, and $acc^{z}_{t}$. Altogether, the real-world
measurement state is formed as $Z_{t}$.

By now, we come to the last stage of the DUKF: update the final state $x_{t}$.

$\displaystyle\bar{\Sigma}^{X,Z}_{t}=\sum\limits_{i=0}^{2n}{w_{c}}^{i}(\bar{X}^{i}_{t}-\bar{\mu}_{t})(\bar{Z}^{i}_{t}-\bar{M_{t}})^{T}$
(14) $\displaystyle K_{t}=\bar{\Sigma}^{X,Z}_{t}{\bar{\Sigma}^{-1}_{t}}$ (15)
$\displaystyle\mu_{t}=\bar{\mu}_{t}+K_{t}(Z_{t}-\bar{M_{t}})$ (16)
$\displaystyle\sigma_{t}=\bar{\sigma}_{t}+K_{t}\bar{\Sigma}_{t}{K_{t}}^{T}$
(17) $\displaystyle x_{t}\sim
N(\mu_{t},\sigma_{t})\triangleq\begin{bmatrix}\alpha_{t},&\beta_{t}\end{bmatrix}^{T}$
(18)

where $\bar{\Sigma}^{X,Z}_{t}$, $K_{t}$ are the cross-covariance and Kalman
Gain respectively, $\mu_{t}$ and $\sigma_{t}$ are the mean and varience of
$x_{t}$ which will also guide the sampling in the following time $t+1$.

###  III-B Differentiable Bird-eye View Projection

To achieve the projection from the front view to the birds-eye view, we first
define the transformation matrix from the IMU to the world frame considering
its pitch, roll, and elevation:

$\displaystyle R^{IMU}_{t}=\begin{bmatrix}cos\beta_{t}&sin\beta_{t}&0\\\
-cos\alpha_{t}sin\beta_{t}&cos\alpha_{t}cos\beta_{t}&sin\alpha_{t}\\\
sin\alpha_{t}sin\beta_{t}&-sin\alpha_{t}cos\beta_{t}&cos\alpha_{t}\end{bmatrix}$
(19) $\displaystyle T^{IMU}_{t}=\begin{bmatrix}0&0&const_{z}\end{bmatrix},$
(20)

the rotation matrix $R^{IMU}_{t}$ of the IMU at time $t$ changes constantly
with $\alpha_{t}$ and $\beta_{t}$ which are estimated by the DUKF. The
translation vector $T^{IMU}_{t}$ of the IMU is known and fixed since it is the
urban road and highway that we are considering in which the height of the IMU
to the ground $const_{z}$ stays unchanged through the journey [21].
Transformation matrices from camera to the world $[R^{CAM}_{t},T^{CAM}_{t}]$
are obtianed with extrinsic calibration of $[R^{IMU}_{t},T^{IMU}_{t}]$.

Besides the unchanged height, one more assumption can be made for a safe trip:
there should be a valuable space in front of the vehicle, forming that the
center bottom part of the front image should represent a plane ground, shown
as the red dashed box in Fig. 2.

With the assumption above, we consider pixels within the red dashed box as a
set of points in the 3D space, and map them from the camera plane to the
ground plane.

To begin with, we generate a set of even points $P_{g}$ on the ground plane
that stand for the location that the BEV should be projected:

$P_{g}\triangleq\\{X^{{}^{\prime}},Y^{{}^{\prime}},Z^{{}^{\prime}}\\}^{T},$
(21)

where $X^{{}^{\prime}}$,$Y^{{}^{\prime}}$, and $Z^{{}^{\prime}}$ are the
coordinates of points. $P_{g}$ is then transformed to the camera coordinate
with:

$P^{t}_{c}={R^{CAM}_{t}}^{T}(P_{g}-T^{CAM}_{t})\triangleq\\{X^{{}^{\prime\prime}}_{t},Y^{{}^{\prime\prime}}_{t},Z^{{}^{\prime\prime}}_{t}\\}^{T},$
(22)

and $P^{t}_{c}$ is projected to the camera plane with the focal lengths
$f_{x}$, $f_{y}$ and principal point $(c_{x},c_{y})$:

$Z^{{}^{\prime\prime}}_{t}\begin{bmatrix}u^{{}^{\prime}}_{t}\\\
v^{{}^{\prime}}_{t}\\\ 1\end{bmatrix}=\begin{bmatrix}f_{x}&0&c_{x}\\\
0&f_{y}&c_{y}\\\ 0&0&1\end{bmatrix}\begin{bmatrix}X^{{}^{\prime\prime}}_{t}\\\
Y^{{}^{\prime\prime}}_{t}\\\ Z^{{}^{\prime\prime}}_{t}\end{bmatrix},$ (23)

where $(u^{{}^{\prime}}_{t},v^{{}^{\prime}}_{t})$ are the indexes on the front
camera image of each corresponding point in $P_{g}$ in time $t$. Therefore, by
remaping the original pixels $(u,v)$ to
$(u^{{}^{\prime}}_{t},v^{{}^{\prime}}_{t})$, we can further retrieve the BEV
projection $I^{bev}_{t}$ out of the red dashed box $I^{f}_{t}$ in Fig. 2.

$I^{bev}_{t}=F^{(u^{{}^{\prime}}_{t},v^{{}^{\prime}}_{t})}_{(u,v)}(I^{f}_{t}),$
(24)

$F^{(u^{{}^{\prime}}_{t},v^{{}^{\prime}}_{t})}_{(u,v)}(\cdot)$ is the
remapping of $\cdot$ from $(u,v)$ to
$(u^{{}^{\prime}}_{t},v^{{}^{\prime}}_{t})$. Since it is parameterized by the
pitch and roll, which is related to the variance of the IMU data, with the
help of DUKF, the gradient can be back-propagated to the variance.

###  III-C Differentiable Phase Correlation

Given two images $I^{bev}_{t-1}$ and $I^{bev}_{t}$ with pose transformations,
a variation of the differentiable phase correlation [4] is utilized to
estimate the the overall relative pose $\xi_{i}$ between $I^{bev}_{t-1}$ and
$I^{bev}_{t}$. The rotation and scale is calculated by

$\displaystyle
p(\xi^{\theta,s}_{i})=\mathfrak{C}(\mathfrak{L}(\mathfrak{F}(I^{bev}_{t})),\mathfrak{L}(\mathfrak{F}(I^{bev}_{t-1})))$
(25)
$\displaystyle(\theta_{t},s_{t})=\mathop{\mathbb{E}}(p(\xi^{\theta,s}_{i})),$
(26)

where $\xi^{\theta,s}_{t}$ is the rotation and scale part of $\xi_{i}$ between
$I^{bev}_{t-1}$ and $I^{bev}_{t}$, $\mathfrak{F}$ is the discrete Fourier
Transform, $\mathfrak{L}$ is the log-polar transform, and $\mathfrak{C}$ is
the phase correlation solver. Fourier Transformation $\mathfrak{F}$ transforms
images into the Fourier frequency domain of which the magnitude has the
property of translational insensitivity, therefore the rotation and scale are
decoupled with displacements and are represented in the magnitude. Log-polar
transformation $\mathfrak{L}$ transforms Cartesian coordinates into log-polar
coordinates so that such rotation and scale in the magnitude of Fourier domain
are remapped into displacement in the new coordinates, making it solvable with
the afterward phase correlation solver. Phase correlation solver
$\mathfrak{C}$ outputs a heatmap indicating the displacements of the two log-
polar images, which eventually stands for the rotation $\theta_{t}$ and scale
$s_{t}$ of the two input images $I^{bev}_{t-1}$ and $I^{bev}_{t}$. To make the
solver differentiable, we use expectation $\mathop{\mathbb{E}}(\cdot)$ as the
estimation of $\cdot$.

Then $I^{bev}_{t}$ is rotated and scaled referring to $\xi^{\theta,s}_{t}$
with the result of $\bar{I}^{bev}_{t}$:

$\bar{I}^{bev}_{t}=\begin{bmatrix}s_{t}\boldsymbol{R}_{\theta_{t}}&0\\\
0&1\end{bmatrix}I^{bev}_{t}.$ (27)

In the same manner, translations $\xi^{\textbf{t}}_{t}$ between
$\bar{I}^{bev}_{t}$ and $I^{bev}_{t-1}$ is calculated

$\displaystyle
p(\xi^{\textbf{t}}_{t})=\mathfrak{C}(\bar{I}^{bev}_{t},I^{bev}_{t-1})$ (28)
$\displaystyle\textbf{t}_{t}=\mathop{\mathbb{E}}(p(\xi^{\textbf{t}}_{t})).$
(29)

With the relative rotation $\theta_{t}$, scale $s_{t}$, and translation
$\textbf{t}_{t}$ between $I^{bev}_{t-1}$ and $I^{bev}_{t}$, and by consulting
with the fixed $Distance\;per\;Pixel$ (DDP), we finally arrive at the
$\mathbb{SIM}(2)$ pose between frame $t-1$ and $t$ of vehicle:

$\boldsymbol{S}_{t}=\begin{bmatrix}\theta_{t},&s_{t},&DDP\cdot\textbf{t}_{t}\end{bmatrix}.$
(30)

Note that the gradient in the part can be back-propagated all the way back to
the variance of the IMU data since the phase correlation is differentiable.

###  III-D Loss Design, Back-Propagation, and Training

With the elaboration above, we define a one-step VIO as a function:

$T=\mathcal{VIO}_{O,Q}(I_{t-1},I_{t},\omega,a)$ (31)

where $\mathcal{VIO}_{O,Q}$ is fully differentiable and interpretable with 4
trainable parameters in $O_{t}$ and $Q_{t}$. $I$, $\omega$, and $a$ are the
front images, gyroscope data, and accelerometer data. Since the network is
fully differentiable, it is trained end to end and supervised only with the
ground truth of the odometry between frames. Such supervision is applied to
optimize the covariance matrices of noise in both the motion model and the
measurement model of DUKF.

Elaborately, we wish to supervise the distribution of the poses
$\xi^{\theta,s}_{t}$ and $\xi^{\textbf{t}}_{t}$ and therefore, we calculate
the Kullback–Leibler Divergence loss $KLD$ between
$\\{p(\xi^{\theta,s}_{t}),p(\xi^{\textbf{t}}_{t})\\}$ and the ground truth:

$\displaystyle\mathcal{L}=KLD(p(\xi^{\theta,s}_{t}),\boldsymbol{1}_{\theta_{GT}})+KLD(p(\xi^{\textbf{t}}_{t}),\boldsymbol{1}_{\textbf{t}_{GT}})$
(32)

where $\boldsymbol{1}_{\theta_{GT}}$ and $\boldsymbol{1}_{\textbf{t}_{GT}}$
are the Gaussian Blurred one-peak distributions centered at the ground truth
$\theta_{GT}$ and $\textbf{t}_{GT}$ respectively.

Training Details: all trainings are conducted on a GPU server with two NVIDIA
RTX 3090 as GPUS and Intel i9 9900k as the CPU. For the optimizer, we choose
the ADAM. We only supervise the model with 4 DoF ground truth (GT) data of the
same frequency as the images, i.e. the high frequentcy GT of the IMU is not
needed.

Plug into Map-Based Localization System: in addition to the odometry, we also
explore the possibility of integrating BEVO with fully differentiable
localization to achieve a drift-free system. However, this is not the focus of
our work and therefore we will not elaborate upon it, but show result as
demonstration of the application.

Figure 3: Top: The qualitative results of the evaluation on KITTI sequence 09,
and 10. Bottom: Noise filtering of the pitch and roll angle through sequence
10. It indicates that even though pitch and roll are not directly supervised,
the low-frequency image-level supervision in 4 DoF only still allows the model
to obtain good results in pitch and roll.

##  IV Experimental Setup

We conduct two main experiments: i) odometry estimation on the KITTI [8]
odometry dataset on which the experiments on ablation study, testing,
generalization are conducted, and ii) heterogeneous localization on CARLA [7]
and the AeroGround Datast [32].

###  IV-A Dataset and Corresponding Experiments

KITTI: We evaluate BEVO on the KITTI odometry dataset using sequences 00-08
for training and 09-10 for the test following the same train/validate/test
split as [16, 14]. In addition to the widely adopted train/validate split
where $20\%$ samples are randomly sampled from the training set for
validation, we also train our model in the first three quarters of each
sequence of 00-08, and evaluate the performance in the last quarter. This
train/validate split is more realistic even though it might present a worse
result. Note that sequence 03 cannot be used for training or evaluation in
this study since the raw data where we grab the IMU data is missing.

CARLA & AG: We further verify the plausibility of applying the BEVO as a
front-end plugin for end-to-end heterogeneous localization following the
routine in [29]. Such experiments are conducted in the simulator CARLA and on
the real-world AeroGround Dataset.

  * •

CARLA: In CARLA, we localize a vehicle on the satellite map in different
weathers. We train and utilize BEVO as the odometry and localize the projected
BEV (from different weathers) on the heterogeneous satellite map.

  * •

AeroGround: The AeroGround (AG) Dataset is collected for multi-robot
collaboration. In this dataset, we train and utilize BEVO as the odometry and
localize the ground robot with its front camera BEV on the heterogeneous map
built by a drone.

###  IV-B Evaluation Metrics

We evaluate our trajectories using $t_{rel}$ and $r_{rel}$ from the standard
KITTI toolkit as in [14]. In our case, we do not estimate the elevation due to
the assumption that the vehicle will not leave the ground. However, we still
recorded the initial elevation and keep it constant through the journey for
comparison, which might degrade the performance in $t_{rel}$ to some extent.

###  IV-C Baselines

For the odometry experiment, we compare to both traditional method VINS[20]
and SOTA depth dependent learning-based methods DeepVIO[9], Li et al. [14],
GeoNet[30], SelfVIO[1], and depth independent learning-based method
TartanVO[27]. For the localization experiments, we compare the performance
with VINS-Mono [20].

TABLE I: Quantitative results of odometry estimation on KITTI dataset from
sequence 00 to 10. Red indicates the best performance and Blue indicates the
second.

| DeepVIO | TartanVO | Li et al. | GeoNet | SelfVIO | VINS | Ours1 | Ours2 | Ours3
---|---|---|---|---|---|---|---|---|---
Seq | frames | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$
00 | 4541 | 11.62 | 2.45 | $\setminus$ | $\setminus$ | 14.21 | 5.93 | 44.08 | 14.89 | 1.24 | 0.45 | 18.83 | 2.49 | 3.26 | 1.44 | 1.55 | 0.83 | 2.01 | 1.02
01 | 1101 | $\setminus$ | $\setminus$ | $\setminus$ | $\setminus$ | 21.36 | 4.62 | 43.21 | 8.42 | $\setminus$ | $\setminus$ | $\setminus$ | $\setminus$ | 5.62 | 3.95 | 1.17 | 1.66 | 1.76 | 2.13
02 | 4661 | 4.52 | 1.44 | $\setminus$ | $\setminus$ | 16.21 | 2.60 | 73.59 | 12.53 | 0.80 | 0.25 | 21.03 | 2.61 | 4.56 | 1.32 | 1.02 | 0.88 | 1.79 | 1.55
04 | 271 | $\setminus$ | $\setminus$ | $\setminus$ | $\setminus$ | 9.08 | 4.41 | 17.91 | 9.95 | $\setminus$ | $\setminus$ | $\setminus$ | $\setminus$ | 3.79 | 1.05 | 2.04 | 0.92 | 3.12 | 1.15
05 | 2761 | 2.86 | 2.32 | $\setminus$ | $\setminus$ | 24.82 | 6.33 | 32.47 | 13.12 | 0.89 | 0.63 | 21.90 | 2.72 | 2.05 | 1.57 | 0.72 | 1.01 | 1.55 | 1.80
06 | 1101 | $\setminus$ | $\setminus$ | 4.72 | 2.95 | 9.77 | 3.58 | 40.28 | 16.68 | $\setminus$ | $\setminus$ | $\setminus$ | $\setminus$ | 3.86 | 2.08 | 1.07 | 1.56 | 2.77 | 1.99
07 | 1101 | 2.71 | 1.66 | 4.32 | 3.41 | 12.85 | 2.30 | 37.13 | 17.20 | 0.91 | 0.49 | 15.39 | 2.42 | 3.58 | 1.47 | 0.88 | 1.35 | 1.15 | 1.79
08 | 4701 | 2.13 | 1.02 | $\setminus$ | $\setminus$ | 27.10 | 7.81 | 11.45 | 62.45 | 1.09 | 0.36 | 32.66 | 3.09 | 1.94 | 1.04 | 0.93 | 0.92 | 1.29 | 1.00
09 | 1591 | 1.38 | 1.12 | 6.0 | 3.11 | 15.21 | 5.28 | 13.02 | 67.06 | 1.95 | 1.15 | 41.47 | 2.41 | 1.22 | 1.05 | 1.29 | 1.12 | 1.32 | 1.31
10 | 1201 | 0.85 | 1.03 | 6.89 | 2.73 | 25.63 | 7.69 | 58.52 | 23.02 | 1.81 | 1.30 | 20.35 | 2.73 | 1.01 | 1.01 | 1.14 | 1.03 | 1.17 | 1.30
Ave 00-10 |  | 3.73 | 1.57 | 5.48 | 3.05 | 16.02 | 4.60 | 33.78 | 22.30 | 1.24 | 0.67 | 21.61 | 2.64 | 2.80 | 1.45 | 1.07 | 1.02 | 1.63 | 1.36
Ave 09-10 |  | 1.12 | 1.08 | 6.74 | 2.92 | 20.42 | 6.49 | 35.77 | 45.04 | 1.88 | 1.23 | 30.91 | 2.57 | 1.12 | 1.03 | 1.22 | 1.07 | 1.24 | 1.31

TABLE II: Quantitative results of ablation study conducted on KITTI. All of
the methods are trained in sequence 00-08 following the train spit in [14] and
tested in sequence 09-10.

| w/o $\boldsymbol{O}_{t}$ | w/o $\boldsymbol{Q}_{t}$ | w/o DUKF | BEVO | w GT
---|---|---|---|---|---
Seq | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$
09 | 4.82 | 4.73 | 5.90 | 5.11 | 22.05 | 71.58 | 1.29 | 1.12 | 1.05 | 1.04
10 | 5.66 | 4.87 | 6.26 | 5.62 | 25.64 | 68.54 | 1.14 | 1.03 | 0.96 | 1.01

##  V Experimental Results

###  V-A Ablation Study

We conduct several ablation studies on KITTI to validate the effectiveness of
the proposed DUKF as well as the learned covariance matrices.

Several experiments conducted in KITTI are designed for the ablation study to
reassure the important role a DUKF plays in the BEVO as well as the role
trainable covariance matrices play in DUKF. First, we replace the learned
covariance matrix $\boldsymbol{O}_{t}$ by one that is empirically set in [24].
The quantitative result in Table. II shows that the accuracy drops drastically
when using a manually tuned matrix, suggesting that the empirical parameter
tuning leads to inaccurate pitch and roll estimation, and subsequently
deteriorates the BEV generation and odometry estimation. Nest, we restore
$\boldsymbol{O}_{t}$ with the learned matrix and replace $\boldsymbol{Q}_{t}$
with an empirical covariance matrix as [24], and note it as “w/o
$\boldsymbol{Q}_{t}$”. The quantitative result also verifies that trainable
covariance matrices for noise help to improve the performance of UKF. However,
even though the performances are degraded with empirical parameters, they are
still competitive due to the dense matching of DPC who seeks global optima.

We further study the importance of DUKF on the BEVO by eliminating the whole
DUKF and fed the raw pitch and roll (directly calculated from the
accelerometer) to the BEV projection. We note it as “w/o DUKF”. The
quantitative result shown in Table. II proves that when the filter is gone,
the BEVO will encounter serious tracking failure in bumpy roads. One more
ablation study is conducted by replacing the output pitch and roll of DUKF
with the ground truth, denoted as “w GT”. The result that BEVO’s performance
is comparable with “w GT” indicates that by learning the covariance matrices,
the error brought by DUKF is too small to influence the odometry than the
error by brought by DPC, showing the effectiveness of the training.

TABLE III: Quantitative results of localization. Note: in CARLA, we localize
the vehicle in different weathers against a sunny satellite map, CARLA-S means
the vehicle runs in a sunny weather, CARLA-F means the vehicle is in a foggy
weather, and CARLA-N means the vehicle runs at night.

| BEVO | VINS-Mono | BEVO+
---|---|---|---
Dataset | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$ | $t_{rel}$ | $r_{rel}$
CARLA-S | 1.05 | 1.82 | 2.06 | 2.58 | 0.21 | 0.33
CARLA-F | 1.11 | 1.38 | 2.92 | 2.14 | 0.24 | 0.31
CARLA-N | 0.98 | 1.03 | 1.82 | 1.47 | 0.19 | 0.28
AG | 1.39 | 1.72 | 2.55 | 2.41 | 0.30 | 0.46

###  V-B Odometry Estimation

We now compare the odometry estimation performance of the proposed BEVO to the
state-of-the-art methods in two training settings from section IV-A: i)
following train/validate splits in [16][14] (Ours1), and ii) split each
sequence by quarters, take the first three quarters of each sequence to train
and the last quarter to validate (Ours2). In addition to these two settings,
we further train BEVO using the first $20\%$ of the training data of each
sequence in Ours2 to evaluate its performance on insufficient data (Ours3).

Fig.3 and Table I show the qualitative and quantitative comparisons,
respectively. Table I shows that when trained in the same settings as other
baselines, Ours2 outperforms most of the baselines in validation results in
sequence 00 to 08. However, regarding the unseen testing sequences, our Ours1
beats all other baseline in both $t_{rel}$ and $r_{rel}$ except for the
$t_{rel}$ of DeepVIO in sequence 10. This may be explained by the fact that
DeepVIO utilizes stereo cameras instead of monocular camera considered in our
method. The result of the testing sequence shows that our method is better at
generalization compared to other learning-based methods.

It is worth addressing that with only 4 trainable parameters and without
employing any nonlinear optimization such as sliding window, our method is
surprisingly comparable with others on sequences 00-08, and is even slightly
better in testing on sequences 09 and 10. The result of Ours3 suggests that,
in contrast to learning-based approaches reliant on a large number of training
images, our method is more data-efficient thanks to the lightweight parameters
and the fully interpretable model.

Figure 4: The demonstration of localizing a front camera’s BEV of a ground
moving robot in drone’s map.

###  V-C Plugin for Heterogeneous Localization

In this part, we show the possibility of enabling an end-to-end localization
with BEVO as the odometry following our previous work [29]. Elaborations on
the implementation details of the fully differentiable drift-free localization
can be found in the supplementary material [33]. By combining BEVO with the
heterogeneous localization RaLL (denoted as BEVO+), the whole pipeline is
fully differentiable, thus allowing for more accurate pose estimation.

We conduct experiments in CARLA and the real-world AG dataset and compare the
localization with BEVO itself as well as the VINS-Mono. Results in Fig. 4 and
Table II proves that the BEVO is able to achieve good performance even when it
acts as a plugin in an end-to-end trainable localization framework. It
indicates a promising future for upcoming end-to-end autonomous driving, e.g.,
auto weeder.

##  VI Conclusion

We present a fully differentiable, interpretable and lightweight model for
monocular VIO, namely, BEVO. Our model comprises a UKF to denoise pitch and
roll via trainable covariance matrices, a BEV projection, and a differentiable
phase correlation to estimate the $\mathbb{SIM}(2)$ pose on the BEV image. The
full pipeline is differentiable, allowing us to train the covariance matrices
end-to-end and thus avoiding empirical parameter tuning. In various
experiments, BEVO presents competitive performance in accuracy and
generalization. We believe our interpretable, simple approach provides an
alternative perspective of traditional methods and should be considered as a
baseline for future works on learning-based monocular VIO.

## References

  * [1] Yasin Almalioglu et al.  “Selfvio: Self-supervised deep monocular visual-inertial odometry and depth estimation”  In _arXiv preprint arXiv:1911.09968_ , 2019
  * [2] Xicheng Ban et al.  “Monocular visual odometry based on depth and optical flow using deep learning”  In _IEEE Transactions on Instrumentation and Measurement_ 70 IEEE, 2020, pp. 1–19
  * [3] Martin Brossard, Silvere Bonnabel and Axel Barrau  “Unscented Kalman filter on Lie groups for visual inertial odometry”  In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 649–655  IEEE
  * [4] Zexi Chen, Xuecheng Xu, Yue Wang and Rong Xiong  “Deep Phase Correlation for End-to-End Heterogeneous Sensor Measurements Matching”  In _arXiv preprint arXiv:2008.09474_ , 2020
  * [5] Cheng Chuanqi, Hao Xiangyang, Zhang Zhenjie and Zhao Mandan  “Monocular visual odometry based on optical flow and feature matching”  In _2017 29th Chinese Control And Decision Conference (CCDC)_ , 2017, pp. 4554–4558  IEEE
  * [6] Ronald Clark et al.  “Vinet: Visual-inertial odometry as a sequence-to-sequence learning problem”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 31.1, 2017
  * [7] Alexey Dosovitskiy et al.  “CARLA: An open urban driving simulator”  In _Conference on robot learning_ , 2017, pp. 1–16  PMLR
  * [8] Andreas Geiger, Philip Lenz, Christoph Stiller and Raquel Urtasun  “Vision meets robotics: The kitti dataset”  In _The International Journal of Robotics Research_ 32.11 Sage Publications Sage UK: London, England, 2013, pp. 1231–1237
  * [9] Liming Han, Yimin Lin, Guoguang Du and Shiguo Lian  “Deepvio: Self-supervised deep learning of monocular visual inertial odometry using 3d geometric constraints”  In _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 6906–6913  IEEE
  * [10] Yijia He et al.  “Pl-vio: Tightly-coupled monocular visual–inertial odometry using point and line features”  In _Sensors_ 18.4 Multidisciplinary Digital Publishing Institute, 2018, pp. 1159
  * [11] Stefan Leutenegger et al.  “Keyframe-based visual–inertial odometry using nonlinear optimization”  In _The International Journal of Robotics Research_ 34.3 SAGE Publications Sage UK: London, England, 2015, pp. 314–334
  * [12] Mingyang Li and Anastasios I Mourikis  “High-precision, consistent EKF-based visual-inertial odometry”  In _The International Journal of Robotics Research_ 32.6 Sage Publications Sage UK: London, England, 2013, pp. 690–711
  * [13] Ruihao Li, Sen Wang, Zhiqiang Long and Dongbing Gu  “Undeepvo: Monocular visual odometry through unsupervised deep learning”  In _2018 IEEE international conference on robotics and automation (ICRA)_ , 2018, pp. 7286–7291  IEEE
  * [14] Shunkai Li et al.  “Self-supervised deep visual odometry with online adaptation”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 6339–6348
  * [15] Xiangyu Li et al.  “Dvonet: unsupervised monocular depth estimation and visual odometry”  In _2019 IEEE Visual Communications and Image Processing (VCIP)_ , 2019, pp. 1–4  IEEE
  * [16] Reza Mahjourian, Martin Wicke and Anelia Angelova  “Unsupervised learning of depth and ego-motion from monocular video using 3d geometric constraints”  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 5667–5675
  * [17] Anastasios I Mourikis and Stergios I Roumeliotis  “A multi-state constraint Kalman filter for vision-aided inertial navigation”  In _Proceedings 2007 IEEE International Conference on Robotics and Automation_ , 2007, pp. 3565–3572  IEEE
  * [18] Peter Muller and Andreas Savakis  “Flowdometry: An optical flow and deep learning based approach to visual odometry”  In _2017 IEEE Winter Conference on Applications of Computer Vision (WACV)_ , 2017, pp. 624–631  IEEE
  * [19] Tejas Pandey, Dexmont Pena, Jonathan Byrne and David Moloney  “Leveraging deep learning for visual odometry using optical flow”  In _Sensors_ 21.4 Multidisciplinary Digital Publishing Institute, 2021, pp. 1313
  * [20] Tong Qin, Peiliang Li and Shaojie Shen  “Vins-mono: A robust and versatile monocular visual-inertial state estimator”  In _IEEE Transactions on Robotics_ 34.4 IEEE, 2018, pp. 1004–1020
  * [21] Tong Qin, Tongqing Chen, Yilun Chen and Qing Su  “Avp-slam: Semantic visual mapping and localization for autonomous vehicles in the parking lot”  In _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020, pp. 5939–5945  IEEE
  * [22] Ethan Rublee, Vincent Rabaud, Kurt Konolige and Gary Bradski  “ORB: An efficient alternative to SIFT or SURF”  In _2011 International conference on computer vision_ , 2011, pp. 2564–2571  Ieee
  * [23] Ke Sun et al.  “Robust stereo visual inertial odometry for fast autonomous flight”  In _IEEE Robotics and Automation Letters_ 3.2 IEEE, 2018, pp. 965–972
  * [24] Eric A Wan, Rudolph Van Der Merwe and Simon Haykin  “The unscented Kalman filter”  In _Kalman filtering and neural networks_ 5.2007 Wiley New York, 2001, pp. 221–280
  * [25] Rui Wang, Stephen M Pizer and Jan-Michael Frahm  “Recurrent neural network for (un-) supervised learning of monocular video visual odometry and depth”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 5555–5564
  * [26] Sen Wang, Ronald Clark, Hongkai Wen and Niki Trigoni  “Deepvo: Towards end-to-end visual odometry with deep recurrent convolutional neural networks”  In _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , 2017, pp. 2043–2050  IEEE
  * [27] Wenshan Wang, Yaoyu Hu and Sebastian Scherer  “TartanVO: A Generalizable Learning-based VO”  In _arXiv preprint arXiv:2011.00359_ , 2020
  * [28] Mingkang Xiong et al.  “Self-supervised Monocular Depth and Visual Odometry Learning with Scale-consistent Geometric Constraints.”  In _IJCAI_ , 2020, pp. 963–969
  * [29] Huan Yin, Runjian Chen, Yue Wang and Rong Xiong  “Rall: end-to-end radar localization on lidar map using differentiable measurement model”  In _IEEE Transactions on Intelligent Transportation Systems_ IEEE, 2021
  * [30] Zhichao Yin and Jianping Shi  “Geonet: Unsupervised learning of dense depth, optical flow and camera pose”  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 1983–1992
  * [31] Huangying Zhan et al.  “Unsupervised learning of monocular depth estimation and visual odometry with deep feature reconstruction”  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 340–349
  * [32] ZJU-Robotics-Lab  “AeroGround Dataset”, 2020  URL: https://github.com/ZJU-Robotics-Lab/OpenDataSet
  * [33] ZJU-Robotics-Lab  “Supplementary Material”, 2020  URL: https://github.com/jessychen1016/BEVO_pytorch
  * [34] Xingxing Zuo et al.  “CodeVIO: Visual-Inertial Odometry with Learned Optimizable Dense Depth”  In _arXiv preprint arXiv:2012.10133_ , 2020


# A scalable network model for electrically tunable ferroelectric domain
structure in twistronic bilayers of two-dimensional semiconductors

V. V. Enaldiev  University of Manchester, School of Physics and Astronomy,
Oxford Road, Manchester M13 9PL, United Kingdom  National Graphene Institute,
University of Manchester, Booth St. E. Manchester M13 9PL, United Kingdom
Kotelnikov Institute of Radio-engineering and Electronics of the RAS,
Mokhovaya 11-7, Moscow 125009, Russia  F. Ferreira  University of Manchester,
School of Physics and Astronomy, Oxford Road, Manchester M13 9PL, United
Kingdom  National Graphene Institute, University of Manchester, Booth St. E.
Manchester M13 9PL, United Kingdom  V.I. Fal’ko  University of Manchester,
School of Physics and Astronomy, Oxford Road, Manchester M13 9PL, United
Kingdom  National Graphene Institute, University of Manchester, Booth St. E.
Manchester M13 9PL, United Kingdom  Henry Royce Institute, University of
Manchester, Booth St. E. Manchester M13 9PL, United Kingdom

###### Abstract

Moiré structures in small-angle-twisted bilayers of two-dimensional (2D)
semiconductors with a broken-symmetry interface form arrays of ferroelectric
(FE) domains with periodically alternating out-of-plane polarization. Here, we
propose a network theory for the tunability of such FE domain structure by
applying an electric field perpendicular to the 2D crystal. Using multiscale
analysis, we derive a fully parametrized string-theory-like description of the
domain wall network (DWN) and show that it undergoes a qualitative change,
after the arcs of partial dislocation (PD) like domain walls merge (near the
network nodes) into streaks of perfect screw dislocations (PSD), which happens
at a threshold displacement field dependent on the DWN period.

Two-dimensional material twistronics, fueled by discoveries of new phenomena
in twisted graphene bilayers [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] and trilayers
[11, 12, 13, 14, 15, 16], has recently expanded onto a broader range of van
der Waals systems [17, 18, 19, 20, 21, 22, 23]. In general, twistronic
structures are associated with geometrical moiré patterns: a periodic
variation local stacking of the two layers. In long-period moiré patterns,
characteristic for small-angle-twisted bilayers, the areas of energetically
preferential stacking expand into mesoscale domains [24, 25, 26, 27], embedded
into a domain wall network (DWN). In particular, by assembling a homobilayer
of two inversion-asymmetric honeycomb monolayers (hBN or transition metal
dichalcogenides (TMD)) with parallel orientation of their unit cells, one
obtains a triangular array of domains with broken mirror and inversion
symmetries [28, 29, 30, 31, 32, 33] and an out-of-plane ferroelectric (FE)
polarization [34, 35, 36], which direction alternates between the neighboring
XM-stacking domains (metal atoms overlaying chalcogens) and their MX twins
(chalcogens over metal atoms) [37, 38].

Figure 1:  Triangular XM/MX DWN (inset) in twisted TMD bilayers with FE
domains areas varied by an out-of-plane displacement field, $D$. At $D\geq
D_{*}$, XM/MX domain walls merge into streaks of perfect screw dislocations
separating - near the network nodes - domains with the same FE polarization.
The remaining arcs of XM/MX domain walls have a universal shape, scaled with
$D_{*}/D$ ratio.

The coupling between the out-of-plane FE polarization and an externally
controlled displacement field, $D$, varies the energies of XM and MX, changing
the ratio between their areas and leading to the deformation of DWN. Here, we
offer a generic theory for the field-tunable FE domain structure in twistronic
TMD bilayers, fully quantified with the help of multiscale modelling approach
[38, 39, 40, 36] for MX2-bilayers (M=Mo,W; X=S,Se). For weak fields
$|D|<D_{*}$, the continuously deforming domain walls retain their partial
dislocation character, however, above the threshold, $|D|/D_{*}\geq 1$, pairs
of partial dislocations (PD) form streaks of perfect (full) screw dislocations
(PSD, of length $\mathcal{L}$) near the network nodes. In Fig. 1, we
illustrate how the arcs of XM/MX domain walls with a universal shape merge and
split apart upon the variation of $D$, where the threshold field and DWN
parameters scale with elastic parameters of a partial dislocation and the DWN
period, $\ell$, as $D_{*}\propto 1/\ell$.

The scenario of the DWN transformation in Fig. 1 is a result of the following
analysis. First, we use density functional theory (DFT) to quantify the
coupling of FE polarization of an asymmetric TMD interface to an external out-
of-plane displacement field, $D$. By taking into account the resulting
coupling in the competition between the inter-layer adhesion and elastic
strain in the layers, which has been used [38] and tested [27] earlier in the
studies of mesoscale lattice relaxation in TMD bilayers, we derive an
effective theory, formulated in terms of the DWN deformations. Finally, we
find an analytical solution for such a ’string-like’ theory, which has a
universal form scaling with the $D/D_{*}$ ratio.

Figure 2:  Top: Solid (dashed) line shows potential energy of electrons
computed in DFT for a double supercell structure of TMD bilayers with 3R
stacking (here, MoS2) with and without external displacement field. A map,
superimposed on the lattice structure, visualises FE charge density,
$\delta\rho$ (averaged over the unit cell are). The two twin bilayers in the
expanded supercell, used in the DFT modelling, represent the actual lattice
structure of MX and XM domains in the reconstructed moire superlattice of a
TMD bilayer, which at $D=0$ have a triangular shape (middle panel). Bottom:
DFT-computed 3R- and 2H-bilayers energy per unit cell as a function of a
displacement field. Solid lines are plotted using Eq.(1) with parameters
listed in Table 1.

Ab initio modelling of FE bilayers in the out-of-plane electric field. Here,
we use two methods to carry out DFT calculations for 3R-bilayers, with the
coinciding results. In the first method, we use Quantum ESPRESSO (QE) [41, 42]
to construct a supercell with a pair of mirror-reflected (MX and XM) bilayers,
separated by a large vacuum gap, as shown in Fig. 2. This choice of the
structure eliminates an issue with periodic boundary condition for potential
of the FE charges [36]. In the second method we construct a single 3R-bilayer
and use Coulomb truncation in the out-of-plane direction [43] with QE and
dipole-correction with VASP[44]. An example of the computed charge transfer,
$\delta\rho$, between the layers of an individual bilayer (which integral,
$\int z\delta\rho(z)dz\equiv P$, determines the areal density of the FE dipole
moment) and a potential drop, $\Delta$, across the double charge layer, is
shown in Fig. 2. To mention, the computed values of $P$ and $\Delta$ are
related as $P=\epsilon_{0}\Delta$, in agreement with the earlier studies of 3D
bulk FE materials [45].

The displacement field, accounted for by adding a triangular potential,
$-D|z|/\epsilon_{0}$, to the input pseudo-potentials, produces a shift,
$\delta U\equiv U(D)-U(0)$ of the bilayer energy. The latter includes a linear
term, accounting for the FE polarization, and a quadratic downward shift due
to the dielectric out-of-plane polarizability of the material:

$\delta U=-\frac{PD}{\chi\epsilon_{0}}-\frac{\alpha_{zz}^{\rm
3R}D^{2}}{2\epsilon_{0}\mathcal{A}}.$ (1)

Here, we parametrise the $P-D$ coupling by a dimensionless parameter $\chi$
and the out-of-plane dielectric polarizability by $\alpha_{zz}^{\rm 3R}$
(where $\mathcal{A}$ is the unit cell area of a TMD monolayer). The DFT-
computed energies of bilayers of four different TMDs are shown in Fig. 2: in
this computation, we used two different supercell periods ($L_{z}=$40 Å and
$L_{z}=$70 Å, which set the length of the vacuum spacer). Using Eq. (1), we
find the polarizability values, which are very close to polarizability,
$\alpha^{2H}_{zz}$, computed for the inversion-symmetric 2H bilayer111For
monolayers and 2H-bilayers we used the Coulomb truncation in the out-of-plane
direction[43].. When recalculated per monolayer, these values also agree with
the separately computed monolayer polarizability, $\alpha_{zz}$
($\alpha_{zz}=\frac{1}{2}\alpha^{3R}_{zz}=\frac{1}{2}\alpha^{2H}_{zz}$, see in
Table 1). This confirms that the dielectric response of a wide band gap van
der Waals materials is determined by the intra-layer polarization of the
constituent atoms 222 Unlike Refs. [55, 56], which use $z$-averaged electric
field for computation of dielectric permitivities of bilayers, we express the
quadratic amendment to the total energy via out-of-plane displacement field,
conserving across every cross-section of the structure. This allows us to
avoid uncertainties in $\alpha_{zz}^{\rm 3R}$, which may appear at averaging
of the electric field in crystals with a few out-of-plane unit cells [45]. We
find that the polarizability, computed for a TMD monolayer, its 2H and 3R
bilayer, and thicker 2H films, linearly scales with with the number of layers.
This observation contradicts some earlier DFT studies of dielectric
susceptibility of TMDs [55, 56] which claimed a pronounce layer-number-
dependence, but agrees with the more recent results [49] published by some of
the authors of Refs [55]. . These values enabled us to estimate the $z$-axis
dielectric susceptibility of bulk TMD crystals [48, 49, 50],
$\epsilon_{zz}=(1-\frac{\alpha_{zz}}{\mathcal{A}d})^{-1}$, arriving at the
values in the range of $\epsilon_{zz}\sim 6-7.5$, listed in Table 1.

Despite a substantial polarizability of monolayers, the data in Fig. 2 are
described well by Eq. (1) with $\chi\approx 1$, pointing towards a decoupling
of the inter-layer FE charge transfer from the intra-layer dielectric
polarizability. Moreover, by comparing the DFT-computed values of the double-
layer potential drop $\Delta$ to the FE coupling with the displacement field,
$D$, in Eq. (2), we find that the linear in $D$ energy shift in MX and XM
domains (which have opposite out-of-plane FE polarization, $\pm P$, and
voltage drop across the double charge layer, $\pm\Delta$) can be described
very well as 333The expression for interaction energy of the FE polarization
with displacement field naturally comes when assuming a local dielectric
permittivity in a continuum medium approximation. Indeed, suppose the FE
charges, with plane-averaged density $\delta\rho(z)$, are placed in the medium
with local dielectric permittivity $\epsilon_{zz}(z)$. From the Poisson
equation and electro-neutrality condition
$\int_{-\infty}^{+\infty}\delta\rho(z)dz=0$
$\partial_{z}(\epsilon_{zz}(z)\partial_{z}\varphi(z))=-\delta\rho(z)/\epsilon_{0}$,
we express the potential drop across the layer of charges as
$\Delta=\int_{-\infty}^{+\infty}\partial_{z}\varphi(z)dz=-\int_{-\infty}^{+\infty}dz\int_{-\infty}^{z}dz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}=\int_{-\infty}^{+\infty}dz\int^{+\infty}_{z}dz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}=\int\int_{z<z^{\prime}}dzdz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}$.
At the same time, interaction energy of these charges with uniform external
out-of-plane displacement field (related to local electric field as
$D=\epsilon_{0}\epsilon_{zz}(z)E(z)$) reads as $\delta
U=-\int_{-\infty}^{+\infty}dz\delta\rho(z)\int_{-\infty}^{z}dz^{\prime}D/\epsilon_{zz}(z^{\prime})\epsilon_{0}=-D\int_{-\infty}^{+\infty}dz^{\prime}\int_{-\infty}^{z^{\prime}}dz\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}=-D\int\int_{z<z^{\prime}}dzdz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}\equiv-D\Delta$.
After changing variables, $z\leftrightarrow z^{\prime}$, at the last step of
that calculation, we arrive at the relation in Eq. (2).

$\delta U_{\rm MX}=-\delta U_{\rm XM}\approx-D\Delta.$ (2)

Below, we will use this coupling to model tunability of the DWN by
displacement field.

Table 1: Dielectric and FE parameters for various TMDs. Left: polarizability
($\alpha_{zz}$) (used to estimate the dielectric permittivity $\epsilon_{zz}$
of 3D-bulk TMDs), and FE coupling parameter ($\chi$) for various 3R-TMD
bilayers$*)$. Right: FE potential drop parameters in Eq. (A scalable network
model for electrically tunable ferroelectric domain structure in twistronic
bilayers of two-dimensional semiconductors).

| $\alpha^{\rm 3R}_{zz}$ | $\alpha^{\rm 2H}_{zz}$ | $\alpha_{zz}$ | $\epsilon_{zz}$ | $\chi$ |  $\Delta$, |  $\Delta_{a}$, |  $q$,
---|---|---|---|---|---|---|---|---
| Å3 | Å3 | Å3 |  |  | mV | mV | nm-1
MoS2 | 89.89 | 89.80 | 44.46 | 6.45 | 1.03 | 68 | 16.4 | 22.152
89.93 | 89.85 | 44.48 | 6.48 | 1.02
WS2 | 88.43 | 88.37 | 43.95 | 5.95 | 1.00 | 62 | 15 | 22.598
88.43 | 88.37 | 43.94 | 5.95 | 1.00
MoSe2 | 105.32 | 105.22 | 52.29 | 7.65 | 1.05 | 66 | 15.7 | 20.520
105.31 | 105.25 | 52.24 | 7.67 | 1.07
WSe2 | 104.37 | 104.33 | 52.03 | 7.39 | 1.04 | 65 | 15.3 | 20.953
104.32 | 104.28 | 51.97 | 7.32 | 1.07

  * $*)$

For $\alpha^{\rm 3R}_{zz}$, $\chi$ and $\epsilon_{zz}$ we show the values
obtained with QE (top row) and - for comparison - with VASP (bottom row). For
QE we used full-relativistic ultra-soft pseudo-potentials and a plane-wave
cutoff energy of 70 Ry. For VASP, we used PAW full-relativistic pseudo-
potentials and 60 Ry cut-off (with a $13\times 13\times 1$ $k$-point grid and
a Perdew-Burke-Ernzerhof (PBE)[52] exchange-correlation functional for both QE
and VASP).

Mesoscale model for lattice reconstruction is formulated in terms of the
bilayer energy dependence on local stacking of the two layers and their strain
[27, 38, 39], incorporated via an interlayer offset,
$\bm{r}_{0}\left(\bm{r}\right)=\theta\hat{z}\times\bm{r}+\bm{u}^{(t)}-\bm{u}^{(b)}$,
which varies across the superlattice, as prescribed by a small-angle,
$\theta$, misalignment between the top (t) and bottom (b) crystals and their
elastic deformations, $\bm{u}^{(t/b)}$. Locally, stacking determines the
interlayer distance, which corresponds to the minimum of the stacking-
dependent adhesion energy, $\mathcal{W}(\bm{r}_{0},d)$, quantified for the
four TMDs using DFT modelling and displayed in Fig. 3. For convenience, the
inter-layer distance ($d$) dependence of the computed
$\mathcal{W}(\bm{r}_{0},d)$ is shown as a function of $Z=d-d_{0}$, counted
from the minimum (at $d=d_{0}$) of the offset-averaged adhesion energy, which
coincides with $\mathcal{W}(\bm{r}_{0}=[-a/3,0],d)$ (one of stacking
configurations shown in Fig. 3). As a result, energy density, characterizing
the relaxation functional, reads:

$\displaystyle\mathcal{E}=$
$\displaystyle\sum_{a=t,b}\left[\frac{1}{2}\lambda\left(u_{ii}^{(a)}\right)^{2}+\mu
u_{ij}^{(a)}u_{ji}^{(a)}\right]$
$\displaystyle+\mathcal{W}\left[\bm{r}_{0}(\bm{r}),Z\right]-D\Delta\left[\bm{r}_{0}(\bm{r}),Z\right];$
(3) $\displaystyle\mathcal{W}(\bm{r}_{0},Z)$ $\displaystyle=\kappa
Z^{2}+\left(\gamma-\gamma^{\prime}Z\right)\sum_{n=1,2,3}{\cos\left(\bm{G}_{n}\cdot\bm{r}_{0}\right)}.$

Here, the first term accounts for strain,
$u_{ij}^{(t/b)}=\frac{1}{2}\left(\partial_{i}u_{j}^{(t/b)}+\partial_{j}u_{i}^{(t/b)}\right)$
($\lambda$ and $\mu$ are the monolayer elastic moduli). The second term
describes adhesion energy between top and bottom layers, where $\kappa$
determines curvature of the adhesion energy of 3R-stacked bilayers, which
characterises frequency of layer breathing mode [38], and $\bm{G}_{1,2,3}$ is
the first star of reciprocal lattice vectors of TMD monolayer,
$\left|\bm{G}_{1,2,3}\right|=4\pi/a\sqrt{3}$ (for details, see Refs. [27, 38,
39]). The last term is responsible for the energy shift due to the external
field.

The next steps in the analysis require choosing an optimal interlayer distance
between the two monolayer, which depends on the encapsulation environment of
the bilayer. For a bilayer in vacuum, or encapsulated into a soft flexible
matrix with a weaker adhesion to TMD than the the interlayer adhesion
$\cal{W}$, we find from Eq. (A scalable network model for electrically tunable
ferroelectric domain structure in twistronic bilayers of two-dimensional
semiconductors) that, locally,
$Z(\bm{r}_{0})=\frac{\gamma^{\prime}}{2\kappa}\sum_{n}\cos(\bm{G}_{n}\cdot\bm{r}_{0})$.
For a strong substrate-TMD coupling, the monolayers would remain flat, and the
optimal interlayer distance would be $Z=0$ across the entire moiŕe pattern.
For either of these two cases, we substitute the corresponding choice of $Z$
in Eq. (A scalable network model for electrically tunable ferroelectric domain
structure in twistronic bilayers of two-dimensional semiconductors), and also
into the local voltage drop across the double layer of charge (in the last
term in Eq. (A scalable network model for electrically tunable ferroelectric
domain structure in twistronic bilayers of two-dimensional semiconductors))
corresponding to the local stacking configuration, which form has been
established earlier [36]:

$\Delta\left[\bm{r}_{0},Z(\bm{r}_{0})\right]=\Delta_{a}e^{-qZ\left(\bm{r}_{0}\right)}\sum_{n=1,2,3}\sin\left(\bm{G}_{n}\cdot\bm{r}_{0}\right).$

For XM and MX stackings, this gives $\Delta\equiv\Delta(\bm{r}_{0}^{\rm
MX})=-\Delta(\bm{r}_{0}^{\rm XM})$, which values, taken from the recent ab
initio simulations [40, 38], are listed in Table 2.

Table 2: Elastic moduli, adhesion and partial dislocation energy density parameters for the studied TMD bilayers.  |  $\gamma$, |  $\gamma^{\prime}$, |  $\kappa$, |  $\mu$ , |  $\lambda$, |  $a$, |  $d_{0}$, | $\bar{w}$ | $\widetilde{w}$ | $u$
---|---|---|---|---|---|---|---|---|---|---
| eV/nm2 | eV/nm3 | eV/nm4 | N/m | N/m | nm | nm | eV/nm | eV/nm | eV/nm
MoS2 | 0.189 | 5.634 | 214 | 70.9 | 83.2 | 0.316 | 0.636 | 0.96 | 0.69 | 2.24
WS2 | 0.212 | 6.327 | 213 | 72.5 | 52.5 | 0.315 | 0.638 | 1.05 | 0.69 | 2.45
MoSe2 | 0.1975 | 5.707 | 189 | 49.6 | 42.3 | 0.329 | 0.670 | 0.79 | 0.54 | 1.84
WSe2 | 0.1514 | 4.381 | 190 | 48.4 | 29.7 | 0.328 | 0.671 | 0.74 | 0.46 | 1.69
Figure 3:  (a-d) TMD bilayer adhesion energies for various offsets and
interlayer distances, $Z=d-d_{0}$, counted from the optimal distance, $d_{0}$,
of stacking configurations averaged adhesion energy. $Z_{\rm 3R}$ labels
interlayer distance used to evaluate the FE parameters in 3R (MX and XM) TMD
bilayers. Details of DFT data (square dots) are described in Ref. [38]. Insets
show stackings with corresponding in-plane offsets between the layers. (e)
Orientation dependence of partial dislocation energy per unit length
calculated in Ref. [38] and fitted with $\bar{w}+\widetilde{w}\sin^{2}\phi$,
where values of parameters are listed in Table 2. $\phi=0$ corresponds to the
armchair direction. (f) Evolution of domain wall shape with displacement
field, given by solution of Eq. (5) for $D<D_{*}$ and $D\geq D_{*}$.

Network model. While the minimization of the energy functional in (A scalable
network model for electrically tunable ferroelectric domain structure in
twistronic bilayers of two-dimensional semiconductors) enables one to find all
mesoscale details of the lattice adjustments of the two crystalline plane of
the bilayer, it is more practical to use another method for the analysis of
small-angle ($\theta\lesssim 0.5^{\circ}$) twisted bilayers. This is because,
in the latter case, most of the areas of the moiré pattern is occupied by the
homogeneous MX and XM stacking domains with a characteristic size of
$\ell\gtrsim 50$ nm, whereas the deformations are concentrated inside narrow
domain walls, which are only only few nanometres in width [38, 27]. Then, we
define the DWN energy,

$\mathcal{E}_{\ell}\equiv\int_{\rm
supercell}d^{2}\bm{r}\left\\{\mathcal{E}-\mathcal{W}\left[\bm{r}_{0}^{\rm
MX},Z\left(\bm{r}_{0}^{\rm MX}\right)\right]\right\\},$

by subtracting the energy of a uniform MX-stacked bilayer at $D=0$ from
$\cal{E}$ in Eq. (A scalable network model for electrically tunable
ferroelectric domain structure in twistronic bilayers of two-dimensional
semiconductors). This reduces the problem to the analysis of
$\mathcal{E}_{\ell}$, where we can treat each domain wall as a string,
characterised by energy per unit length dependent on a crystallographic
orientation of the domain wall axis [38], with - as shown in Fig. 3(e) - a
pronounced minimum at the armchair direction in the TMD crystal.

A finite displacement field acts as an external drive for increasing areas of
domains with the preferential ferroelectric polarization. This leads to the
bending of the partial dislocation (PD) domain walls, hence, increasing their
energy due to their elongation and the PD axis deviation from the armchair
axis. The exact form on the deformed DWN can be found by solving a ’string
theory’ model, expressed in terms of a deflection, $y(x)$, of PD segments from
the closest armchair direction, formalised using an energy functional,

$\mathcal{E}_{\ell}[y\left(x\right),\delta]=3\int_{\delta}^{\ell-\delta}\left[\left(\bar{w}+\widetilde{w}\frac{y^{\prime
2}}{1+y^{\prime 2}}\right)\sqrt{1+y^{\prime
2}}-2yD\Delta\right]dx+2\sqrt{3}\left[u-D\Delta\delta\right]\delta.$ (4)

Here, the first term in the integral describes the orientation-dependence of
the PD energy, $\bar{w}+\widetilde{w}\sin^{2}{\phi}$, and its stretching,
accounted by a factor $\sqrt{1+{y^{\prime}}^{2}}$ with $\phi$ standing for an
angle between the dislocation axis and the closest armchair direction in the
crystal, so that $\sin^{2}\phi=y^{\prime 2}/(1+y^{\prime 2})$. Two ’stiffness’
parameters, $\bar{w}$ and $\widetilde{w}$, were determined using the data from
Ref. [38], see Table 2 for their values for various TMDs. The second term in
the integral stands for the energy gain from a bigger area of the
energetically preferential FE polarization domain. The last two terms in
$\mathcal{E}_{\ell}$ accounts for two PDs merging, near each DWN node, into
streaks of PSDs with a length $\mathcal{L}=\delta/\sqrt{3}$ (projected onto
$0\leq x\leq\delta$ and $\delta\leq x\leq\delta-\ell$ intervals), which
energetically preferable orientation is along zigzag axis in the crystal and
energy per unit length is $u$ (Table 2).

Using variational principle, we obtain an equation for the the shape of each
string,

$\begin{split}y^{\prime\prime}&=-\frac{2D\Delta\left(1+y^{\prime
2}\right)^{\frac{5}{2}}}{\bar{w}+2\widetilde{w}+(\bar{w}-\widetilde{w})y^{\prime
2}},\\\ y&(\delta)=y(\ell-\delta)=\frac{\delta}{\sqrt{3}},\end{split}$ (5)

where $\delta$ appears as an additional variable. Note that zero values of
$y(x)$ at the network nodes requires vanishing of derivative in the middle of
interval, i.e. providing a symmetrical shape of the PD, with
$y^{\prime}(\ell/2)=0$.

For a small $D$, energetically favourable domains grow due to bending of PDs,
with their ends fixed at network nodes, $\delta=0$. This elastic stretching
has no threshold in $D$ and the string form is described by 444At small $D$
(when $y^{\prime 2}\ll 1$) approximate solution of Eq. (5) is given by
parabola: $y(x)=2D\Delta(\ell-x)x/(\bar{w}+2\widetilde{w})$.,

$y(x)=\int_{0}^{x}\frac{f(x^{\prime})}{\sqrt{1-f^{2}(x^{\prime})}}dx^{\prime},\quad
0\leq x\leq\ell,$ (6)

where $f(x)$ is a real root of

$f^{3}-\left(\frac{\bar{w}}{\widetilde{w}}+2\right)f-2\frac{D\Delta}{\widetilde{w}}\left(x-\frac{\ell}{2}\right)=0.$
(7)

Note that $f(\ell/2)=0$ and that $f(x)\equiv 0$ for $D=0$.

Bending, described by Eq. (6), takes a new qualitative form after each pair of
PDs near the network nodes ($x=0;\ x=\ell$) touch each other, see Fig. 4(j).
This condition corresponds to $\phi^{\circ}(0)=-\phi^{\circ}(l)=30^{\circ}$
(i.e., $y^{\prime}\left(0\right)=-y^{\prime}\left(\ell\right)=1/\sqrt{3}$),
which, together with Eq. (6) determines a threshold displacement field,

$\displaystyle D_{\ast}$ $\displaystyle=\epsilon_{0}\frac{\mathcal{V}}{\ell},$
(8) $\displaystyle\mathcal{V}=$
$\displaystyle\frac{\left(\frac{\bar{w}}{2}+\frac{7\widetilde{w}}{8}\right)}{\epsilon_{0}\Delta}\Rightarrow\begin{cases}\mathcal{V}_{\rm
MoS_{2}}=287\,{\rm V},\mathcal{V}_{\rm WS_{2}}=310\,{\rm V},\\\
\mathcal{V}_{\rm MoSe_{2}}=236\,{\rm V},\mathcal{V}_{\rm WSe_{2}}=214\,{\rm
V}.\end{cases}$

It is important to note that pairs of touching PDs cannot annihilate, as a sum
of their Burgers vectors $\bm{b}_{1}+\bm{b}_{2}=\bm{b}$
($\bm{b}_{1}=a(1/\sqrt{3},0)$, $\bm{b}_{2}=a(1/2\sqrt{3},1/2)$,
$\bm{b}=a(\sqrt{3}/2,1/2)$) corresponds to the Burgers vector ($|\bm{b}|=a$)
of a PSD, required to maintain an overall twist between two monolayers in the
bilayer. This means that they actually form streaks of PSDs, which lengths
grow with a further increase of the electric drive, as

$\mathcal{L}=\frac{\ell}{\sqrt{3}}\left(1-\frac{D_{*}}{D}\right),$ (9)

The latter scaling law follows from the solution of Eq. (5) on the interval
$\delta\leq x\leq\ell-\delta$ with the boundary conditions,
$y^{\prime}\left(\delta\right)=-y^{\prime}\left(\ell-\delta\right)=1/\sqrt{3}$,
where $\delta=\sqrt{3}\mathcal{L}/2$ is a projection of $\mathcal{L}$ onto the
$x$-axis. Hence, above the threshold, DWN is composed of (A) PSD streaks near
each network node aligned with zigzag directions in TMD, and separating the
adjacent expanded energetically favourable 3R-domains and (B) PD arcs
splitting up at the PSD ends, $x=\delta$ and $x=\ell-\delta$, and touching
each other at the splitting points.

To find the shape of the remaining PD arcs for any $D>D_{*}$, we rewrite Eq.
(5) in dimensionless coordinates, $\widetilde{x}$ and $\widetilde{y}$ defined
by $x=\ell\widetilde{x}D_{\ast}/D+\delta$,
$y=\ell\widetilde{y}D_{\ast}/D+\delta/\sqrt{3}$,

$\displaystyle\widetilde{y}^{\prime\prime}=$
$\displaystyle-\frac{\bar{w}+\frac{7}{4}\widetilde{w}}{\bar{w}+2\widetilde{w}}\frac{\left(1+\widetilde{y}^{\prime
2}\right)^{\frac{5}{2}}}{1+\frac{\bar{w}-\widetilde{w}}{\bar{w}+2\widetilde{w}}\widetilde{y}^{\prime
2}},\quad 0\leq\widetilde{x}\leq 1;$
$\displaystyle\widetilde{y}^{\prime}(0)=-\widetilde{y}^{\prime}(1)=\frac{1}{\sqrt{3}}.$

Its solution has the same structure for any $D\geq D_{*}$,

$\displaystyle\widetilde{y}\left(\widetilde{x}\right)=\int_{0}^{\widetilde{x}}{\frac{\widetilde{f}\left(\widetilde{x}^{\prime}\right)}{\sqrt{1-\widetilde{f}^{2}\left(\widetilde{x}^{\prime}\right)}}d\widetilde{x}^{\prime}},\qquad\widetilde{f}(0)=\frac{1}{2};$
(10)
$\displaystyle\widetilde{f}^{3}-\left(\frac{\bar{w}}{\widetilde{w}}+2\right)\widetilde{f}-\left(\frac{\bar{w}}{\widetilde{w}}+\frac{7}{4}\right)\left(\widetilde{x}-\frac{1}{2}\right)=0,$

describing the universal shape of partial dislocation arcs, scaled by the
ratio $D_{*}/D$, as shown in Fig. 3(f).

Finally, to validate the accuracy of the obtained analytical solution,
obtained using the network model, we also performed the mesoscale lattice
relaxation for a MoS2 bilayer with a $0.1^{\circ}$ interlayer twist
(corresponding to $\ell\approx 180$ nm), by minimising the full energy
functional in Eq. (A scalable network model for electrically tunable
ferroelectric domain structure in twistronic bilayers of two-dimensional
semiconductors). For this, we analysed the Euler-Lagrange equations for the
displacement fields $\bm{u}^{(t/b)}$, obtained by applying variational
principle to $\mathcal{E}\left[\bm{u}^{(t/b)}\right]$. We solved those
numerically, seeking for periodic (with the the period $\ell$) solutions on a
sufficiently dense grid using interior point method implemented in GEKKO
Optimization Suite package [54]. The computed fields were used to obtain the
local values of the inter-layer offset,
$\bm{r}_{0}\left(\bm{r}\right)=\theta\hat{z}\times\bm{r}+\bm{u}^{(t)}-\bm{u}^{(b)}$,
which we substituted into Eqs. (A scalable network model for electrically
tunable ferroelectric domain structure in twistronic bilayers of two-
dimensional semiconductors) to find the interlayer distance, $Z$, and to map
the interlayer potential drop across the domain structure. The results for the
latter quantity are shown in Fig. 4 for various values of displacement field,
displaying a close agreement with the analytical solution of the network model
shown by yellow lines.

Figure 4:  (a-d) Maps of the interlayer potential drop, $\Delta$ in MoS2
bilayer with a $\theta=0.1^{\circ}$ twist, computed by the minimisation of
energy in Eq. (A scalable network model for electrically tunable ferroelectric
domain structure in twistronic bilayers of two-dimensional semiconductors) for
various values of $D$ below and above the threshold. These maps display the
deformation of DWN from a triangular web of partial dislocations - see Figs. 1
and 2, with a $\ell\approx 180$ nm period. Yellow lines are the streaks of PSD
and arcs of PD described by Eqs. (6)-(10).

To conclude, the developed effective network model gives an efficient
description of the long-period domain structure in small-angle-twisted TMD
bilayers with ferroelectric interface and its deformation by an out-of-plane
electric field. Its analytical solution in Eqs. (6) - (10) gives a simple tool
for the interpretation of experimentally measured variations of the domains’
shapes, such as in the recently studied bilayers with ferroelectric properties
[30, 29, 31, 32, 33].

Acknowledgments. We thank Andre Geim, Roman Gorbachev, Philip Kim, and David
Vanderbilt for useful discussions. This study has been supported by the
European Graphene Flagship Core3 Project, EU Quantum Flagship Project 2D-SIPC,
and EPSRC grants P/W006502/1, EP/V007033/1, EP/S030719/1. The reported ab
initio simulations were performed using EPSRC ARCHER facility and the
University of Manchester CSF.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and et al, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018b).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Cao  _et al._ [2021] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, Nematicity and competing orders in superconducting magic-angle graphene, Science 372, 264 (2021).
  * Kerelsky  _et al._ [2021] A. Kerelsky, C. Rubio-Verdú, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wang, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, D. N. Basov, A. Rubio, and A. N. Pasupathy, Moiréless correlations in abca graphene, PNAS 118 (2021).
  * Xu  _et al._ [2019] S. G. Xu, A. I. Berdyugin, P. Kumaravadivel, F. Guinea, R. Krishna Kumar, D. A. Bandurin, S. V. Morozov, W. Kuang, B. Tsim, S. Liu, J. H. Edgar, I. V. Grigorieva, V. I. Fal’ko, M. Kim, and A. K. Geim, Giant oscillations in a triangular network of one-dimensional states in marginally twisted graphene, Nature Communications 10, 4008 (2019).
  * Gadelha  _et al._ [2021] A. C. Gadelha, D. A. A. Ohlberg, C. Rabelo, E. G. S. Neto, T. L. Vasconcelos, J. L. Campos, J. S. Lemos, V. Ornelas, D. Miranda, R. Nadas, F. C. Santana, K. Watanabe, T. Taniguchi, B. van Troeye, M. Lamparski, V. Meunier, V.-H. Nguyen, D. Paszko, J.-C. Charlier, L. C. Campos, L. G. Cançado, G. Medeiros-Ribeiro, and A. Jorio, Localization of lattice dynamics in low-angle twisted bilayer graphene, Nature 590, 405 (2021).
  * Kazmierczak  _et al._ [2021] N. P. Kazmierczak, M. Van Winkle, C. Ophus, K. C. Bustillo, S. Carr, H. G. Brown, J. Ciston, T. Taniguchi, K. Watanabe, and D. K. Bediako, Strain fields in twisted bilayer graphene, Nature Materials 20, 956 (2021).
  * Polshyn  _et al._ [2020] H. Polshyn, J. Zhu, M. A. Kumar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. H. MacDonald, _et al._ , Electrical switching of magnetic order in an orbital chern insulator, Nature 588, 66 (2020).
  * Xu  _et al._ [2021] S. Xu, M. M. Al Ezzi, N. Balakrishnan, A. Garcia-Ruiz, B. Tsim, C. Mullan, J. Barrier, N. Xin, B. A. Piot, T. Taniguchi, K. Watanabe, A. Carvalho, A. Mishchenko, A. K. Geim, V. I. Fal’ko, S. Adam, A. H. C. Neto, K. S. Novoselov, and Y. Shi, Tunable van hove singularities and correlated states in twisted monolayer–bilayer graphene, Nature Physics 17, 619 (2021).
  * Chen  _et al._ [2019] G. Chen, L. Jiang, S. Wu, B. Lyu, H. Li, B. L. Chittari, K. Watanabe, T. Taniguchi, Z. Shi, J. Jung, Y. Zhang, and F. Wang, Evidence of a gate-tunable mott insulator in a trilayer graphene moiré superlattice, Nature Physics 15, 237 (2019).
  * Shen  _et al._ [2020] C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, _et al._ , Correlated states in twisted double bilayer graphene, Nature Physics 16, 520 (2020).
  * Chen  _et al._ [2021] S. Chen, M. He, Y.-H. Zhang, V. Hsieh, Z. Fei, K. Watanabe, T. Taniguchi, D. H. Cobden, X. Xu, C. R. Dean, _et al._ , Electrically tunable correlated and topological states in twisted monolayer–bilayer graphene, Nature Physics 17, 374 (2021).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021).
  * Kunstmann  _et al._ [2018] J. Kunstmann, F. Mooshammer, P. Nagler, A. Chaves, F. Stein, N. Paradiso, G. Plechinger, C. Strunk, C. Schüller, G. Seifert, and et al, Momentum-space indirect interlayer excitons in transition-metal dichalcogenide van der Waals heterostructures, Nature Physics 14, 801 (2018).
  * Rivera  _et al._ [2018] P. Rivera, H. Yu, K. L. Seyler, N. P. Wilson, W. Yao, and X. Xu, Interlayer valley excitons in heterobilayers of transition metal dichalcogenides, Nature Nanotechnology 13, 1004 (2018).
  * Nayak  _et al._ [2017] P. K. Nayak, Y. Horbatenko, S. Ahn, G. Kim, J.-U. Lee, K. Y. Ma, A.-R. Jang, H. Lim, D. Kim, S. Ryu, and et al, Probing evolution of twist-angle-dependent interlayer excitons in MoSe2/WSe2 van der Waals heterostructures, ACS Nano 11, 4041 (2017).
  * Alexeev  _et al._ [2019] E. M. Alexeev, D. A. Ruiz-Tijerina, M. Danovich, M. J. Hamer, D. J. Terry, P. K. Nayak, S. Ahn, S. Pak, J. Lee, J. I. Sohn, and et al, Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures, Nature 567, 81 (2019).
  * McGilly  _et al._ [2020] L. J. McGilly, A. Kerelsky, N. R. Finney, K. Shapovalov, E.-M. Shih, A. Ghiotto, Y. Zeng, S. L. Moore, W. Wu, Y. Bai, K. Watanabe, T. Taniguchi, M. Stengel, L. Zhou, J. Hone, X. Zhu, D. N. Basov, C. Dean, C. E. Dreyer, and A. N. Pasupathy, Visualization of moiré superlattices, Nature Nanotechnology 15, 580 (2020).
  * Zhang  _et al._ [2020] Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc, and B. J. LeRoy, Flat bands in twisted bilayer transition metal dichalcogenides, Nature Physics 16, 1093 (2020).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Correlated electronic phases in twisted bilayer transition metal dichalcogenides, Nature Materials 19, 861 (2020).
  * Alden  _et al._ [2013] J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. Brown, J. Park, D. A. Muller, and P. L. McEuen, Strain solitons and topological defects in bilayer graphene, PNAS 110, 11256 (2013).
  * Yoo  _et al._ [2019] H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, and et al, Atomic and electronic reconstruction at the van der Waals interface in twisted bilayer graphene, Nature materials 18, 448 (2019).
  * Rosenberger  _et al._ [2020] M. R. Rosenberger, H.-J. Chuang, M. Phillips, V. P. Oleshko, K. M. McCreary, S. V. Sivaram, C. S. Hellberg, and B. T. Jonker, Twist angle-dependent atomic reconstruction and moiré patterns in transition metal dichalcogenide heterostructures, ACS Nano 14, 4550 (2020).
  * Weston  _et al._ [2020] A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, J. Zultak, D. Hopkinson, A. Barinov, T. H. Bointon, A. Kretinin, N. R. Wilson, P. H. Beton, V. I. Fal’ko, S. J. Haigh, and R. Gorbachev, Atomic reconstruction in twisted bilayers of transition metal dichalcogenides, Nature Nanotechnology 15, 592 (2020).
  * Sung  _et al._ [2020] J. Sung, Y. Zhou, G. Scuri, V. Zólyomi, T. I. Andersen, H. Yoo, D. S. Wild, A. Y. Joe, R. J. Gelly, H. Heo, S. J. Magorrian, D. Bérubé, A. M. M. Valdivia, T. Taniguchi, K. Watanabe, M. D. Lukin, P. Kim, V. I. Fal’ko, and H. Park, Broken mirror symmetry in excitonic response of reconstructed domains in twisted MoSe2/MoSe2 bilayers, Nature Nanotechnology 15, 750 (2020).
  * Woods  _et al._ [2021] C. Woods, P. Ares, H. Nevison-Andrews, M. Holwill, R. Fabregas, F. Guinea, A. Geim, K. Novoselov, N. Walet, and L. Fumagalli, Charge-polarized interfacial superlattices in marginally twisted hexagonal boron nitride, Nature communications 12, 1 (2021).
  * Yasuda  _et al._ [2021] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Stacking-engineered ferroelectricity in bilayer boron nitride, Science 372, 1458 (2021).
  * Stern  _et al._ [2021] M. V. Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. B. Shalom, Interfacial ferroelectricity by van der Waals sliding, Science 372, 1462 (2021).
  * Weston  _et al._ [2021] A. Weston, E. G. Castanon, V. Enaldiev, F. Ferreira, S. Bhattacharjee, S. Xu, H. Corte-Leon, Z. Wu, N. Clark, A. Summerfield, _et al._ , Interfacial ferroelectricity in marginally twisted 2D semiconductors, arXiv:2108.06489 (2021).
  * Wang  _et al._ [2021] 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, arXiv:2108.07659 (2021).
  * 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).
  * Tong  _et al._ [2020] Q. Tong, M. Chen, F. Xiao, H. Yu, and W. Yao, Interferences of electrostatic moiré potentials and bichromatic superlattices of electrons and excitons in transition metal dichalcogenides, 2D Materials 8, 025007 (2020).
  * Ferreira  _et al._ [2021a] F. Ferreira, V. V. Enaldiev, V. I. Fal’ko, and S. J. Magorrian, Weak ferroelectric charge transfer in layer-asymmetric bilayers of 2D semiconductors, Scientific Reports 11, 13422 (2021a).
  * Carr  _et al._ [2018] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Relaxation and domain formation in incommensurate two-dimensional heterostructures, Phys. Rev. B 98, 224102 (2018).
  * Enaldiev  _et al._ [2020] V. V. Enaldiev, V. Zólyomi, C. Yelgel, S. J. Magorrian, and V. I. Fal’ko, Stacking domains and dislocation networks in marginally twisted bilayers of transition metal dichalcogenides, Phys. Rev. Lett. 124, 206101 (2020).
  * Enaldiev  _et al._ [2021] V. V. Enaldiev, F. Ferreira, S. J. Magorrian, and V. I. Fal’ko, Piezoelectric networks and ferroelectric domains in twistronic superlattices in WS2/MoS2 and WSe2/MoSe2 bilayers, 2D Materials 8, 025030 (2021).
  * Ferreira  _et al._ [2021b] F. Ferreira, S. J. Magorrian, V. V. Enaldiev, D. A. Ruiz-Tijerina, and V. I. Fal’ko, Band energy landscapes in twisted homobilayers of transition metal dichalcogenides, Applied Physics Letters 118, 241602 (2021b).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Sohier  _et al._ [2017] T. Sohier, M. Calandra, and F. Mauri, Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene, Phys. Rev. B 96, 075448 (2017).
  * 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, Physical Review B 54, 11169 (1996).
  * Meyer and Vanderbilt [2001] B. Meyer and D. Vanderbilt, Ab initio study of BaTiO3 and PbTiO3 surfaces in external electric fields, Phys. Rev. B 63, 205426 (2001).
  * Note [1] For monolayers and 2H-bilayers we used the Coulomb truncation in the out-of-plane direction[43].
  * Note [2] Unlike Refs. [55, 56], which use $z$-averaged electric field for computation of dielectric permitivities of bilayers, we express the quadratic amendment to the total energy via out-of-plane displacement field, conserving across every cross-section of the structure. This allows us to avoid uncertainties in $\alpha_{zz}^{\rm 3R}$, which may appear at averaging of the electric field in crystals with a few out-of-plane unit cells [45]. We find that the polarizability, computed for a TMD monolayer, its 2H and 3R bilayer, and thicker 2H films, linearly scales with with the number of layers. This observation contradicts some earlier DFT studies of dielectric susceptibility of TMDs [55, 56] which claimed a pronounce layer-number-dependence, but agrees with the more recent results [49] published by some of the authors of Refs [55].
  * Slizovskiy  _et al._ [2019] S. Slizovskiy, A. Garcia-Ruiz, A. I. Berdyugin, N. Xin, T. Taniguchi, K. Watanabe, A. K. Geim, N. D. Drummond, and V. I. Fal’ko, Out-of-plane dielectric susceptibility of graphene in twistronic and bernal bilayers, arXiv:1912.10067 (2019).
  * Tian  _et al._ [2020] T. Tian, D. Scullion, D. Hughes, L. H. Li, C.-J. Shih, J. Coleman, M. Chhowalla, and E. J. G. Santos, Electronic polarizability as the fundamental variable in the dielectric properties of two-dimensional materials, Nano Letters 20, 841 (2020).
  * Slizovskiy  _et al._ [2021] S. Slizovskiy, A. Garcia-Ruiz, A. I. Berdyugin, N. Xin, T. Taniguchi, K. Watanabe, A. K. Geim, N. D. Drummond, and V. I. Fal’ko, Out-of-plane dielectric susceptibility of graphene in twistronic and bernal bilayers, Nano Letters 21, 6678 (2021).
  * Note [3] The expression for interaction energy of the FE polarization with displacement field naturally comes when assuming a local dielectric permittivity in a continuum medium approximation. Indeed, suppose the FE charges, with plane-averaged density $\delta\rho(z)$, are placed in the medium with local dielectric permittivity $\epsilon_{zz}(z)$. From the Poisson equation and electro-neutrality condition $\intop\ilimits@_{-\infty}^{+\infty}\delta\rho(z)dz=0$ $\partial_{z}(\epsilon_{zz}(z)\partial_{z}\varphi(z))=-\delta\rho(z)/\epsilon_{0}$, we express the potential drop across the layer of charges as $\Delta=\intop\ilimits@_{-\infty}^{+\infty}\partial_{z}\varphi(z)dz=-\intop\ilimits@_{-\infty}^{+\infty}dz\intop\ilimits@_{-\infty}^{z}dz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}=\intop\ilimits@_{-\infty}^{+\infty}dz\intop\ilimits@^{+\infty}_{z}dz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}=\intop\ilimits@\intop\ilimits@_{z<z^{\prime}}dzdz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}$. At the same time, interaction energy of these charges with uniform external out-of-plane displacement field (related to local electric field as $D=\epsilon_{0}\epsilon_{zz}(z)E(z)$) reads as $\delta U=-\intop\ilimits@_{-\infty}^{+\infty}dz\delta\rho(z)\intop\ilimits@_{-\infty}^{z}dz^{\prime}D/\epsilon_{zz}(z^{\prime})\epsilon_{0}=-D\intop\ilimits@_{-\infty}^{+\infty}dz^{\prime}\intop\ilimits@_{-\infty}^{z^{\prime}}dz\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}=-D\intop\ilimits@\intop\ilimits@_{z<z^{\prime}}dzdz^{\prime}\delta\rho(z^{\prime})/\epsilon_{zz}(z)\epsilon_{0}\equiv-D\Delta$. After changing variables, $z\leftrightarrow z^{\prime}$, at the last step of that calculation, we arrive at the relation in Eq. (2).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981).
  * Note [4] At small $D$ (when $y^{\prime 2}\ll 1$) approximate solution of Eq. (5) is given by parabola: $y(x)=2D\Delta(\ell-x)x/(\bar{w}+2\widetilde{w})$.
  * Beal  _et al._ [2018] L. D. Beal, D. C. Hill, R. A. Martin, and J. D. Hedengren, Gekko optimization suite, Processes 6, 106 (2018).
  * Santos and Kaxiras [2013] E. J. Santos and E. Kaxiras, Electrically driven tuning of the dielectric constant in MoS2 layers, ACS Nano 7, 10741 (2013).
  * Laturia  _et al._ [2020] A. Laturia, M. L. Van de Put, and W. G. Vandenberghe, Dielectric properties of hexagonal boron nitride and transition metal dichalcogenides: from monolayer to bulk, npj 2D Materials and Applications 4, 28 (2020).