††thanks: Co-first author††thanks: Co-first author

# Spontaneous interlayer compression in commensurately stacked van der Waals
heterostructures

Nicholas A. Pike<EMAIL_ADDRESS>nanomat/Q-MAT/CESAM and European
Theoretical Spectroscopy Facility, Université de Liège, B-4000 Liège, Belgium
Antoine Dewandre  nanomat/Q-MAT/CESAM and European Theoretical Spectroscopy
Facility, Université de Liège, B-4000 Liège, Belgium  François Chaltin
Department of Chemical Engineering, Université de Liège (B6a), B-4000 Liège,
Belgium.  Laura Garcia  Department of Physics, Université de Liège (B5),
B-4000 Liège, Belgium.  Salvatore Pillitteri  Department of Physics,
Université de Liège (B5), B-4000 Liège, Belgium.  Thomas Ratz  Solid State
Physics, Interfaces and Nanostructures/Q-MAT/CESAM, Université de Liège,
B-4000 Liège, Belgium  Matthieu J. Verstraete  nanomat/Q-MAT/CESAM and
European Theoretical Spectroscopy Facility, Université de Liège, B-4000 Liège,
Belgium  Catalan Institute of Nanoscience and Nanotechnology (ICN2), Campus
UAB, Bellaterra, 08193 Barcelona, Spain.

###### Abstract

Interest in layered two dimensional materials, particularly stacked
heterostructures of transition metal dichalcogenides, has led to the need for
a better understanding of the structural and electronic changes induced by
stacking. Here, we investigate the effects of idealized heterostructuring,
with periodic commensurate stacking, on the structural, electronic, and
vibrational properties, when compared to the counterpart bulk transition metal
dichalcogenide. We find that in heterostructures with dissimilar chalcogen
species there is a strong compression of the inter-layer spacing, compared to
the bulk compounds. This compression of the heterostructure is caused by an
increase in the strength of the induced polarization interaction between the
layers, but not a full charge transfer. We argue that this effect is real, not
due to the imposed commensurability, and should be observable in
heterostructures combining different chalcogens. Interestingly, we find that
incommensurate stacking of Ti-based dichalcogenides leads to the stabilization
of the charge density wave phonon mode, which is unstable in the 1T phase at
low temperature. Mixed Ti- and Zr- heterostructures are still unstable, but
with a charge density wave or a ferroelectric instability.

##  I Introduction

The idea that an ordered stacking of two-dimensional (2D) materials leads to
exotic electronic properties has been around since the discovery of graphene
Novoselov  _et al._ (2004). 2D materials are often used as building-blocks for
"materials-by-design" Novoselov  _et al._ (2016): each layer displays
characteristics that depend on a number of controllable factors including the
chemical composition Pike  _et al._ (2018, 2019), the number of layers, and
the rotation angle between individual layers Wilson  _et al._ (2017); Wang
_et al._ (2019). With the increasing interest in heterostructure devices, it
is important to understand, from a first-principles approach, how their
structural, electronic, and vibrational properties change due to stacking.
This is especially true as experimental limitations are overcome in producing
large-area, transferable, monolayers (ML) of material Jawaid  _et al._ (2017)
and heterostructure devices Iannaccone  _et al._ (2018).

Figure 1: (Color Online) Upper left: Sketches of bulk h-TMDs and t-TMDs with
the transition metal atoms as blue/green spheres, chalcogen atoms in red
spheres, and unit cell coordinate system (middle left). The parallelograms in
black represent the unit cells for each structure. The difference in stacking
sequences for the TMDs can be observed in the lower left of the figure. H-TMDs
have AB stacking in the $z$ direction while the t-TMDs have AA stacking in the
$z$ direction. The geometric thickness ($t$) and vdW gap ($d$) is shown in
relation to the out-of-plane lattice constant ($c$) for all cases. Right side:
Cartoon representation of our calculation process. Given two bulk structures,
A and B in the upper right, we generate the heterostructure (A+B, lower right)
and use the relaxed lattice parameters of A+B to generate the constrained bulk
structures (middle right). Orange vertical bars indicate that the in-plane
lattice parameters of the heterostructure and constrained calculations are
identical.

Experiments combine (preferably air-stable) ML materials (graphene, black
phosphorus, hexagonal boron-nitride, MoS2, etc.) into heterostructure devices,
either by creating lateral junctions via in-plane covalent bonding Wei  _et
al._ (2015) or by growing or stacking them vertically Chen  _et al._ (2019).
These devices are used as biosensors Bolotsky  _et al._ (2019), field-effect
transistors Zhou  _et al._ (2019), photo-detectors Li  _et al._ (2016),
electro-optical devices Wang  _et al._ (2018a); Rivera  _et al._ (2015);
Taghinejad  _et al._ (2019); Lattyak  _et al._ (2019). Cao et. al. Cao  _et
al._ (2018) even showed that at very small angles twisted bilayer graphene can
become superconducting. Stacked heterostructures with different transition
metal dichalcogenides (TMDs) (e.g. MoSe2 on WSe2) have long-lived inter-layer
excitons, resulting in spatially separated electrons and holes Rivera  _et
al._ (2015); Wang  _et al._ (2018b); Torun  _et al._ (2018), an ideal feature
for photovoltaic applications. Stacking of non-TMD layered materials into
incommensurate heterostructures has also been shown to lead to enhancements in
the electronic specific heat of the heterostructure and stabilization of
native charge density waves Bao  _et al._ (2020). Low frequency Raman
measurements of heterostructure materials Lui  _et al._ (2015); Amin  _et al._
(2016); Zhang  _et al._ (2016) can be used to determine the strength of layer
interactions in a heterostructure lattice. The presence and frequency shifts
of in-plane shear modes show commensurate stacking and strain effects. Even
for incommensurate materials, the interlayer breathing mode is always present
and its frequency and width reveal details of the layer distance and bonding
strength.

Recent high-throughput calculations, see for example Refs. 25; 26; 27 and 28,
have discovered numerous two-dimensional materials that could be used to
create devices with unique electric and optical properties. The combinatorial
phase space for 2D heterostructures is immense and there are many unique and
exotic phenomena yet to be explored both experimentally and theoretically. Ab
initio simulations determining the properties of stacked heterostructures have
expanded significantly in the past few years, but are often limited to
commensurate cases, despite the large sizes of the super-lattices that are
observed experimentally.

Recent work by Van Troeye et. al. Troeye  _et al._ (2019) on van der Waals
(vdW) heterostructures used first-principles ingredients to predict the
coherence of a heterostructure from the individual ML elastic constants and
structural parameters. Work by Pizzi et. al. Pizzi  _et al._ (2020) on layer-
dependent interactions in stacked materials composed of the same ML used
symmetry arguments and a spring model to understand layer dependent
vibrational properties in vdW structures. In the hexagonal TMDs, Phillips et.
al. Phillips and Hellberg (2019) and Terrones et. al. Terrones  _et al._
(2013) calculated the electronic structure of experimentally realized bilayer
heterostructures in both AA and AB stacking assuming a commensurate structure
and found both direct and indirect electronic band gaps.

In the following, we investigate heterostructures of vertically stacked TMD
materials using first-principles calculations. This approach ensures a
consistent interpretation of the results across a wide variety of physical
properties. In section II, we describe the structure of our TMD stacks, for
the two crystal symmetries used here. Section III describes our first-
principles methods, and Section IV details and discusses the results of our
calculations for all the commensurate heterostructures, in particular the
interlayer compression. Finally, in Section V, we provide concluding remarks.

##  II Crystal Structure

In their bulk forms, the TMD compounds are stacked layers bound by vdW forces.
They belong to two main symmetry classes: hexagonal (P63/mmc, space group 194)
and trigonal (P$\overline{3}$m1, space group 164). A few metals, such as Re,
produce a closely related triclinic structures (P$\overline{1}$, space group
2). The structures are similar, with layers of trigonal or octrahedral prisms,
but give different stacking orders in the bulk. An example of these two
materials is shown in Fig. 1 in which the bulk structures of generic hexagonal
and trigonal TMDs are shown on the left and the commensurately stacked
heterostructures are on the right. The hexagonal TMDs have AB stacking, as a
result of the 180∘ rotation between the layers, and the trigonal TMDs have AA
stacking. We denote the structures here as, e.g. MoSe2/WSe2, in the notation
of Ref. 33, omitting the rotation angle between the layers. In our case this
angle is always 0, but in reality many intermediate angles can be chosen,
resulting in long range moiré patterns Kang  _et al._ (2013); van Wijk  _et
al._ (2015).

Our heterostructures are composed of two different TMD layers repeating
periodically in the out-of-plane direction. The number of possible
combinations is restricted by studying only combinations of compounds with the
same symmetry. It is important to note that the hexagonal TMD MLs do not
possess inversion symmetry (only the bulk crystal does), whereas trigonal TMD
MLs do, and this carries over to the heterostructures they compose. Thus, for
hexagonal MLs stacked as a heterostructure, the space group of the
heterostructure becomes P$\overline{6}$m2 (space group 187).

Figure 2: Calculated out-of-plane lattice parameter (full height of each bar)
for each of the heterostructures considered here. Each bar is broken down into
contributions from the geometric thickness ($t$) and vdW gap ($d$) of each ML
as described in Figure 1. For each set of three bars, Left: sum of the bulk
values of $t$ and $d$ (Ref. 3). Right: fully relaxed heterostructures. Middle
bar: sum of constrained calculations, for each material with in-plane lattice
matching the relaxed heterostructure. The inset white numbers correspond to
the amount of compression or expansion (in percent) of the individual
thickness/spacing, when compared to bulk values. A and B correspond to the
first and second material given in the title of each triplet of bars. There is
a systematic trend to global compression, which is much stronger (up to 10%
and more) in the mixed chalcogen cases (group II, right 3 sets).

##  III Calculation Methods

We employ density functional theory (DFT) Martin (2004) and density functional
perturbation theory (DFPT) Gonze (1997); Baroni  _et al._ (2001); Gonze and
Lee (1997) calculations, as implemented in the Abinit  software package Gonze
_et al._ (2005, 2009, 2016); Romero  _et al._ (2020), with a plane wave basis
set and norm-conserving pseudopotentials. The pseudopotentials used here are
generated use the fhi98pp code Fuchs and Scheffler (1999), except for W and Ti
which produced inaccurate lattice parameters (errors over 5%) or overly large
band gap energies. For W, we use a pseudopotential generated with the OPIUM
code Walter and Rappe  and for Ti the pseudopotential is generated using the
ONCVPSP software Hamann (2013). In all cases, the pseudopotentials use the
GGA-PBE exchange-correlation functional Perdew  _et al._ (1996) without spin-
orbit coupling. As shown in Ref. 3 the inclusion of spin-orbit interactions is
not necessary in these TMD compounds for the structural and vibrational
properties.

Our calculations of the heterostructures follow careful convergence studies of
the energy cut-off (Ecut) and of the reciprocal space $k$-point mesh for each
of the compounds. A plane wave cutoff energy of E${}_{cut}=$ 50 Ha and a
k-point grid size of $8\times 8\times 8$ ($16\times 16\times 8$ for the
metallic TMD systems) guarantee total energies to within 3 meV per six atom
unit cell for all the compounds. The self-consistent cycle is converged to
total energy differences less than 10-10Ha. For the relaxation of the
structural parameters we use the Broyden-Fletcher-Goldfarb-Shanno minimization
procedure Fletcher (1987) in which we allow both the position of the atoms and
cell shape to change simultaneously, while imposing symmetry. Convergence of
this relaxation procedure stopped when the forces are lower than 1
$\mu$Ha/Bohr. Our DFPT calculations for the inter-atomic force constants used
the same optimized $E_{cut}$ and the relaxed geometries. Static atomic charges
are determined from a Bader atom-in-molecule approach Bader. (1985) within
Abinit .

Bulk | Eg(eV)
---|---
MoS2 | $0.886$
MoSe2 | $0.853$
WS2 | $1.029$
WSe2 | $0.938$
TiS2 | –
TiSe2 | –
ZrS2 | $1.041$
ZrSe2 | $1.004$

Heterostructure | Eavg(eV) | Eg(eV)
---|---|---
MoSe2/WSe2 | $0.895$ | $0.879$
|  |  $1.12$111Ref. 32 Theory
WS2/MoS2 | $0.958$ | $1.049$
|  |  $1.23$222Ref. 31 Theory
TiS2/ZrS2 | – | $0.766$
TiSe2/ZrSe2 | – | $0.159$
WS2/MoSe2 | $0.941$ | $0.722$
|  |  $1.05$11footnotemark: 1
|  |  $1.58$333Ref. 23 Theory
WS2/WSe2 | $0.983$ | $0.817$
|  |  $1.07$11footnotemark: 1
TiS2/TiSe2 | – | $0.059$

Table 1: Calculated electronic band gap energies for each of our
heterostructures and their bulk counterparts. For the heterostructures, we
also report the average of the host material band gaps. – indicates that the
calculated structure is metallic. NB: The calculated values of Refs. 31 and 32
are for commensurate bilayers, not periodic heterostructures.

To take into account the long-range electron-electron correlation, we use the
dispersion scheme given by Grimme Grimme  _et al._ (2010) which is known to
reproduce the interactions of 2D materials quite well VanTroeye  _et al._
(2016); Pike  _et al._ (2018); Sabatini  _et al._ (2016). This dispersion
scheme, known as DFT-D, is based on simple atomic pair-wise terms, with
environment-dependent dispersion coefficients tabulated beforehand using time-
dependent DFT. We showed previously that DFT-D3 in particular is very accurate
compared to experiment in the 2D TMDs Pike  _et al._ (2018). More generally,
the performance of the most-popular vdW functional and dispersion schemes have
been assessed recently Tawfik  _et al._ (2018) and DFT-D3 performs well
compared to more elaborate methods.

For all the vibrational property calculations we numerically converge the
ground state wave function to a relative variance less than $10^{-18}$, and
the first order wave functions to less than $10^{-10}$. The dynamical matrices
are calculated on an irreducible Brillouin Zone wedge of $q$-points
corresponding to an unshifted $4\times 4\times 4$ mesh.

The overall calculation procedure for each heterostructure is as follows:
First, we generate and fully relax the periodically stacked commensurate
heterostructures A/B (third column of data in figure 2). These heterostructure
geometries differ from the bulk due to two main effects: the juxtaposition of
the layers (both short-range “chemical” and long-range dielectric
interaction), and the strain imposed by the artificial commensuration. Second,
to separate these effects, we generate constrained bulk structures for each
component, e.g. A, constraining the in-plane lattice parameters of a bulk unit
cell of A to those of the heterostructure, and allowing the out-of-plane
lattice parameter to relax (middle bars in figure 2). After relaxation we
calculate the electronic and vibrational properties of each system.

##  IV Results

Comparing the heterostructure to isolated MLs, the relaxation of a multi-
layered system must account for the lattice mismatch between the individual
crystalline layers. This mismatch adds a strain energy to the interface, which
generically will compress one layer and expand the other. In reality, 2D vdW
systems often release this strain energy by rotating one layer with respect to
the other, forming commensurate long range relaxation patterns (moiré) or
incommensurate interfaces Troeye  _et al._ (2019). Given that our structures
are forced to be commensurate, the strain energy leads to an intermediate in-
plane lattice parameter, indicated in Table 1 of the supplemental material
(SM) sup , which integrates the elasticity of each layer (the softer material
will accommodate its lattice more). The table also compares the out-of-plane
lattice parameters of our heterostructures to their bulk counterparts Pike
_et al._ (2018) and to measurements of the vdW gap Wilson and Yoffe (1969).

Figure 3: Phonon dispersion relations for the hexagonal heterostructures
projected by layer (sum of projected eigendisplacement weight on the atoms in
the layer). The color indicates the fraction of the phonon mode on each layer.
Blue: 100% first TMD; Red: 100% second TMD; Magenta: hybridized mode. The
optical modes are very well separated by layer, whereas the acoustic manifold
is more mixed, in particular for WS2/MoSe2. The high-symmetry path in momentum
space used for our calculations is shown in Figure 4.

Figure 2 shows a bar graph comparing the out-of-plane lattice parameters in
different conditions: sum of the bulk compounds, of the constrained
heterostructures, and the fully relaxed heterostructure calculations. In each
case, the calculated out-of-plane lattice parameter, $c$ can be subdivided
into contributions from each layer’s geometric thickness and vdW spacing. For
layers A and B, $c=t_{geom,A}+d_{vdW,A}+t_{geom,B}+d_{vdW,B}$ (in the
heterostructures the vdW spacing’s are equivalent and “shared” by both
layers).

Comparing the height of the relaxed heterostructure (right bar in each set of
figure 2) to the sum of the heights of the bulk compounds first relaxed (left
column) then epitaxially constrained (middle) we find there are two groups of
materials. In group I (first 4 sets in figure 2), there is only a small
compression of the out-of-plane lattice parameter (maximum 5%) and most of it
is due to strain and the unequal elastic responses of the two materials: for
idealized equal and isotropic elastic responses, all of the expansion in A
would be a contraction in B and the three bars would be identical.

In group II, there is a strong contraction of the out-of-plane lattice
parameter (around 10% on average), and most of the effect is due to the
interaction between the layers (change from middle to right bars). The main
difference between the two groups is that, in group I, the chalcogens of the
two layers are identical (MoSe2/WSe2, WS2/MoS2, TiS2/ZrS2, and TiSe2/ZrSe2)
whereas, in group II, the chalcogens are different (WS2/MoSe2, WS2/WSe2, and
TiS2/TiSe2). Figure 2 shows that the change in out-of-plane lattice parameter
is mainly located in a contraction of the vdW gap ($d_{vdw}$, blue segments).
Changes in the geometric thickness ($t$ red and orange segments), on the other
hand, are mainly caused by strain. The reason for the vdW gap contraction in
Group II will be explored through several derived properties below.

The compression of the out-of-plane lattice parameter, through the reduction
of the vdW gap, leads to important changes in the electronic and vibrational
properties of the system. Plots of the Kohn-Sham band structures for each of
the compounds are given in SM Figures 1 and 2 sup . To quantify the change, we
compare in Table 1 the calculated DFT electronic band gap energies for the
heterostructures with the average band gap energy of the heterostructure
components (bulk compound references listed as well). We report band gap
energies for the constrained structures in SM Table 2 sup . Interestingly,
within group II, the resulting band gap energies for all semiconducting
heterostructures are smaller than the average and the individual band gaps of
their bulk counterparts, whereas the TiS2/TiSe2 heterostructure becomes
semiconducting. Within group I, the average bulk and heterostructure band gaps
are comparable. In all of the hexagonal cases, the $Q$ valley between $\Gamma$
and $K$ is lower (indirect band gap SM Figure 1), as in the bulk TMD
compounds. This will strongly affect the transport and optical properties of
the stack, but will also be very sensitive to the precise geometry: a bilayer
may be very different from a periodic continuous heterostructure, and strain
can invert the valley ordering.

Figure 4: Phonon dispersion relations of the trigonal heterostructures
projected by layer (sum of projected eigendisplacement weight on the atoms in
the layer). The color indicates the fraction of the phonon mode on each layer.
Blue: 100% first TMD; Red: 100% second TMD; Magenta: hybridized mode. Lower
right corner: High-symmetry path in momentum space used for our band structure
plots from Ref. 56.

Given the changes in the lattice parameters of the heterostructures compared
to their bulk counterparts, we expect changes in the vibrational structure of
these compounds. We use DFPT to calculate the vibrational properties of our
systems, determining the inter-atomic force constants as well as responses to
homogeneous electric fields, giving access to the dielectric constant and the
Born effective charges (BEC). The inter-atomic force constants allow us to
determine the phonon dispersion relations and related properties (e.g. free
energies).

###  IV.1 Vibrational Properties

In Figures 3 and 4, we show the phonon dispersion curves of the
heterostructures. We determine the contribution from each layer by projecting
the phonon mode displacement vectors onto the atoms in each layer, and color
the phonon band structure accordingly. In Figure 3, we show the phonon
dispersion relations for hexagonal heterostructures WS2/MoS2, WS2/MoSe2,
MoSe2/WSe2, and WS2/WSe2. There is a clear separation between the optical
modes from each layer, which is a signature of weak vdW bonding between
layers. The only noticeable mixing between the layers occurs in certain
acoustics modes near $\Gamma$, and avoided crossings which appear with
exchanges of layer character (blue/red switching). In WS2/MoSe2 (group II), a
significantly stronger mixing occurs in the acoustic branches along the in-
plane momentum path. Furthermore, the weak interaction between the layers is
supported by small dispersion in the out-of-plane direction (segments
$\Gamma$-$A$, $K$-$H$ and $L$-$M$).

A number of studies have examined the phonon modes of TMD bilayers such as 22;
23; 24; 57 (our list is not exhaustive and focuses on works with dissimilar
chalcogens). In particular Ref. 22 pioneered the field by tracking the Raman
shifts in different hetero-bilayer combinations, including WSe2/MoS2 and
MoS2/MoSe2. The shear mode disappears in the Raman spectra due to the layer
incommensurability (it survives in our calculations due to the enforced
epitaxy), but the layer breathing mode is always present, with a frequency
between those of the two component bulks, showing strong mode mixing as we see
in this work. Table 2 compares our calculated interlayer breathing mode
frequencies to available literature data. Our frequencies are systematically
higher, due to the imposed commensuration and to the periodic nature of our
structures out of plane (3D periodicity instead of free bilayers), which will
increase the effective interlayer spring constant Liang  _et al._ (2017).

Heterostructure |  $\omega_{LBM}$(cm-1)
---|---
MoS2/MoSe2 |  $\approx$ 32 444Ref. 22 experiment on bilayer
MoSe2/WSe2 | 39.8
WS2/MoS2 | 49.3
| 31.2 555Ref. 24 experiment on bilayer
TiS2/ZrS2 | 33.9
TiSe2/ZrSe2 | 50.1
WS2/MoSe2 | 64.7
| 39 666Ref. 23 theory on bilayer
WSe2/MoS2 |  $\approx$ 32 11footnotemark: 1
WS2/WSe2 | 61.1
TiS2/TiSe2 | 94.0
Table 2: Frequencies of the interlayer breathing mode ($\omega_{LBM}$) for our
commensurate heterostructures divided into group I and group II structures.
Experimental Raman data for two related structures are shown for comparison.

In Figure 4, we present the trigonal heterostructures (TiS2/TiSe2, TiS2/ZrS2,
and TiSe2/ZrSe2) which show a qualitatively different and more complex
behavior. The trigonal heterostructure phonon modes are more strongly
hybridized than the hexagonal ones (including some optical modes as well),
visible in the more numerous mixed magenta bands. There is no splitting
between the acoustic and optical manifolds, which is the case in the bulk as
well: the bonding is softer than in hexagonal TMDs.

Bulk TiS2 and TiSe2 are challenging to simulate because they present a charge-
density wave (CDW), leading to a phonon instability that is very sensitive to
unit-cell volume, pseudopotential, and strain Dolui and Sanvito (2016); Guster
_et al._ (2018); Chen  _et al._ (2015). Experimentally and theoretically, the
phonon instability in TiS2 and TiSe2 appears at the M special point for both
bulk and ML compounds.

For TiS2/TiSe2, we find that stacking dynamically stabilizes the system: the
heterostructure stays metallic, but there are no unstable phonon modes. Recent
work by Liao et. al. Liao  _et al._ (2020) shows that, in these CDW materials,
induced strain is not enough to stabilize the phonon mode, and that charge
transfer is most likely necessary for system stabilization. This is the case
for our Ti based heterostructures as well: both layers become doped as shown
in Table 3, and there is a net charge transfer between the layers.
Stabilization via charge transfer is also thought to occur in non-TMD layered
materials Bao  _et al._ (2020).

We find that TiS2/ZrS2 and TiSe2/ZrSe2 are dynamically unstable: in
TiSe2/ZrSe2 the system remains semi-metallic, and the phonon dynamics resemble
the CDW in TiSe2 Guster  _et al._ (2018), with an instability at M. For
TiS2/ZrS2, the band structure becomes semiconducting, and the layers are
roughly charge neutral. Here the stronger instability has a different origin,
with a Ti-layer optical branch which is imaginary in the full Brillouin zone.
At $\Gamma$ the two lowest modes are polar in-plane displacements of the Ti
against S, i.e. ferroelectric. Given the experimental demonstration of CDW
suppression in heterostructures Bao  _et al._ (2020), it would be very
interesting if this behavior could be confirmed experimentally as well: 2D
ferroelectrics are highly sought after for ultrathin and flexible sensors,
actuators and memory devices. The phonon displacement vectors at M for the
lowest imaginary mode are in-plane, with the two S atoms moving in opposite
directions along a (110) high-symmetry axis. This is similar to the CDW
displacements in bulk, and to the TiSe2/ZrSe2 instability at M. For the second
unstable mode at M, the eigendisplacements follow the (210) in-plane
direction. There is a significant out-of-plane component, and a net movement
of the TiSe2 layer against the ZrSe2 layer, which produces local Peierls-like
dimerization and global buckling of the layers.

###  IV.2 Origin of the vdW contraction

A natural explanation for the contraction of the vdW gap would be charge
transfer between layers, leading to Coulombic attraction which is not present
in the bulk. A strong charge transfer would make the two layers metallic, with
a net electron(hole) excess on layer A(B). This is not the case for the
semiconducting TMDs, as the electron band structures show no metallization
(see SM Figures 1 and 2 for band structure plots sup ).

A weaker effect would be charge polarization in each layer. To quantify charge
transfer and displacement in our heterostructures, we calculate the dynamical
(Born) and static (Bader) charges for each atom, compared to the bulk. The
dynamical charges are shown in Table 3, and systematically increase in the
heterostructure environment compared to the bulk. We have highlighted the
cells in which there are significant differences (> 5%) with respect to the
reported bulk BEC in Ref. 3, which happens almost systematically in the out of
plane direction. We even find that the BEC of WS2 in mixed group II
heterostructures become larger out-of-plane than in-plane (yellow cells). In
hexagonal systems, the counterintuitive sign and abnormal magnitude of the BEC
observed in bulk is still apparent even under compression of the unit cell and
heterostructure stacking Pike  _et al._ (2017).

For the trigonal compounds, one can not compare to the (metallic) bulk, as the
calculated BEC are not meaningful. In semiconducting TiS2/TiSe2, the large
values for Ti atoms are similar to those found, e.g., in oxide perovskites
Ghosez  _et al._ (1998). The large BEC are linked to the small electronic gap
and the sizable dielectric polarizability.

| atom | Z${}^{*}_{xx}$ | Z${}^{*}_{zz}$ | qB,L(e)
---|---|---|---|---
|  | $het$ | $bulk$ | $het$ | $bulk$ |
WS2/MoS2 | Mo | $\pagecolor{grey}-1.18$ | $-1.09$ | $-0.63$ | $-0.63$ | 0.003
S | $\pagecolor{grey}0.59$ | $0.54$ | $\pagecolor{grey}0.33$ | $0.31$
W | $\pagecolor{grey}-0.52$ | $-0.49$ | $\pagecolor{grey}-0.47$ | $-0.43$ | 0.007
S | $0.26$ | $0.25$ | $0.22$ | $0.21$
WS2/WSe2 | W | $-0.48$ | $-0.49$ | $\pagecolor{yellow}-0.92$ | $-0.11$ | -0.027
S | $0.26$ | $0.25$ | $\pagecolor{yellow}0.53$ | $0.21$
W | $\pagecolor{grey}-1.32$ | $-1.24$ | $\pagecolor{grey}-0.81$ | $-0.30$ | 0.035
Se | $0.65$ | $0.62$ | $\pagecolor{grey}0.33$ | $0.39$
WS2/MoSe2 | Mo | $\pagecolor{grey}-2.02$ | $-1.91$ | $\pagecolor{grey}-1.00$ | $-0.95$ | 0.030
Se | $0.96$ | $0.95$ | $\pagecolor{grey}0.44$ | $0.48$
W | $-0.49$ | $-0.49$ | $\pagecolor{yellow}-1.00$ | $-0.11$ | -0.025
S | $0.26$ | $0.25$ | $\pagecolor{yellow}0.54$ | $0.21$
MoSe2/WSe2 | Mo | $-1.99$ | $-1.91$ | $-0.99$ | $-0.95$ | 0.000
Se | $\pagecolor{grey}1.00$ | $0.95$ | $\pagecolor{grey}0.52$ | $0.48$
W | $-1.25$ | $-1.24$ | $\pagecolor{grey}-0.80$ | $-0.30$ | 0.008
Se | $0.62$ | $0.62$ | $0.38$ | $0.39$
TiS2/TiSe2 | Ti | – | – | – | – | -0.032
S | – | – | – | –
Ti | – | – | – | – | 0.030
Se | – | – | – | –
TiS2/ZrS2 | Ti | – | – | – | – | 0.002
S | – | – | – | –
Zr | – | $6.19$ | – | $1.82$ | -0.006
S | – | $-3.17$ | – | $-0.60$
TiSe2/ZrSe2 | Ti | $6.63$ | – | $1.03$ | – | 0.011
Se | $-3.24$ | – | $-0.50$ | –
Zr | $\pagecolor{grey}7.66$ | $6.76$ | $\pagecolor{grey}1.89$ | $1.70$ | -0.007
Se | $\pagecolor{grey}-3.90$ | $-4.12$ | $\pagecolor{grey}-0.94$ | $-0.55$
Table 3: Comparison of dynamical and static charges in the TMD
heterostructures. The calculated BEC for each of the heterostructures are
compared to their bulk counterpart, both in-plane (xx) and out-of-plane (zz).
The first column for each BEC direction corresponds to the heterostructure and
the second column to the bulk. Shaded cells highlight heterostructure values
that differ by more than 5% from the bulk, and yellow cells represent
compounds with a larger out-of-plane than an in-plane value of the BEC. –
indicates that the calculated structure is metallic. Most of the zz values
increase perceptibly, and none of the values decrease: heterostructuring
always increases polarizability with respect to a chemically homogeneous
reference bulk. The last column contains the static Bader charges, summed over
a layer (L) (bulk reference is 0). The heterostructures with identical
chalcogens have negligible charge transfer, within the numerical noise of the
Bader charge integration routine, whereas the hetero-chalcogen cases show
transfers of 0.02 to 0.03 electrons.

Our calculated values of the static Bader charges, for both bulk and
heterostructure systems, are given in Table 3 of the SM sup , and the layer
charges (sum of atomic charges for layers A and B), which are shown in the
last column of Table 3. The charges show small changes with respect to the
bulk compounds, but there is a clear difference between the homo-chalcogen
cases (few milli-electron charges for the layers) and the hetero-chalcogens
(an order of magnitude more), in favor of the sulfide layers.

At this point we have ingredients to understand the shrinking of the
interlayer gaps presenting different chalcogens. Changing the bulk environment
systematically induces a greater polarizability in both materials, generating
additional attractive forces between the layers. The modified charge density
gives rise to an induced polarization within the layers, both static (Bader
charge imbalance) and dynamic (larger BEC, especially out of plane). The
latter is not exclusive to the hetero-chalcogen case, and would mainly be
visible in a modification of the phonon frequencies (LO/TO splitting) compared
to the bulk. The main effect is therefore the static charge redistribution,
polarizing the two layers in opposite ways, and leading to a net compression
of the interlayer spacing. The (dispersive) vdW interaction between differing
chemical species is not expected to yield this kind of anomaly, and we confirm
numerically that the S-Se interaction is simply the average of the S-S and Se-
Se (by construction in the Grimme D3 scheme).

A number of recent works (e.g. References 64 and 65) have examined Janus TMD
structures, where the chalcogen on one side of a layer has been substituted.
These can be made experimentally Zhang  _et al._ (2020) by exposing (e.g.)
MoS2 to a hydrogen plasma and Se vapor. If there are islands or flakes of MoS2
on top of a MoS2 surface, the S against S interface is preserved, and the
Janus flake now presents an intrinsic dipole. Epitaxy is enforced by the
synthesis process. The interlayer distance is also reduced by this
arrangement, as in our calculations. The authors do not specify the mechanism
for the layer attraction, but it seems natural to expect it is electrostatic
as in our work, and not simply vdW interactions. In Ref. 64 bilayers of Janus
TMDs are studied theoretically and again the intrinsic dipole leads to an
electrostatic attraction/repulsion between the layers: head to tail dipoles
yield a smaller interlayer spacing. Again the origin of the contraction
observed here and in Janus bilayers is the same.

##  V Conclusions

In summary, we have calculated the structural, electronic, and vibrational
properties of ideal commensurately-stacked heterostructures containing MoS2,
MoSe2, WS2, WSe2, TiS2, TiSe2, ZrS2, and ZrSe2 using the Abinit  software
package. For dissimilar chalcogen atoms (group II), we find there is a strong
contraction of the c lattice parameter, localized in the van der Waals gap
distances between the layers. We trace this contraction back to an induced
electrostatic effect and a polarization of the layers with a partial spatial
charge transfer, leading to a net attraction between the layers. This effect
is disentangled from the artificially imposed epitaxy, by comparing with an
intermediate bulk structure, which is constrained in-plane to the
heterostructure lattice constant (and relaxed out of plane). For identical
chalcogens (group I) there is only a small compression, which is mainly due to
the epitaxy.

Calculations of the dynamical stability of these compounds, via the inter-
atomic force constants and the phonon spectra, indicate that the
semiconducting heterostructures are stable. In the hexagonal heterostructures
there is a clear separation of the optical vibrational modes for each layer,
indicating a weaker inter-layer interaction in these materials. The acoustic
modes mix more, especially in the hetero-chalcogen group II cases. In the
trigonal heterostructures, we find that the charge density wave instability
from the Ti based compounds is affected strongly: Counter-intuitively, when
one generates a heterostructure of TiS2 and TiSe2, the resulting compound is
dynamically stable due to a small charge transfer between the layers. The
TiS2/ZrS2 combination is semiconducting with a small gap, but dynamically
unstable, with an unstable polar mode at $\Gamma$. TiSe2/ZrSe2 is semi-
metallic, and retains the CDW instability at the M point of the Brillouin
zone. These three evolutions are linked to interlayer charge transfers and not
only strain, and may be observable in experimental bilayers. Ferroelectricity
in TiS2/ZrS2 is a particularly exciting possibility.

While these commensurate heterostructures will not all be the experimental
ground state (preferring to twist and create moiré patterns), they provide
insight into heterostructure property modifications, and guidance to
understand how the individual layers change due to stacking in
heterostructures. We note several effects which should be measurable:
interlayer distances, phonon frequency shifts, stabilization of charge density
waves (TiS2/TiSe2), and the appearance of ferroelectricity in TiS2/ZrS2. The
proximity and layer polarization should also have an impact on transport
through the layers.

###### Acknowledgements.

N.A.P. and M.J.V. gratefully acknowledge funding from the Belgian Fonds
National de la Recherche Scientifique (FNRS) under grants PDR T.1077.15-1/7
and T.0103.19-ALPS. M.J.V. acknowledges a FNRS sabbatical “OUT” grant at ICN2
Barcelona, as well as ULiège and the Fédération Wallonie-Bruxelles (ARC AIMED
G.A. 15/19-09). Computational resources have been provided by the Consortium
des Equipements de Calcul Intensif (CECI), funded by FRS-FNRS G.A. 2.5020.11;
the Zenobe Tier-1 supercomputer funded by the Gouvernement Wallon G.A.
1117545; and by a PRACE-3IP DECI grants 2DSpin and Pylight on Beskow, and
RemEPI on Archer (G.A. 653838 of H2020). This publication is based upon work
of the MELODICA project, funded by the EU FLAG-ERA_JTC2017 call. N.A.P. and
A.D. contributed equally to this work.

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353, 6298 (2016).
  * Pike  _et al._ (2018) N. A. Pike, A. Dewandre, B. Van Troeye, X. Gonze, and M. J. Verstraete, Physical Review Materials 2 (2018), 10.1103/PhysRevMaterials.2.063608.
  * Pike  _et al._ (2019) N. A. Pike, A. Dewandre, B. VanTroeye, X. Gonze, and M. J. Verstraete, Phys. Rev. Mat. 3, 074009 (2019).
  * Wilson  _et al._ (2017) N. R. Wilson, P. V. Nguyen, K. Seyler, P. Rivera, A. J. Marsden, Z. P. L. Laker, G. C. Constantinescu, V. Kandyba, A. Barinov, N. D. M. Hine, X. Xu, and D. H. Cobden, Science 3, 1601832 (2017).
  * Wang  _et al._ (2019) L. Wang, S. Zihlmann, M.-H. Liu, P. Makk, K. Watanabe, T. Taniguchi, A. Baumgartner, and C. Schönenberger, Nano Letters 19, 2371 (2019).
  * Jawaid  _et al._ (2017) A. Jawaid, J. Che, L. F. Drummy, J. Bultman, A. Waite, M.-S. Hsiao, and R. A. Vaia, ACS Nano 11, 635 (2017).
  * Iannaccone  _et al._ (2018) G. Iannaccone, F. Bonaccorso, L. Colombo, and G. Fiori, Nature Nanotechnology 13, 183 (2018).
  * Wei  _et al._ (2015) W. Wei, Y. Dai, Q. Sun, N. Yin, S. Han, B. Huang, and T. Jacob, Phys. Chem. Chem. Phys. 17, 29380 (2015).
  * Chen  _et al._ (2019) F. Chen, Y. Wang, W. Su, S. Ding, and L. Fu, The Journal of Physical Chemistry C 123, 30519 (2019).
  * Bolotsky  _et al._ (2019) A. Bolotsky, D. Butler, C. Dong, K. Gerace, N. R. Glavin, C. Muratore, J. A. Robinson, and A. Ebrahimi, ACS Nano 13, 9781 (2019).
  * Zhou  _et al._ (2019) W. Zhou, J. Chen, P. Bai, S. Guo, S. Zhang, X. Song, L. Tao, and H. Zeng, Research 2019, 21 (2019).
  * Li  _et al._ (2016) M.-Y. Li, C.-H. Chen, Y. Shi, and L.-J. Li, Materials Today 19, 322 (2016).
  * Wang  _et al._ (2018a) S. Wang, H. Tian, C. Ren, J. Yu, and M. Sun, Scientific Reports 8, 12009 (2018a).
  * Rivera  _et al._ (2015) P. Rivera, J. R. Schaibley, A. M. Jones, J. S. Ross, S. Wu, G. Aivazian, P. Klement, K. Seyler, G. Clark, N. J. Ghimire, J. Yan, D. G. Mandrus, W. Yao, and X. Xu, Nature Communications 6, 6242 (2015).
  * Taghinejad  _et al._ (2019) H. Taghinejad, A. A. Eftekhar, and A. Adibi, Opt. Mater. Express 9, 1590 (2019).
  * Lattyak  _et al._ (2019) C. Lattyak, V. Steenhoff, N. Osterthun, M. Vehse, D. Pareek, M. A. Gonzalez, L. Gütay, S. Schäfer, and C. Agert, in  _2019 IEEE 46th Photovoltaic Specialists Conference (PVSC)_ (2019) pp. 1760–1762.
  * Cao  _et al._ (2018) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * Wang  _et al._ (2018b) H. Wang, W. Wei, F. Li, B. Huang, and Y. Dai, Phys. Chem. Chem. Phys. 20, 25000 (2018b).
  * Torun  _et al._ (2018) E. Torun, H. P. C. Miranda, A. Molina-Sánchez, and L. Wirtz, Phys. Rev. B 97, 245427 (2018).
  * Bao  _et al._ (2020) J.-K. Bao, C. D. Malliakas, C. Zhang, S. Cai, H. Chen, A. J. E. Rettie, B. L. Fisher, D. Y. Chung, V. P. Dravid, and M. G. Kanatzidis, “Quasi-two-dimensional heterostructures (k$m_{1-x}$te)(late3) ($m$ = mn, zn) with charge density waves,” (2020), arXiv:2011.11068 [cond-mat.mtrl-sci] .
  * Lui  _et al._ (2015) C. H. Lui, Z. Ye, C. Ji, K.-C. Chiu, C.-T. Chou, T. I. Andersen, C. Means-Shively, H. Anderson, J.-M. Wu, T. Kidd, Y.-H. Lee, and R. He, Phys. Rev. B 91, 165403 (2015).
  * Amin  _et al._ (2016) B. Amin, T. P. Kaloni, G. Schreckenbach, and M. S. Freund, Applied Physics Letters 108, 036105 (2016).
  * Zhang  _et al._ (2016) J. Zhang, J. Wang, P. Chen, Y. Sun, S. Wu, Z. Jia, X. Lu, H. Yu, W. Chen, J. Zhu, G. Xie, R. Yang, D. Shi, X. Xu, J. Xiang, K. Liu, and G. Zhang, Advanced Materials 28, 1950 (2016).
  * Mounet  _et al._ (2018) N. Mounet, M. Gibertini, P. Schwaller, A. Merkys, I. E. Castelli, A. Cepellotti, G. Pizzi, and N. Marzari, Nature Nanotechnology 13, 246 (2018).
  * Choudhary  _et al._ (2017) K. Choudhary, I. Kalish, R. Beams, and F. Tavazza, Sci. Repts. 7, 5179 (2017).
  * Zhang  _et al._ (2019) X. Zhang, A. Chen, and Z. Zhou, WIREs Computational Molecular Science 9, e1385 (2019).
  * Curtarolo  _et al._ (2012) S. Curtarolo, W. Setyawan, G. L. W. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Mehl, H. T. Stokes, D. O. Demochenko, and D. Morgan, Comp. Mater. Sci. 58, 218 (2012).
  * Troeye  _et al._ (2019) B. V. Troeye, A. Lherbier, S. M. M. Dubois, J.-C. Charlier, and X. Gonze, “First-principles prediction of lattice coherency in van der waals heterostructures,” (2019), arXiv:1911.06846 [cond-mat.mtrl-sci] .
  * Pizzi  _et al._ (2020) G. Pizzi, S. Milana, A. C. Ferrari, N. Marzari, and M. Gibertini, “Shear and breathing modes of layered materials,” (2020), arXiv:2011.14681 [cond-mat.mtrl-sci] .
  * Phillips and Hellberg (2019) M. Phillips and C. S. Hellberg, “Commensurate structures in twisted transition metal dichalcogenide heterobilayers,” (2019), arXiv:1909.02495 [cond-mat.mtrl-sci] .
  * Terrones  _et al._ (2013) H. Terrones, F. Lopez-Urias, and M. Terrones, Sci. Repts. 3, 1549 (2013).
  * Tritsaris  _et al._ (2020) G. A. Tritsaris, Y. Xie, A. M. Rush, S. Carr, M. Mattheakis, and E. Kaxiras, Journal of Chemical Information and Modeling 60, 3457 (2020).
  * Kang  _et al._ (2013) J. Kang, J. Li, S.-S. Li, J.-B. Xia, and L.-W. Wang, Nano Letters 13, 5485 (2013).
  * van Wijk  _et al._ (2015) M. M. van Wijk, A. Schuring, M. I. Katsnelson, and A. Fasolino, 2D Materials 2, 034010 (2015).
  * Martin (2004) R. M. Martin,  _Electronic Structure: Basic Theory and Practice Methods_ (Cambridge University Press, 2004).
  * Gonze (1997) X. Gonze, Phys. Rev. B 55, 10337 (1997).
  * Baroni  _et al._ (2001) S. Baroni, S. de Gironcoli, A. D. Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Gonze and Lee (1997) X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
  * Gonze  _et al._ (2005) X. Gonze, G.-M. Rignanese, M. Verstraete, J.-M. Beuken, Y. Pouillon, R. Caracas, F. Jollet, M. Torrent, G. Zerah, M. Mikami, P. Ghosez, M. Veithen, J.-Y. Raty, V. Olevano, F. Bruneval, L. Reining, R. Godby, G. Onida, D. Hamann, and D. Allan, Z. Kristallogr. 220, 558 (2005).
  * Gonze  _et al._ (2009) X. Gonze, B. Amadon, P. M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Cote, T. Deutsch, L. Genovese, P. Ghosez, M. Giantomassi, S. Goedecker, D. Hamann, P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, M. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rignanese, D. Sangalli, R. Shaltaf, M. Torrent, M. Verstraete, G. Zérah, and J. Zwanziger, Comp. Phys. Comm. 180, 2582 (2009).
  * Gonze  _et al._ (2016) X. Gonze, F. Jollet, F. A. Araujo, D. Adams, B. Amadon, T.Applencourt, C. Audouze, J.-M. Beuken, J. Bieder, A. Bokhanchuk, E. Bousquet, F. Bruneval, D. Caliste, M. Cote, F. Dahm, F. D. Pieve, M. Delaveau, M. D. Gennaro, B. Dorado, C.Espejo, G. Geneste, L. Genovese, A. Gerossier, M. Giantomassi, Y. Gillet, D. Hamann, L.He, G. Jomard, J. L. Janssen, S. L. Roux, A. Levitt, A. Lherbier, F. Liu, I. Lukacevic, A. Martin, C. Martins, M. Oliveira, S. Ponce, Y. Pouillon, T. Rangel, G.-M. Rignanese, A. Romero, B. Rousseau, O. Rubel, A. Shukri, M. Stankovski, M. Torrent, M. VanSetten, B. V. Troeye, M. Verstraete, D. Waroquiers, J. Wiktor, B. Xu, A. Zhou, and J. W. Zwanziger, Comp. Phys. Comm. 205, 106 (2016).
  * Romero  _et al._ (2020) A. H. Romero, D. C. Allan, B. Amadon, G. Antonius, T. Applecourt, L. Baguet, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, F. Bruneval, G. Brunin, D. Caliste, M. Cote, J. Denier, C. Dreyer, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, D. R. Hamann, G. Hautier, F. Jollet, G. Jomard, A. Martin, H. P. C. Miranda, F. Naccarato, G. Petretto, N. A. Pike, V. Planes, S. Prokhorenko, T. Rangel, F. Ricci, G.-M. Rignanese, M. Royo, M. Stegel, M. Torrent, M. J. van Setten, B. V. Troeye, M. J. Verstraete, J. Wiktor, J. W. Zwanziger, and X. Gonze, J. Chem. Phys. 152, 124102 (2020).
  * Fuchs and Scheffler (1999) M. Fuchs and M. Scheffler, Comp. Phys. Comm. 119, 67 (1999).
  * (45) E. J. Walter and A. Rappe, http://opium.sourceforge.net/ .
  * Hamann (2013) D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Fletcher (1987) R. Fletcher,  _Practical Methods of Optimization_ , 2nd ed. (John Wiley & Sons, New York, NY, USA, 1987).
  * Bader. (1985) R. F. W. Bader., Acc. Chem. Res. 18, 9 (1985).
  * Grimme  _et al._ (2010) S. Grimme, S. Ehrlich, and L. Goerigk, J. Chem. Phys. 132, 154104 (2010).
  * VanTroeye  _et al._ (2016) B. VanTroeye, M. Torrent, and X. Gonze, Phys. Rev. B 93, 144304 (2016).
  * Sabatini  _et al._ (2016) R. Sabatini, E. Kucukbenli, C. H. Pham, and S. de Gironcoli, Phys. Rev. B 93, 235120 (2016).
  * Tawfik  _et al._ (2018) S. A. Tawfik, T. Gould, C. Stampfl, and M. J. Ford, Phys. Rev. Materials 2, 034005 (2018).
  * (54) _Supplemental material found at [URL] which includes references3; 55; 49; 62; 23, electronic band structure diagrams, lattice parameters, electronic band gap values, and Bader charges._
  * Wilson and Yoffe (1969) J. Wilson and A. Yoffe, Advances in Physics 18, 193 (1969).
  * Setyawan and Curtarolo (2010) W. Setyawan and S. Curtarolo, Compt. Mat. Sci 49, 299 (2010).
  * Liang  _et al._ (2017) L. Liang, J. Zhang, B. G. Sumpter, Q. H. Tan, P. H. Tan, and V. Meunier, ACS Nano 11 (2017), 10.1021/acsnano.7b06551.
  * Dolui and Sanvito (2016) K. Dolui and S. Sanvito, E. Phys. Lett. 115, 47001 (2016).
  * Guster  _et al._ (2018) I. B. Guster, E. Canadell, M. Pruneda, and P. Ordejon, 2D materials 5, 025024 (2018).
  * Chen  _et al._ (2015) P. Chen, Y.-H. Chan, X.-Y. Fang, Y. Zhang, M. Y. Chou, S.-K. Mo, Z. Hussain, A. V. Fedorov, and T. C. Chaing, Nat. Commun. 6, 8943 (2015).
  * Liao  _et al._ (2020) J.-H. Liao, Y.-C. Zhao, Y.-J. Zhao, X.-B. Yang, and Y. Chen, Journal of Applied Physics 127, 044301 (2020).
  * Pike  _et al._ (2017) N. A. Pike, B. VanTroeye, A. Dewandre, G. Petretto, X. Gonze, G. M. Rignanese, and M. J. Verstraete, Phys. Rev. B 95, 201106R (2017).
  * Ghosez  _et al._ (1998) P. Ghosez, J.-P. Michenaud, and X. Gonze, Phys. Rev. B 58, 6224 (1998).
  * Idrees  _et al._ (2019) M. Idrees, H. U. Din, R. Ali, G. Rehman, T. Hussain, C. V. Nguyen, I. Ahmad, and B. Amin, Phys. Chem. Chem. Phys. 21, 18612 (2019).
  * Zhang  _et al._ (2020) K. Zhang, Y. Guo, Q. Ji, A.-Y. Lu, C. Su, H. Wang, A. A. Puretzky, D. B. Geohegan, X. Qian, S. Fang, E. Kaxiras, J. Kong, and S. Huang, Journal of the American Chemical Society 142, 17499 (2020).


11institutetext: Graduate School of Arts and Sciences, University of Tokyo,
Japan
11email<EMAIL_ADDRESS>22institutetext: Department of
Engineering Mathematics, University of Bristol, United Kingdom
22email<EMAIL_ADDRESS>
# The Benefits of Interaction Constraints
in Distributed Autonomous Systems

Michael Crosscombe 11 Jonathan Lawry 22

###### Abstract

The design of distributed autonomous systems often omits consideration of the
underlying network dynamics. Recent works in multi-agent systems and swarm
robotics alike have highlighted the impact that the interactions between
agents have on the collective behaviours exhibited by the system. In this
paper, we seek to highlight the role that the underlying interaction network
plays in determining the performance of the collective behaviour of a system,
comparing its impact with that of the physical network. We contextualise this
by defining a collective learning problem in which agents must reach a
consensus about their environment in the presence of noisy information. We
show that the physical connectivity of the agents plays a less important role
than when an interaction network of limited connectivity is imposed on the
system to constrain agent communication. Constraining agent interactions in
this way drastically improves the performance of the system in a collective
learning context. Additionally, we provide further evidence for the idea that
‘less is more’ when it comes to propagating information in distributed
autonomous systems for the purpose of collective learning.

###### keywords:

collective learning, interaction networks, multi-agent systems, robot swarms

##  1 Introduction

During the development of autonomous systems, it is only natural that we focus
our attention on the aspects of a system which are responsible for producing
the desired behaviours. In distributed autonomous systems the focus is often
on how the interactions between individual agents can lead to effective macro-
level collective behaviours. What may be oversimplified or even forgotten,
however, are the underlying networks which underpin those interactions. A
typical assumption to make is that the agent’s _interaction network_ is merely
a product of the spatial properties of the system, i.e. agents interact
randomly with other agents based only on whether or not they are within the
required distance for interaction to occur; often this takes the form of a
communication range in swarm robotics using platforms such as the Kilobots
with infrared transceivers that are limited to a physical communication radius
of $10$ cm [16, 20]. However, in this paper we argue that we ought to consider
which agents should be allowed to interact, rather than merely considering how
those agents interact.

A common modelling assumption in distributed multi-agent and multi-robot
systems is the “well-stirred system” assumption: In a well-stirred system each
agent is equally likely to encounter any other agent in the system and,
therefore, each interaction is regarded as an independent event [15]. This
equates to a system in which the interaction network of the agents is totally
connected and edges are selected at random to simulate the random interactions
of different pairs of agents. However, this assumption does not consider that
intentional constraints placed on agent interactions may yield better
performance than if the interaction network is assumed to be totally connected
and interaction depends only upon the physical proximity of the agents. More
recently, [18] found that “less is more” when it comes to communication in a
robot swarm and that, by reducing the physical communication radius of the
robots (or the density of the swarm), the system-level performance improves as
a result of there being less communication overall. While this is the most
recent study to highlight the problem of too much communication, [9, Sec. 3.2]
previously demonstrated a similar problem in which agents became “too
cooperative”. This general problem can be traced further back to [14] in which
the authors examined the effect of network structure on system-level
performance. A general survey of recent literature that considers the
connectivity of a multi-agent/robot system conducting collective exploration
is provided in [13].

In previous work we have studied the impact of applying constraints to the
interaction network in a non-physical multi-agent system within a collective
learning setting [5]. We showed that the “less is more” hypothesis holds true
for collective learning when the underlying interaction network’s topology is
a small-world network and that interaction networks with low connectivity lead
to optimal performance for collective learning. However, this study was
limited to studying interaction networks in abstract simulations, without the
ability to consider the physical network layer. Therefore, in this work we
propose to study both the physical and interaction networks in combination, to
determine which layer of connectivity is more important for performance in the
development of distributed autonomous systems.

An outline of the remainder of this paper is as follows: We begin by
describing collective learning as the context of our investigation, before
defining a three-valued model for collective learning. We present results for
a totally-connected interaction network but where the physical network
connectivity is varied, so as to understand the extent to which the physical
network has an impact on system performance. Then, we introduce constraints
over the interaction network by restricting its topology to that of a regular
lattice network, before studying system performance with respect to different
levels of connectivity. Finally, we conclude our findings with some
discussions and highlight our main contributions before detailing directions
for future work.

##  2 Collective Learning in Multi-Agent Systems

Social learning broadly describes the observed phenomena in animals where an
individual’s learning is influenced by observing and interacting with other
individuals or their products [10]. An example of this is the use of
pheromones to facilitate teaching in the tandem-running ant Temnothorax
albipennis [8]. This style of learning forms the basis of many distributed
autonomous systems which, by design, rely on the interactions of individual
agents to produce behaviours more complex than could be demonstrated by any
one individual.

In the last decade, distributed decision-making has been identified as a
foundational problem in swarm systems due to the need for large groups of
agents to cooperate on a shared task by functioning as a collective [2, 17].
In swarm systems, the ability for agents to make decisions in a _distributed_
manner is crucial for their real-world deployment, as robot swarms are
appropriate platforms for deployment in difficult terrain and hazardous
environments, such as for the detection of wildfires [12]. Due to their large
scale, swarm systems often cannot be supervised at the individual level and
must instead rely on human supervision at the swarm level [11], while making
its own decisions and taking actions autonomously, often by reaching a
consensus [19] or through quorum sensing [1]. Much of the literature on
collective decision-making in robot swarms centres on the best-of-$n$ problem;
a bio-inspired problem in which the swarm must decide which is the best option
of $n$ possible alternatives [15]. However, this class of decision-making is
limited to scenarios in which each option is associated with a quality value.
Instead, we focus on a state-of-the-world problem whereby agents must learn
the current state of the world, in order to make informed decisions about
their assigned task.

In this paper we adopt the term collective learning to describe social
learning in a population of artificial agents attempting to reach a consensus
about the state of their environment [4]. In this context, collective learning
is used specifically to describe the combination of two distinct processes:
direct learning based on evidence obtained directly from an agent’s
environment, and indirect learning in the form of an agent fusing its belief
with the belief of another agent. The effects of combining these two processes
has been studied from an epistemological perspective by [7] who argue that
communication between agents acts both to propagate information through the
population, and to correct for errors that occur whilst gathering evidence.
The connection between these two processes has been further investigated more
recently in [3, 6].

##  3 A Three-Valued Model for Collective Learning

We consider a collective learning problem in which a distributed system of
autonomous agents attempts to learn the true state of their environment.
Suppose that the environment can be described by a set of $n$ propositional
variables $\mathcal{P}=\\{p_{1},\dots,p_{n}\\}$ with each relating to a
discrete location, then an agent’s belief about the environment is an
allocation of truth values to each of the propositional variables
$b:\mathcal{P}\rightarrow\\{0,\frac{1}{2},1\\}^{n}$. In our model we adopt a
three-valued propositional logic with the truth values $0$ and $1$ denoting
false and true, respectively, whilst the third truth-value $\frac{1}{2}$ is
used to denote unknown and allows agents to express uncertainty about the
environment. For notational convenience, we can instead represent an agent’s
belief $b$ by the $n$-tuple $B=\langle B(p_{1}),\ldots,B(p_{n})\rangle$. By
doing so, an agent can express uncertainty about a location in the environment
by assigning the truth value $\frac{1}{2}$ to any of the propositional
variables in $\mathcal{P}$. As an example, let us suppose that an environment
can be described by $n=2$ propositional variables, then the belief $B=\langle
0,1\rangle$ would represent that the agent believes $p_{1}$ to be false, while
$p_{2}$ is true. Now let us suppose another agent holds the belief
$B^{\prime}=\langle\frac{1}{2},\frac{1}{2}\rangle$, then the agent is
expressing uncertainty about both of the propositions $p_{1}$ and $p_{2}$,
thereby indicating that the agent is totally uncertain about the environment.

We now define the two processes of direct and indirect learning in the context
of our proposed three-valued model and describe how they combine to produce
collective learning at the macro level in distributed autonomous systems.

#### Evidential updating.

The learning process by which an agent obtains evidence directly from its
environment, e.g., through onboard sensory capabilities, and updates its
current belief to reflect this evidence is referred to as evidential updating.
In our model, an agent first decides which proposition to investigate by
selecting a single proposition at random from those propositions about which
the agent is uncertain, i.e. where $B(p_{i})=\frac{1}{2}$. Then, having
selected a proposition $p_{i}$ to investigate, the agent travels to the
corresponding location in the environment where it gathers evidence. Evidence
takes the form of an assertion about the truth value of the chosen proposition
$p_{i}$ such that
$E=\langle\frac{1}{2},\ldots,S^{*}(p_{i}),\ldots,\frac{1}{2}\rangle$, where
$S^{*}:\mathcal{P}\rightarrow\\{0,1\\}$ denotes the ground truth – often
referred to as the true state of the world – and that the truth value of each
proposition is certainly true or false. Upon obtaining evidence $E$ the agent
then updates its belief to $B|E$ as follows:

$\displaystyle\begin{split}B|E&=\langle B(p_{1})\odot
E(p_{1}),\dots,B(p_{n})\odot E(p_{n})\rangle\\\ &=B\odot E\end{split}$ (1)

Where $\odot$ is a binary operator defined on $\\{0,\frac{1}{2},1\\}$ given in
Table 1. We will elaborate on the reasoning behind our choice of fusion
operator when we introduce the belief fusion process but, for now, notice that
the operator works for evidential updating by preserving certainty over
uncertainty and that updating in this manner does not, therefore, alter the
truth value of any proposition $p_{j}\in\mathcal{P}$ where $p_{j}\neq p_{i}$
because $E(p_{j})=\frac{1}{2}$. An agent repeatedly gathers evidence about
uncertain propositions until it has become totally certain about the
environment, at which point it ceases to look for evidence.

In real-world environments the act of gathering evidence is often subject to
noise, stemming either from the use of noisy sensors – such as thermal sensors
used by UAVs [21] – or from the environment itself. In this model evidence
shall therefore take the following form:

$\displaystyle E(p_{i})=\begin{cases}~{}S^{*}(p_{i})&:~{}\text{ with
probability }1-\epsilon\\\ ~{}1-S^{*}(p_{i})&:~{}\text{ with probability
}\epsilon\end{cases}$ (2)

where $\epsilon\in[0,0.5]$ is a noise parameter denoting the probability that
the evidence is erroneous.

#### Belief fusion.

In collective learning the agents do not exist in isolation. Instead,
individual agents can benefit from being able to interact with, and learn
from, other agents in a shared environment. To this end, we now describe the
indirect learning process by which pairs of agents fuse their beliefs in order
to improve the system’s ability to learn the state of the world.

|  $B^{\prime}(p_{i})$
---|---
$B(p_{i})$ | $\odot$ | $\boldsymbol{0}$ | $\boldsymbol{\frac{1}{2}}$ | $\boldsymbol{1}$
$\boldsymbol{0}$ | 0 | 0 | $\frac{1}{2}$
$\boldsymbol{\frac{1}{2}}$ | 0 | $\frac{1}{2}$ | 1
$\boldsymbol{1}$ | $\frac{1}{2}$ | 1 | 1
Table 1: The fusion operator $\odot$ applied to beliefs $B$ and $B^{\prime}$.

The belief fusion process provides two crucial benefits: Firstly, agents are
able to propagate information they obtain directly from the environment
through the system by fusing their beliefs with those of other agents. We
propose to combine beliefs using the three-valued fusion operator in Table 1
which is applied element-wise to all of the propositional variables
$p_{i}\in\mathcal{P}$. Given two beliefs $B$ and $B^{\prime}$, corresponding
to the beliefs of two agents, the fused belief is then given by:

$\displaystyle B\odot B^{\prime}=\langle B(p_{1})\odot
B^{\prime}(p_{1}),\dots,B(p_{n})\odot B^{\prime}(p_{n})\rangle$ (3)

Secondly, the belief fusion operator helps to correct for errors that occur
when the evidence obtained by the agents is subject to noise. This error-
correcting effect arises when two agents have gathered evidence about the same
proposition $p_{i}$, independently of one another, and encounter a
disagreement. When a pair of agents differ in their belief about the truth
value of $p_{i}$, then the belief fusion process leads them both to become
uncertain with regards to $p_{i}$. As an example, suppose that two agents hold
the beliefs $B_{1}(p_{i})=1$ and $B_{2}(p_{i})=0$, then upon the agents fusing
their beliefs such that $B_{1}\odot B_{2}(p_{i})=\frac{1}{2}$, both agents
will attempt to seek additional evidence about proposition $p_{i}$, either
directly from the environment or indirectly via fusion with other agents,
having become uncertain about the truth value of $p_{i}$. Notice that the
belief fusion process can lead to agents becoming less certain when the fusing
agents disagree about the truth value of a given proposition, while evidential
updating only leads to agents becoming more certain.

#### Measuring performance.

In the context of collective learning agents are attempting to learn a
representation of the true state of the world underlying their environment. As
such, we can measure the accuracy of the average belief of the population with
respect to the true state of the world in the following way: Given a
population of $m$ agents, the normalised difference between each agent’s
belief $B$ and the true state of the world $S^{*}$, averaged across the
population, is the average error of the system;

$\displaystyle\frac{1}{m}\frac{1}{n}\sum_{j=1}^{m}\sum_{i=1}^{n}\left|B_{j}(p_{i})-S^{*}(p_{i})\right|$
(4)

This shall be our primary metric of performance.

##  4 Simulation Environment

(a) At initialisation.

(b) In progress.

(c) At convergence.

Figure 1:  Simulation environment depicting a noisy collective learning
scenario in which $m=20$ agents attempt to learn about $n=126$ locations when
the environment noise $\epsilon=0.2$. Each agent has a communication radius
$C_{r}=20$.

We now introduce the physical simulation environment used to study collective
learning in a multi-agent system. In Figure 1 the simulation environment is
depicted with $20$ agents and $126$ hexagonal locations for the agents to
investigate111We found that $20$ agents were sufficient to complete the task
with performance similar to that of $50$ agents.. Agents are represented by
black Unmanned Aerial Vehicles (UAVs) with either a grey or a blue dashed
circle surrounding them, indicating their physical communication radius
$C_{r}$. If two agents are within each other’s physical radii of communication
then communication between those two agents is _possible_. If the radius is
blue, then an agent is said to be in a communicating state, whilst grey
indicates that it is not communicating. At initialisation (fig. 1(a)), the
environment is greyed out to represent the collective lack of information that
the system holds about its environment. During the collective learning process
(fig. 1(b)), agents begin to build an internal representation of the state of
the world, which can be averaged across the system to give us a general view
of performance. To show this, the colours of the hexagons increase in
intensity representing the collective certainty about each location as
contained in the belief of the average agent. The different colours are used
to represent ocean (blue), beach (yellow), and low and high terrain (light and
dark green), but are purely visual in this context. The brown hexagon is the
initial launch location from which all agents begin exploring. In the case
that the system collectively adopts the incorrect truth value for a given
proposition, the hexagon at its associated location is coloured red, with its
intensity still representing the average certainty with which the agents have
(incorrectly) determined its truth value. At the end of the simulation (fig.
1(c)) the agents have reached a consensus about the state of the world and the
map is in full view, albeit with some errors depicted by red hexagons.

##  5 The Networks Underlying Agent Interactions

(a)

(b)

(c)

(d)

(e)

Figure 2: A snapshot of the simulation during an experiment. In this model two
different types of networks are formed which govern the agent interactions in
the system. (a) A snapshot of the simulation environment during runtime. (b) A
zoomed-in view of 7 agents with a system-wide communication radius $C_{r}=20$.
(c) A representation of the physical network that forms from the agents’
communication radii. (d) A representation of the pre-assigned interaction
network between the visible agents. (e) The resulting communication links that
are available according to edges which are present in both the physical and
interaction networks.

When we consider the design of distributed autonomous systems (e.g. multi-
agent systems or robot swarms), the system-level performance is driven solely
by the interactions of its constituent agents with their environment and with
each other. The consequence of this fact is that we focus our effort on
designing the rules and processes which define how the agents interact, but
often we neglect to consider which agents ought to interact. In the past, we
have often assumed that it is beneficial for agents to be able to interact
freely with any other agent in the system, but recent works have shown this to
be false [5, 18]. To understand the extent to which agents interact in our
simulation, we now distinguish between the physical network and the
interaction network as they occur in our simulation environment.

In Figure 2 we show a snapshot of the simulation during runtime in which we
overlay visual representations of the different networks which are formed
during a single time step $t$. Assuming each agent acts as a node in both a
physical network and an interaction network, then the red edges connect pairs
of agents which are physically capable of communicating, based on the agents
being within their communication range with radius $C_{r}$. Edges in the
physical network are added or removed during each time step $t$ of the
simulation according to which pairs of agents are within communication range
of one another. The yellow edges, meanwhile, depict the connectivity of the
interaction network as determined at the start of each experiment. Unlike in
the physical network, the pairs of agents which are allowed to interact are
immutable, such that no new edges are added, nor existing edges removed after
the network is initialised. The green edges depict which edges exist in both
networks at a given time step $t$. Each edge then corresponds to a pair of
agents which are both able and allowed to interact with one another.

An interaction network can be constructed in whatever manner leads to optimal
system performance, provided that the agents remain capable of interacting
physically in the environment during their task. The default interaction
network in many multi-agent and multi-robot systems is a totally-connected
network. In this paper we seek to constrain interactions by reducing the
connectivity of the network whilst improving system performance. To do this we
shall use regular lattice networks, having discovered in previous work that
such a topology is more optimal even than small-world networks in this setting
[5, 22].

##  6 Multi-Agent Simulation Experiments

For each experimental condition, the system proceeds in the following way: We
initialise $m=20$ agents at the same starting location and with total
uncertainty about the $n=126$ locations, i.e. each agent begins with the
belief $B(p_{i})=\frac{1}{2}$ for $p_{i}\in\mathcal{P}$. Every agent is also
initialised with the same communication radius $C_{r}\in\\{20,\ldots,100\\}$
and a communication frequency $C_{f}\in[0,1]$.

After initialisation, each agent selects a proposition $p_{i}$ at random about
which they are uncertain and begins travelling to the associated location in
search of evidence about $p_{i}$. When they arrive at the location the agent
will observe their environment and obtain evidence $E$ according to eq. 2,
before updating its belief by eq. 1. The agent will then select another
uncertain proposition $p_{i}$ to investigate and begin to travel in the
direction of the new location. At this moment, the agent will also enter a
communicating state with probability $C_{f}$. If the agent does not enter the
communicating state then it continues to its destination to obtain evidence
and the process repeats. If the agent does enter the communicating state,
however, then the agent begins broadcasting its belief to any agents that come
within its communication radii $C_{r}$ as it travels to its destination. The
agent remains in the communicating state until it encounters one or more other
communicating agents within communication range. Upon receiving the belief(s)
of the other agent(s), the agent will cease communicating and update its
belief $B$ by selecting another belief $B\prime$ at random before adopting the
fused belief $B\odot B\prime$ by eq. 3. If the agent remains uncertain about
the location to which it is travelling, then the agent continues to that
location as normal. However if the agent, having fused its belief, becomes
certain about the location during travel then it simply repeats the process of
searching for new evidence by first selecting a new uncertain location from
its belief. Finally, if instead the agent holds a totally certain belief, such
that $B(p_{i})\in\\{0,1\\}$ for all $p_{i}\in\mathcal{P}$, then the agent
repeatedly enters a communicating state in an attempt to reach a consensus
with the rest of the system.

Each experiment runs for a maximum of $30,000$ time steps or until the
population has reached a unanimous consensus about their environment. We
repeat each experiment $50$ times for each parameter combination $C_{r}$ and
$C_{f}$ and for three different noise levels $\epsilon\in\\{0,0.1,0.3\\}$
representing noise-free, low- and moderate-noise scenarios, respectively. The
results are then averages across the $50$ runs. We now present results for
simulation experiments where the underlying interaction network is totally-
connected, in keeping with the “well-stirred system” assumption, in which we
study how changes to the connectivity of the physical network affect the
macro-level dynamics of the system.

###  6.1 Convergence Results for Physical Networks

Figure 3: Average error at steady state against communication radius $C_{r}$
for increasing communication frequency $C_{f}=[0,1]$. Each dotted line shows
the noise $\epsilon\in\\{0,0.1,0.3\\}$. The shaded regions show a $95\%$
confidence interval.

In this section we show convergence results for a static totally-connected
interaction network and focus our attention on the dynamic physical network
for which the connectivity changes based on the chosen communication radius
$C_{r}$. In all figures, the shaded regions show a $95\%$ confidence interval.

In Figure 3 we show the average error of the population (eq. 4) at steady
state against communication radius $C_{r}$. Each solid line shows the results
for different noise $\epsilon$, while the dashed lines of the same colour show
the noise $\epsilon$. Each plot within this figure represents a different
communication frequency $C_{f}$ increasing from $0$ (no communication) to $1$
(unrestricted communication). For $C_{f}=0$ we present a special case in which
agents do not communicate with one another (referred to as “asocial learning”
by [10]). As a result of this total absence of belief fusion, the agents rely
solely upon evidential updating to inform their beliefs and, in the case where
the evidential updating process is noisy, the agents are unable to learn an
accurate representation of their environment. The result is the system
reaching an average error around $\epsilon$ for all communication radii
$C_{r}$. Therefore, $\epsilon$ serves as a simple benchmark of performance and
any reduction in the average error to below $\epsilon$ demonstrated improved
performance.

As we increase $C_{f}$ and agents are able to fuse their beliefs, we see an
immediate and significant reduction in the average error for $C_{f}\geq 0.001$
such that for all noisy scenarios (i.e. for $\epsilon=0.1$ and $0.3$) the
average error is reduced to below the level of noise in each scenario. In the
noise-free scenario (green line) the populations converge to an average error
above $0$ when the communication frequency is too low (i.e. for $C_{f}<0.1$,
suggesting that while the population has converged to a highly accurate
average belief, the population does not reach a consensus across the system.
Here we see that the best performance is measured for a communication
frequency $C_{f}=0.025$, but increasing $C_{f}$ further results in an
increased average error and, therefore, worse performance in the system.

Perhaps surprisingly, increasing the communication radius $C_{r}$ of the
system, and thereby the connectivity of the physical network over which agents
communicate, has little effect on the macro-level performance as measured by
the average error of the system at steady state. The lines in Figure 3 remain
relatively flat across different communication radii $C_{r}$, with the
performance of the system being almost entirely determined by the
communication frequency $C_{f}$. This suggests that the level of connectivity
of the physical network has very little impact on the performance of the
system for the purpose of collective learning when agents can communicate
freely provided they are in range.

##  7 Constraining the Interaction Network

(a)

(b)

(c)

Figure 4: Network diagrams with different levels of connectivity. (a) Regular
lattice network with $k=2$. (b) Regular lattice network with $k=4$.(c)
Totally-connected network.

We now investigate how applying constraints to the interaction network, i.e.
limiting the micro-level dynamics, affects the macro-level performance of our
model for collective learning. In the remainder of this paper, we will study
interaction networks which take the form of regular lattice networks with
different levels of connectivity. A regular lattice network can be defined by
its connectivity $k\in[2,m-1]$ where $k$ determines how many “nearest
neighbours” each agent is connected to when the agents are arranged around a
ring. Examples of lattice networks for $k=2$ and $k=4$ are shown in Figures
4(a) and 4(b), respectively, and are compared with a totally-connected network
equivalent to $k=19$ in Figure 4(c).

###  7.1 Convergence Results for Constrained Interaction Networks

Results are presented for interaction networks with differing levels of
connectivity $k$. We continue to vary the connectivity of the physical network
via the communication radius to observe whether the physical layer, in
conjunction with the interaction layer, has some effect on system performance
that wasn’t observed previously.

(a) Noise $\epsilon=0.1$.

(b) Noise $\epsilon=0.3$.

Figure 5: Average error at steady state against communication radius $C_{r}$
for increasing communication frequency $C_{f}=[0,1]$. Each dotted line shows a
different lattice network with connectivity $k\in\\{2,4,8,16\\}$. The shaded
regions show a $95\%$ confidence interval.

In Figure 5 we show a comparison of the average error of the system at steady
state against communication radius $C_{r}$ for different communication
frequencies $C_{f}$. Each solid line represents a different regular lattice
network of connectivity $k$. In addition, we compare results for noise
$\epsilon=0.1$ (top) against $\epsilon=0.3$ (bottom), and show the noise level
as a red, dotted line. Beginning with Figure 5(a) we can see that the low-
noise scenario converges to an average error at or around $0$ for
communication frequencies $C_{f}<0.2$ with the best performance observed for
$C_{f}=0.1$. We also see that the physical connectivity of the system has
little impact when comparing across communication radius $C_{r}$. In the high-
noise scenario of Figure 5(b), however, we see that the average error remains
above $0$ for all values of $C_{r}$ and $C_{f}$ and that, when $C_{f}\geq
0.1$, there is a sudden increase in average error for $k\geq 4$. Most
importantly, we observe that the networks with the highest connectivity
exhibit the largest increases in average error as $C_{f}$ increases, such that
$k=16$ always results in worse performance than $k=8$, and so forth. An
exception occurs in the high-noise scenario with very low communication
frequencies $C_{f}<0.1$, where a network with $k=2$ performs slightly worse
than more connected networks. However, for $C_{f}\geq 0.1$, $k=2$ performs on
par or better than more connected networks. For example, at the extreme with a
communication frequency $C_{f}=1$, a network of $k=2$ has $\frac{1}{5}$ the
average error compared with $k=16$ for the same communication radius
$C_{r}=20$. When we increase the radius $C_{r}=100$, a network of $k=2$
achieves less than $\frac{1}{2}$ the average error of a network with $k=16$.
Furthermore, the impact of the physical network also becomes more apparent in
the high-noise scenario where, for $C_{f}=1$, we see a clear increase in
average error as the communication radius $C_{r}$ is increased.

To better understand the dynamics of our model and how the underlying
connectivity affects time to convergence, we present trajectories of the
average error against time in Figure 6. Each plot represents a network of
different connectivity $k$, and each line shows that network’s performance for
different noise $\epsilon$. We have fixed $C_{r}=20$ and $C_{f}=0.2$ as the
best performing parameters in this comparison. From left to right, as
connectivity increases, we see that for a high-noise scenario the time to
convergence is much slower, with the $k=2$ network converging after around
$25,000$ time steps. This is reduced to less than $15,000$ time steps for
$k=4$ but with slightly higher average error. Then, as we continue to increase
the connectivity $k$, the network continues to converge in less time but with
greater average error. Therefore, increasing the connectivity worsens
performance when the noise is high.

Figure 6: Average error against time $t$ for $C_{r}=20$ and $C_{f}=0.2$ for
different regular lattice networks with connectivity $k\in\\{2,4,8,16\\}$ and
a totally-connected network (i.e. $k=19$). Each dotted line shows the noise
$\epsilon\in\\{0,0.1,0.3\\}$. The shaded regions show a $95\%$ confidence
interval.

These results suggest that, for higher communication frequencies, the
combination of both the interaction network and the physical network is an
important consideration. It is also clear that as the interaction network
becomes increasingly connected, the connectivity of the physical network has
less of an impact, and that the inverse is true where, for a less connected
interaction network, the physical network has greater impact overall.
Therefore, both networks are important for considering optimum performance in
distributed autonomous systems.

##  8 Conclusion and Future Work

In this paper we have proposed to study the effects of both the physical and
interaction networks on the collective learning performance of a multi-agent
system. We have demonstrated that, when the underlying interaction network is
totally connected, the physical network has minimal impact on performance in
our model. However, when the connectivity of the interaction network is
greatly reduced the system performs better, achieving lower average error when
compared with highly-connected interaction networks under the same
experimental conditions. Limiting the interactions between agents also results
in improved robustness to noise, where increased connectivity at the physical
layer only worsens performance for the same limited interaction network.
Overall, performance is best when the connectivity of both the physical and
interaction networks is severely reduced, e.g., for a connectivity of $k=2$
and a communication radius $C_{r}=20$.

To our knowledge, this is the first such study explicitly investigating the
impact of both the physical and interaction networks in combination in a swarm
setting. While it is difficult to draw strong, broad conclusions from this
study, we hope that this will encourage other researchers working on
collective behaviours to consider the importance of the underlying interaction
network and its level of connectivity. Obviously, this study is limited to
simulation studies of the effects of the interplay between the physical and
interaction network, but a more theoretical study of this interplay is
important to understand the reason behind the observed dynamics.

Up to now we have only considered regular lattice networks and small-world
networks as the basis for interaction networks. These networks possess certain
desirable properties which may be more conducive to the processes involved in
collective learning than may other, different, types of networks. Furthermore,
we have only considered the effect of network connectivity in the context of
collective learning. This effect may not necessarily be observed for other
kinds of collective behaviours. As such, we intend to further explore how
network connectivity, and information sharing more generally, impacts
different kinds of collective behaviours in multi-agent systems. We are also
interested in understanding how a more dynamic environment might change the
optimal level of connectivity in the networks, as less connectivity may lead
to outdated information being retained by the population which is likely to
affect performance over time.

##  9 Acknowledgements

This work was funded and delivered in partnership between Thales Group,
University of Bristol and with the support of the UK Engineering and Physical
Sciences Research Council, ref. EP/R004757/1 entitled “Thales-Bristol
Partnership in Hybrid Autonomous Systems Engineering (T-B PHASE).”

## References

  * [1] Bechon, P., Slotine, J.J.: Synchronization and quorum sensing in a swarm of humanoid robots.  arXiv preprint arXiv:1205.2952 (2012)
  * [2] Brambilla, M., Ferrante, E., Birattari, M., Dorigo, M.: Swarm robotics: a review from the swarm engineering perspective.  Swarm Intelligence 7(1), 1–41 (2013)
  * [3] Crosscombe, M., Lawry, J.: A model of multi-agent consensus for vague and uncertain beliefs.  Adaptive Behavior 24(4), 249–260 (2016)
  * [4] Crosscombe, M., Lawry, J.: Collective preference learning in the best-of-n problem.  Swarm Intelligence 15(1), 145–170 (2021)
  * [5] Crosscombe, M., Lawry, J.: The impact of network connectivity on collective learning.  In: International Symposium Distributed Autonomous Robotic Systems, pp. 82–94. Springer (2021)
  * [6] Douven, I.: Optimizing group learning: An evolutionary computing approach.  Artificial Intelligence 275, 235–251 (2019)
  * [7] Douven, I., Kelp, C.: Truth approximation, social epistemology, and opinion dynamics.  Erkenntnis 75(2), 271–283 (2011)
  * [8] Franks, N.R., Richardson, T.: Teaching in tandem-running ants.  Nature 439(7073), 153–153 (2006)
  * [9] Hamann, H.: Superlinear scalability in parallel computing and multi-robot systems: Shared resources, collaboration, and network topology.  In: Architecture of Computing Systems–ARCS 2018: 31st International Conference, Braunschweig, Germany, April 9–12, 2018, Proceedings 31, pp. 31–42. Springer (2018)
  * [10] Heyes, C.M.: Social learning in animals: categories and mechanisms.  Biological Reviews 69(2), 207–231 (1994)
  * [11] Hogg, E., Harvey, D., Hauert, S., Richards, A.: Evolving robust supervisors for robot swarms in uncertain complex environments.  In: F. Matsuno, S.i. Azuma, M. Yamamoto (eds.) Distributed Autonomous Robotic Systems, pp. 120–133. Springer International Publishing, Cham (2022)
  * [12] Innocente, M.S., Grasso, P.: Self-organising swarms of firefighting drones: Harnessing the power of collective intelligence in decentralised multi-robot systems.  Journal of Computational Science 34, 80–101 (2019)
  * [13] Kwa, H.L., Leong Kit, J., Bouffanais, R.: Balancing collective exploration and exploitation in multi-agent and multi-robot systems: A review.  Frontiers in Robotics and AI 8 (2022)
  * [14] Lazer, D., Friedman, A.: The network structure of exploration and exploitation.  Administrative science quarterly 52(4), 667–694 (2007)
  * [15] Parker, C.A., Zhang, H.: Cooperative decision-making in decentralized multiple-robot systems: The best-of-n problem.  IEEE/ASME Transactions on Mechatronics 14(2), 240–251 (2009)
  * [16] Rubenstein, M., Ahler, C., Hoff, N., Cabrera, A., Nagpal, R.: Kilobot: A low cost robot with scalable operations designed for collective behaviors.  Robotics and Autonomous Systems 62(7), 966–975 (2014).  Reconfigurable Modular Robotics
  * [17] Schranz, M., Umlauft, M., Sende, M., Elmenreich, W.: Swarm robotic behaviors and current applications.  Frontiers in Robotics and AI 7 (2020)
  * [18] Talamali, M.S., Saha, A., Marshall, J.A., Reina, A.: When less is more: robot swarms adapt better to changes with constrained communication.  Science Robotics 6(56), eabf1416 (2021)
  * [19] Valentini, G.: Achieving consensus in robot swarms.  Studies in computational intelligence 706 (2017)
  * [20] Valentini, G., Ferrante, E., Dorigo, M.: The best-of-n problem in robot swarms: Formalization, state of the art, and novel perspectives.  Frontiers in Robotics and AI 4 (2017)
  * [21] de Vries, J.: Image processing and noise reduction techniques for thermographic images from large-scale industrial fires.  In: Proc. of the 12th International Conference on Quantitative InfraRed Thermography (QIRT) (2014)
  * [22] Watts, D., Strogatz, S.: Collective dynamics of ‘small-world’ networks.  Nature 393, 440–442 (1998)