#  Symmetry-based computational search for novel binary and ternary 2D
materials

Hai-Chen Wang  Institut für Physik, Martin-Luther-Universität Halle-
Wittenberg, D-06099 Halle, Germany  Jonathan Schmidt  Institut für Physik,
Martin-Luther-Universität Halle-Wittenberg, D-06099 Halle, Germany  Miguel A.
L. Marques<EMAIL_ADDRESS>Institut für Physik, Martin-
Luther-Universität Halle-Wittenberg, D-06099 Halle, Germany  Ludger Wirtz
Department of Physics and Materials Science, University of Luxembourg, 162a
avenue de la Faïencerie, L-1511 Luxembourg, Luxembourg  Aldo H. Romero
<EMAIL_ADDRESS>Department of Physics, West Virginia University,
Morgantown, WV 26506, USA  Department of Physics and Materials Science,
University of Luxembourg, 162a avenue de la Faïencerie, L-1511 Luxembourg,
Luxembourg

###### Abstract

We present a symmetry-based exhaustive approach to explore the structural and
compositional richness of two-dimensional materials. We use a “combinatorial
engine” that constructs potential compounds by occupying all possible Wyckoff
positions for a certain space group with combinations of chemical elements.
These combinations are restricted by imposing charge neutrality and the
Pauling test for electronegativities. The structures are then pre-optimized
with a specially crafted universal neural-network force-field, before a final
step of geometry optimization using density-functional theory is performed. In
this way we unveil an unprecedented variety of two-dimensional materials,
covering the whole periodic table in more than 30 different stoichiometries of
form AnBm or AnBmCk. Among the found structures we find examples that can be
built by decorating nearly all Platonic and Archimedean tesselations as well
as their dual Laves or Catalan tilings. We also obtain a rich, and unexpected,
polymorphism for some specific compounds. We further accelerate the
exploration of the chemical space of two-dimensional materials by employing
machine-learning-accelerated prototype search, based on the structural types
discovered in the exhaustive search. In total, we obtain around 6500
compounds, not present in previous available databases of 2D materials, with
an energy of less than 250 meV/atom above the convex hull of thermodynamic
stability.

##  I Introduction

Since the synthesis of single graphene layers Novoselov  _et al._ (2004), two-
dimensional (2D) materials have attracted a significant interest from the
community. Their relevance extends through different research fields, such as
catalysis, electronic transport, optical properties, and topological
properties. However, even though great effort has been spent on investigating
the vast chemical space for bulk, three-dimensional (3D) compounds, the
chemical space for 2D materials is still relatively unexplored. Experimental
synthesis efforts have focused only on a few dozens of structures, mostly
obtained by exfoliation of known 3D layered materials.Novoselov  _et al._
(2005)

On the computational side, we can find a few online databases of 2D materials,
such as Materials Cloud (MC2D) Mounet  _et al._ (2018), V2DB Sorkun  _et al._
(2020), 2DMatpedia Zhou  _et al._ (2019), and the Computational 2D Materials
Database (C2DB) Haastrup  _et al._ (2018); Gjerding  _et al._ (2021); Lyngby
and Thygesen (2022). These databases were built starting from 3D databases, by
theoretically exfoliating single-layers from layered, van der Waals bound
compounds. At the moment the vast majority of known 2D materials correspond to
binaries Sorkun  _et al._ (2020); Song  _et al._ (2021). An exception is the
very recent addition of materials discovered via a crystal diffusion
variational autoencoder in Ref. Lyngby and Thygesen (2022). Due to the recency
of this addition, we still perform comparisons with respect to the previous
version of the C2DBGjerding  _et al._ (2021).

These 2D databases are newer, and considerably smaller, than their three-
dimensional counterparts, e.g. Materials Project Jain  _et al._ (2013), the
crystallographic open database (COD) Gražulis  _et al._ (2009), the Cambridge
structural database Allen (2002), the NIST crystallographic database Van Hove
_et al._ (2002), OQMD Saal  _et al._ (2013), AFLOW Curtarolo  _et al._ (2012),
materials cloud Talirz  _et al._ (2020), and many others.

All these initiatives were seeded by experimental crystal structures stored in
the Inorganic Crystal Structure Database (ICSD) and other experimental
databases. In fact, the creation of the ICSD in 1912 Bergerhoff and Brown
(1987); Zagorac  _et al._ (2019); Belsky  _et al._ (2002), paved the way to
the systematic study of the relationship between crystal structure and
materials properties. To complete experimental data, many databases (both 2D
and 3D) also contain results from high-throughput studies (often accelerated
by machine learning).

High-throughput searches are responsible for a majority of the calculations in
the large theoretical databases like AFLOW Curtarolo  _et al._ (2012), OQMD
Saal  _et al._ (2013) and DCGATSchmidt  _et al._ (2022a). Traditional high-
throughput searches solely rely on simple empirical rules to select candidate
materials for evaluation with density functional theory (DFT) and consequently
result in a large number of highly unstable systems. A particularly popular
approach is prototype search, where new materials are hypothesized by
transmuting the chemical elements in a known crystal structure (often stemming
from ICSD). In some cases, all combinations of chemical elements are taken
into account, while in other cases arguments based on charge neutrality,
atomic or ionic radii, etc. are used to circumvent the combinatorial nature of
the problem. In comparison to these rule-based selections, machine learning
algorithms have allowed us to speed up the scanning of known prototypes by a
factor of up to $\sim$ 2000Park and Wolverton (2020); Schmidt  _et al._
(2022b) by predicting the stability of all respective compounds and filtering
out the most unstable predictions.

No matter the selection rule, the next step consists in structural
optimization of each proposed compound, invariantly using DFT as the workhorse
method Hohenberg and Kohn (1964); Kohn and Sham (1965). Analysis of
thermodynamic stability can then be achieved by computing the formation energy
or the distance to the convex hull. With this or similar approaches, databases
have grown considerably and can now sometimes reach millions of crystal
structures.

While the success of using chemical combinatorics is recognized, it has been a
substantial handicap for predicting new 2D materials. The number of known
materials prototypes is small, which reduces the number of possibilities.
Different research groups have considered different strategies to address this
issue, often resorting to machine learning methods. This paper presents an
entirely different approach, which is not based on motifs or chemical
substitutions. Instead, we create all possible species combinations for binary
(and ternary) systems for specific two-dimensional space groups. Therefore,
based on symmetries and chemical criteria, we can arrive at a sizeable two-
dimensional crystal structure database and a diverse set of structural
prototypes. The crystal shapes show a variety of bondings and forms absent in
existing databases. As our search is systematic and exhaustive, crystal
structures contain a large number of different chemical formulae as well as
almost all possible Wyckoff positions allowed by the space groups. Here we
focus on two-dimensional materials but our approach is general such that it
can also be applied to three-, one- or even zero-dimensional structures.

This article is structured as follows. We start by discussing in detail the
systematic approach to discover 2D crystal structures and our strategies to
accelerate the search. We then present an overview of the materials we
discover, giving a few examples of structural diversity and polymorphism. In
the following, we discuss machine-learning accelerated prototype search based
on the wealth of prototypes obtained. In the Appendix following the
conclusions, we give details on our methodology.

##  II Strategy

Figure 1: Summarized flowchart for choosing potential materials using group
symmetries and their Wyckoff positions.

Figure 1 and Fig. 2 summarize our approach for the generation of materials.
The first step corresponds to a combinatorial workflow that creates
hypothetical compounds. The initial input parameter is the number of desired
species in the particular material. We consider most elements within the
periodic table, starting from Li until Bi, including first-row rare earth
elements. We exclude, however, radioactive elements and rare gases, namely At,
Tc, Pr, Pm, He, Ne, Ar, Kr, Xe, Rn. We then generate all possible combinations
of the different elements selected from the periodic table. For example, a
total of 2701 combinations are obtained with no repeated elements for a binary
compound.

The second parameter is the two-dimensional space group. We use the table of
layer group symmetries, created by considering the wallpaper group and adding
reflections in the perpendicular direction. A full description of the possible
groups is given in Refs. Litvin and Kopsky, 2002; Aroyo  _et al._ , 2011,
2006. To generate the atomic positions provided by the layer space group for
all possible Wyckoff positions (WPs), we used the package PyxTal Fredericks
_et al._ (2021). A list of the possible WPs can be found in the Bilbao
Crystallographic Server. For this publication, we selected the space groups
summarized in Table 1 (other space groups are in progress).

The next step is the creation of all possible combinations of the WPs for each
chemical element in the combined list. For example, if the number of WPs is
four, we get four factorial posibilities. This strategy allows different WPs
for the same species and therefore broadens the number of possible
stoichiometries (i.e., a chemical species can occupy more than one WP). We
then create a product of this list associated with the number of selected
species. This selection leads to the definition of a chemical composition
based on the occupation of the different WPs for each species in the compound.
In parallel, we create a list of possible oxidation states of the considered
species. We make all possible combinations without replacement for each
element from this list, and we create a product of the different list
elements. We used the experimentally most common oxidation states, as they
will have considerably larger potential to be synthesized (the selected
oxidation states are included in the Supplementary Information). Finally, a
compound is created from the provided number of species, the combination of
WPs, and the oxidation state.

Table 1: Crystallographic summary of the layer groups considered in this work: the space group symbol (and number in parenthesis), the Wyckoff positions (and site symmetries in parenthesis). We also show the number of binary ($N^{(2)}_{\text{tot}}$) and ternary systems ($N^{(3)}_{\text{tot}}$) generated by our combinatorial engine and the number of entries that were found below 250 meV/atom from the convex hull of stability ($N^{(2,3)}_{<0.25}$). Space Group | Wyckoff Positions | $N^{(2)}_{\text{tot}}$ | $N^{(2)}_{<0.25}$ | $N^{(3)}_{\text{tot}}$ | $N^{(3)}_{<0.25}$
---|---|---|---|---|---
$p1$ (01) | 1a (1) | 225 | 5 |  |
$p$11m (04) | 2b (1), 1a (..m) | 1321 | 339 |  |
$p$11a (05) | 2a (1) | 225 | 47 | 1944 | 307
$p211$ (08) | 2c(1), 1b(2..), 1a(2..) | 6645 | 623 |  |
$p2_{1}11$ (09) | 2a (1) | 225 | 40 | 1944 | 273
$c211$ (10) | 4b (1), 2a (2..) | 1321 | 153 |  |
$p$b11 (12) | 2a (1) | 225 | 15 | 1944 | 448
$c$m11 (13) | 4b (1), 2a (m..) | 1321 | 268 |  |
$p2_{1}$/b11 (17) | 4c (1), 2b (-1), 2a (-1) | 3129 | 220 |  |
$p2_{1}2_{1}2$ (21) | 4c(1), 2b(..2),2a(..2) | 6645 | 398 |  |
$p$b$2_{1}$m (29) | 4b (1), 2a (..m) | 1321 | 140 |  |
$p$b2b (30) | 4c (1), 2b (.2.), 2a (.2.) | 6645 | 228 |  |
$p$m2a (31) | 4c (1), 2b (m..), 2a (.2.) | 6645 | 478 |  |
$p$m$2_{1}$n (32) | 4b (1), 2a (m..) | 1321 | 148 |  |
$p$b$2_{1}$a (33) | 4a (1) | 225 | 19 | 1944 | 139
$p$b2n (34) | 4b (1), 2a (.2.) | 1321 | 74 |  |
$c$m2e (36) | 8c (1), 4b (m..), 4a (.2.) | 6645 | 383 |  |
p31m (70) | 6d (1), 3c (..m), | 15728 | 783 |  |
2b (3..) 1a (3.m) |  |

After the material is obtained from the previous step, and before we perform a
complete electronic structure calculation, we conduct a materials screening,
which allows us to reduce the number of compounds to be fully considered. For
the screening, we used rules implemented in SMACT Davies  _et al._ (2016). In
this package, decisions are made based on stoichiometry. The first rule is to
have only neutral compounds, which can be easily computed from the
stoichiometry and the oxidation states. Next, we use the so-called Pauling
test for neutral materials, which requires that positive ions are lower in
electronegativity than negative ions.

This process is summarized in Fig. 1. The approach followed here does not
consider intermetallics of unusual oxidation states as a possible outcome in
the search, but besides that, it is fully exhaustive. We show in the last
section how we can go beyond this limitation efficiently using machine
learning.

After screening, we use the main properties of a given material, such as
oxidation states, stoichiometry, and Wyckoff positions to generate the
potential compounds. With this information, we use the PyxTal utility
Fredericks  _et al._ (2021) to create a unit cell, the chosen dimensionality,
and the number of species. When WPs have internal degrees of freedom, PyxTal
tries to create a unit cell with the provided symmetry constraints. First, the
cell directions are selected according to the space group. Then, the WPs are
generated from the symmetry operations, and, if there are internal degrees of
freedom, they are set randomly. Next, the cell parameters and the volume are
determined, assuming that each atom has a radius equal to the covalent bond
radius. Finally, a density is obtained with the cell volume and atomic masses,
which is compared with a minimum density. If the cell density is smaller than
a minimum density (set to 0.75 in scaled units), the package attempts first to
define the atomic positions randomly (setting the number of attempts to 40),
and, in case this fails, it tries to change the cell parameters (10 times) and
repeat the generation of the atomic position. If a cell can not be defined,
the structure generation is considered unsuccessful, and a new material is
considered. A summary of the pipeline followed in this step is represented in
figure 2.

Figure 2: Pipeline for the materials screening, from materials chosen
following flowchart 1.

Based on the previous discussed steps, we generate systematically two
dimensional structures for the space groups shown in Table 1. In this table we
also include the corresponding Wyckoff positions, the site symmetry and the
number of different compounds generated for each space group.

The next step is the geometry optimization. Unfortunately, the initial
structures are usually very far away from equilibrium, making structural
optimization with DFT cumbersome. To increase the efficiency of our workflow
we decided, therefore, to perform an intermediate geometry optimization step
using a universal neural-network force-field Chen and Ong (2022). At this
point we remove from our dataset materials that are too thick (using a
threshold of 7.5 Å), predicted to be too unstable by the machine learning
model (more than 600 meV/atom from the hull, corresponding approximately to
twice the MAE of the original model), or that were already included in C2DB.

The use of machine learning force fields solves several technical problems:
The pre-converged geometries are, in most cases, already quite good, only
requiring a few steps of geometry optimization using DFT. They also allow us
to discard many repeated and very high-energy structures. After the DFT
geometry optimization we evaluate the distance to the convex hull of
stability. We use the convex hull of Ref. Schmidt  _et al._ , 2022a, b that is
considerably larger than the one of the Materials Project Jain  _et al._
(2013), in particular in what concerns the ternary (and quaternary) sector.
Consequently, our distances to the hull are sometimes larger than in other 2D
databases.

##  III 2D Materials

(a)

(b)

Figure 3: Periodic tables showing the frequency of the chemical elements in
binary compounds below 250 meV/atom in (a) our work and (b) C2DB. We emphasize
that our entries do not include, by construction, any of the compounds of
C2DB. Figure 4: Pie chart of the frequency of different general chemical
formula for compounds with distances to the convex hull smaller than 250
meV/atom in our database. The category other includes 29 general formulas.

It turns out that our workflow was able to arrive at the large majority of
systems already present in C2DB. This is particularly true for binary systems,
as these were more extensively investigated than ternaries (see Table 1).
This, in our opinion, fully validates our workflow.

Figure 3 presents a comparison of the binary materials present in our database
(excluding the ones found in C2DB) compared to C2DB. For the discussion, we
only took into account the materials that are within a distance of 250
meV/atom from the convex hull of stability, that corresponds loosely to the
definition of “high-stability” in C2DB Gjerding  _et al._ (2021). Note that
for consistency we have reoptimized the C2DB structures using our convergence
criteria and our selected set of pseudopotentials.

Figure 5: Diversity of structures obtained by our procedure. We label in bold
the convex tiling from which the structure can be derived. Next to the
chemical formula, in parenthesis, is the distance to the convex hull (in
meV/atom) of the specific compound. Note that the structure depicted is not
necessarily the lowest 2D phase for the composition we found in our search.

We find 2D compounds across the whole periodic table, including some with
lanthanides that have been up to now excluded from previous works. Not
surprisingly, the majority of compounds includes a non-metal element (due to
the requirement of charge neutrality), leading to the pronounced peaks for O,
S, Se, Te, F, Cl, Br, I, etc. The figure also reveals some differences in the
prevalence of certain elements between our dataset and C2DB. For example, we
find considerably more compounds with F than with O, while in C2DB we observe
the opposite behavior. As our approach is to a large extent systematic in what
concerns chemical compositions and geometries, we believe that the differences
are explained by a bias already present in ICSD and other databases that were
used to seed the 2D databases. For example, it is well know that oxides are
overrepresented in experimental works as they can be more easily synthesized
and are often stable in air.

Other conclusions can be drawn from Fig. 3. For example, it can be seen
clearly that first-row non-metals have much more difficulty in forming low-
energy compounds than other non-metals in the same group. This is a
consequence of the Singularity Principle Restrepo  _et al._ (2004), i.e., that
the chemistry of the these elements is often different to the later members of
their respective groups. Furthermore, elements like N, O, C, and F form very
strong directional covalent bonds that leave little room for small adjustments
of the structure. In what concerns the metallic elements, the fourth-period
transition elements from Ti to Cu (and in particularly this last one),
together with late group III-A (In and Tl) seem to form easily 2D compounds.

The diversity of stoichiometries is illustrated in Fig. 4. As our largest
emphasis was in binary compounds it is not surprising to that the most
represented stoichiometries are binary. Among these, the simple AB2, AB3,
A2B3, etc. dominate the low-energy structures. This fact can be easily
understood by the requirement of charge neutrality, and the fact that most
non-metals have oxidation states of -I, -II, or -III. As such, the same
situation can be found for bulk, 3D semiconductors and insulators. However, we
do find a long list of other stoichiometries (more than 30), and these often
reveal very interesting and unexpected structures.

In Fig. 5 we give a glimpse of the diversity of structural motifs found by our
method. Note that this is far from a complete list of all 2D structures found,
and that we concentrate on unusual arrangements that go way beyond the most
common square and hexagonal lattices. We emphasize that these motifs appeared
naturally in our workflow and were by no means constructed by hand.
Interestingly, we easily found examples that can be derived from the majority
of the different Euclidean uniform tilings, both Platonic and Archimedean as
well as their dual Laves or Catalan tilings. Moreover, many of these tilings
seem to be unique to the two-dimensional world, as no layered 3D material is
known to possess them.

The first two structures can be derived from a truncated square and a rhombic
tiling. In the first case, Cs2Br2 squares are connected forming regular empty
octahedra, leading to a rather open lattice. In the second, Se3O3 rectangular
units form bonds along the corners, leading to flattened octahedra. We then
present an example of a Pythagorean tiling, a motif that is composed of two
different squares that share one side, and can be found all over the world in
kitchen or garage floors. Interestingly, it was proposed recently that
elementary, two-dimensional Cl, Br, and I might be able to adopt this
arrangement Huran  _et al._ (2021a). We also find a series of Cairo pentagonal
tilings. In this example, AgSe2 forms overlaying two tessellations of the
plane by irregular hexagons, where each of the hexagons is formed by four
identical pentagons. At the center of the hexagons we find a Se–Se bond. One
of the possible structures of Cu3S2 consists on a trihexagonal tiling (that is
often called the Kagome lattice due to its use in traditional Japanese
basketry) of the plane by Cu atoms, decorated by a S atom in the middle of the
triangles. Triakis triangular lattices appear quite commonly in our data. The
example in Fig. 5 can be seen as composed of Ce equilateral triangles
decorated with a Se atom at its center.

The following three examples are derived from snub-square lattices. In the
first two, the metal forms this interesting square-triangle lattice and the
non-metal decorates the squares. In the case of CoTe, Te-atoms can be found
above and below the plane of the Co atoms, while in Cu2Se Se-atoms alternate
above and below the plane of the Cu-atoms. We note that this specific lattice
was recently proposed for some noble metal chalcogenide monolayers Gao  _et
al._ (2021) and for certain Ba and Ti oxides Huran  _et al._ (2021b). The
third example is more complex, as both In and Tl form a distorted version of
the snub-square lattice, with further Tl-atoms alternating above and below the
plane. Note, however, that this curious structure is almost at the limit of
our energy threshold. Finally, we present six examples of rhombitrihexagonal
lattices, where the metal atoms form the triangle-square-hexagon lattice that
is then decorated (mostly) by the non-metals. We found a very diverse number
of different decorations, allowing for many different stoichiometries, ranging
from the simple AB and AB2 to the more unconventional A2B5 and A2B7. A very
interesting possibility raised by the finding of all these snub-square and
rhombitrihexagonal structures is that these can be easily inflated by a
recursive approach to generate quasi-crystalline systems Förster  _et al._
(2013, 2020). This is, however, only possible for structures not including
out-of-plane alternating atoms, as this induces frustration in the system
reducing its stability Huran  _et al._ (2021b).

Figure 6: Polymorphism in CuI. The vertical axis denotes increasing energy
distance to the convex hull of stability in meV/atom. Copper atoms are
depicted in gold, while iodine is in violet. For each structure we show a top
and a side view.

We have given several examples of the different stoichiometries in our dataset
and of the structural variety stemming from then. Now we look at the issue of
polymorphism, i.e., the different phases possible for a specific chemical
composition. Not surprisingly, polymorphism depends strongly on the chemical
elements present in the compound. For example, for BN we found a single
structure below 250 meV/atom, the well-known honeycomb lattice, while for
other compounds we have an extraordinary variety in the same energy range.

As an example we show, in Fig. 6, a selection of the crystal structures we
found for CuI. We recall that zincblende CuI is at the moment the most
promising $p$-type transparent conducting semiconductor Grundmann  _et al._
(2013). However, CuI has a number of polymorphs, including a a couple of
trigonal phases Kurdyumova and Baranova (1961); Sakuma (1988); Villars and
Cenzual (2009); Akopyan  _et al._ (2010) that are layered, with a bonding
pattern rather different from the $\gamma$-phase. At a quick glance, Fig. 6
shows us that also in the 2D case we find a large variety of structure and of
bonding patterns.

As the lowest-energy 2D layer we find a covalently bounded hexagonal double-
layer (i), that is essentially on the convex hull of thermodynamic stability.
The (buckled) single layer (x) and the van-der-Waals bounded double flat-layer
also appears in the energy spectrum but considerably higher at more than 70
meV/atom. The second most stable structure is surprisingly a rectangular
lattice of Cu–I (ii), with the I-atoms alternating above and below the plane
of the Cu-atoms, with a related lattice (iv) appearing just a few meV/atom
above. Structures (iii), that is only 6 meV/atom above the hull, and (ix) are
in fact arrangements of one-dimensional objects. The first exhibits nanowires
with a triangular section arranged in an alternating fashion as depicted in
Fig. 6. The latter (ix) is a periodic arrangements of nanostripes.
(Incidentally, higher in energy, at 161 meV/atom we even find a molecular
crystal of Cu4I4 pyramidal clusters). All these systems turn out to be
semiconducting, with calculated (PBE) band gaps ranging from around 0.5 eV to
more than 2.1 eV.

##  IV Prototype search

The largest advantage of the workflow presented above is that is (i)
systematic and (ii) unbiased in what concerns the structural variety.
Unfortunately, the price to pay for these advantages is efficiency, in the
sense that it is computationally expensive to go through all possible
compositions and space groups and that many of the possibilities turn out to
be highly unstable or lead to thick slabs. It is, however, possible to
accelerate considerably the exploration of the 2D material space by using the
structural prototypes discovered by our approach, and combining them with a
machine-learning model appropriate for prototype search Schmidt  _et al._
(2019); Kulik  _et al._ (2022). Of course, in this way we will not discover
new structural motifs, but we can explore the whole compositional space very
efficiently.

Our approach follows the same basic principles as V2DB Sorkun  _et al._
(2020), but goes beyond it in a number of different directions. First, we
perform transfer learning from a 3D machine, which allows to transfer many of
the chemical principles that govern atomic bonding. Second, we use a much
larger training set, increasing the accuracy of the machine. Third, we lift
several constrains (like charge neutrality or electronegativity rules) used in
V2DB and in our systematic search, and we expand the possible chemical
elements to the whole periodic table. This allows us to discover a variety of
intermetallics and compounds combining elements with unusual oxidation states.
We furthermore perform machine-learning predictions for all two-dimensional
prototypes, either already present in C2DB or stemming from our systematic
search. Finally, we perform validation DFT calculations for some of the
predictions, specifically for the binary stoichiometries A2B5, A2B7, and the
ternary ABC2, ABC3, AB2C2, A2B2C3.

To keep the number of structures manageable, we asked the machine to output
all structures that it found below 200 meV/atom for the binaries and 50
meV/atom for the ternaries. In total, we obtained 1023 candidates that were
pre-optimized with our neural-network force-field and then optimized with DFT.
From these 638 were found to be below 250 meV/atom from the hull, yielding an
exceptional success rate of around 62%. The lowest success rate, of only 9%,
was found for the A2B7 stoichiometry: as these compounds were sparsely present
in the training set, the machine could not learn the specificity of those
structures. The problem can be, of course, solved by adding further samples to
the dataset, in order to remove the structural (and compositional) bias, as
previously shown for bulk systems in Ref. Schmidt  _et al._ , 2022b. We are
currently performing DFT calculations for $\sim 40000$ more materials,
resulting from 238 million machine-learning predictions, that will be
available in the next release of our dataset.

##  V Conclusions

We have presented a systematic, exhaustive approach to explore the structural
and compositional diversity that is possible in the chemical space of 2D
materials. The main advantage of this approach is that it is not based on a
specific number of structural prototypes. This is particularly important for
2D materials, as the space of possible structures is still rather unexplored
and only few prototype structures, mainly from exfoliation of layered 3D
materials, are known so far. In this way we have discovered thousands of
unexpected phases that have no counterpart in the world of layered three-
dimensional materials. We expect that such unusual bonding patterns and
geometrical patterns will also lead to unique mechanical, electronic, optical,
and magnetic properties.

Our method relies heavily on the use of machine learning. We note that our
extensive use of neural networks in several parts of our workflow is self-
accelerating. In fact, the faster we generate more data for two-dimensional
systems, the larger will be our training sets, resulting in even more accurate
machine learning models. This leads to a virtuous cycle that, in our opinion,
will pave the way for a rather complete exploration of the binary, ternary,
and eventually also the quaternary two-dimensional phases in the near future.

Finally, an important question is how many of the phases in our dataset can be
synthesized. We chose to filter our results to include only compounds with an
energy less than 250 meV/atom above the convex hull, as these have higher
stability and therefore higher probability to be synthesized. (For comparison,
silicene, that has been experimentally synthesizedVogt  _et al._ (2012), is
more than 600 meV/atom above the hull in its free-standing form.) However,
besides these thermodynamic arguments, a key factor will be the choice of
suitable substrates that stabilize the two-dimensional layers, and, of course,
the ingenuity of experimental physicists to design targeted synthesis
strategies.

##  VI Appendix: Methods

###  VI.1 DFT calculations

We performed all geometry optimizations and total energy calculations with the
code vasp Kresse and Furthmüller (1996, 1996). The 2D Brillouin zones were
sampled by uniform $\Gamma$-centered k-point grids with a density of 6
$k$-points/Å-2. We performed spin-polarized calculations starting from a
ferromagnetic state, and used the projector augmented wave setups Blöchl
(1994); Kresse and Joubert (1999) of vasp version 5.2 with a cutoff of 520 eV.
We converged the calculations to forces smaller than 0.005 eV/Å. As exchange-
correlation functional we used the Perdew-Burke-Ernzerhof Perdew  _et al._
(1996) functional with on-site corrections for oxides, fluorides containing
Co, Cr, Fe, Mn, Mo, Ni, V, and W. The repulsive on-site corrections to the
d-states were respectively 3.32, 3.7, 5.3, 3.9, 4.38, 6.2, 3.25, and 6.2 eV.
These parameters were chosen to be compatible with the Materials Project
database Jain  _et al._ (2013). We imposed a vacuum region of 15 Å, and
systems that resulted in structures with a thickness greater than 7.5 Å were
automatically discarded. Finally, as it is common in this kind of approaches,
some of the calculations did not converge due to a multitude of reasons. The
corresponding phases were then simply eliminated from the dataset.

Distances to the convex hull were evaluated using pymatgen Ong  _et al._
(2013) using the large complex hull of Ref. Schmidt  _et al._ , 2022b
corresponding to the dataset available in the Materials Cloud repository
Schmidt  _et al._ (2022a).

###  VI.2 m3gnet

We employed the universal neural-network force-field m3gnet Chen and Ong
(2022) that was developed to reproduce the energies and the forces of bulk
structures with remarkable results. As a starting point we used the pretrained
network distributed with m3gnet. We tested this model on 1300 of our systems
by measuring the difference between the energy calculated with m3gnet at the
m3gnet relaxed structure and the energy calculated with DFT at the DFT relaxed
structure. We arrived at a mean absolute error (MAE) of 320 meV/atom and a
median absolute error of 223 meV/atom. These numbers are already rather small,
especially when we consider that the training set of m3gnet did not include 2D
systems, that can exhibit very different bonding patterns than their bulk
counterparts. As such, as soon as enough data was available from our own
simulations, we used transfer learning techniques to retrain m3gnet for 2D
materials (see Section VI). Specifically, we build a dataset comprising
energies, forces, and stresses from structures calculated during the geometry
optimization steps. Structures with extremely high forces above 50 eV/Å were
removed from the data as were structures with no neighboring atoms inside the
cutoff radius to avoid errors during training. To balance the training set,
for systems with more than 4 recorded geometry optimization steps only the
first, last and N${}_{\text{steps}}$/3 step were used. The final training set
for m3gnet contained 11612 geometry relaxations corresponding to 34944
energies and structures. The resulting network had a validation MAE of 61
meV/atom for direct energy predictions after training. The test errors for
geometry optimizations on the same dataset as the pretrained model were 198
meV/atom for the MAE and 96 meV/atom for the median absolute error proving the
efficiency of our transfer learning strategy. Of course, we expect these
errors to decrease further simply by adding more data to the training set. The
models were trained with the base hyperparameters from m3gnet and by setting
the loss function of the stress in the non-periodic direction to zero.

###  VI.3 Crystal-graph attention networks

We used the crystal-graph attention neural networks developed in Ref. Schmidt
_et al._ , 2021 as they were specifically crafted for prototype searches. In
particular, they only require as input the (unrelaxed) structural prototype
and not accurate relaxed structures. Of course, this model was trained on bulk
3D structures, so we do not expect it to perform accurately in our case.
However, many of the bonding patterns present in our 2D materials can already
be found in the 3D world. To take advantage of this, we performed transfer
learning of the 3D model, using the 2D structures in our dataset as training
data. We used a dataset of DFT calculations with 22007 entries, 80% of which
were used for training 10% for validation and 10% for testing. Evaluating both
models on the test set, we arrive at an MAE of 222 meV/atom for the original
model and 86 meV/atom for the model transferred to the 2D data.

Figure 7: Distribution of the distances to the convex hull (calculated with
DFT) for the compounds obtained through our systematic approach and stemming
from the machine-learning assisted prototype search compared to the entries of
C2DB.

In Fig. 7 we present a histogram with the distances to the convex hull of
stability calculated with DFT. The C2DB data is in light green, and is highly
peaked at zero, decaying slowly for larger energies. This is expected, as C2DB
was seeded with stable 3D, van der Waals bonded structures from ICSD. In blue
we depict the structures obtained through our systematic approach. These form
a continuous distribution with a peak at around 300 meV/atom, and extending
beyond 1.5 eV. Knowing that such distribution for random compounds can extend
beyond 4 eV, we see how the charge and electronegativity constrains lead to
relatively stable compounds (at the price of overlooking intermetallics or
compounds with unusual oxidation states. In light blue we show the machine-
learning binaries predicted to be within 200 meV/atom from the hull. This
shows a peak at around that value, as expected form the cutoff, then decaying
similarly to the C2DB data. The ternary entries, displayed in green, are
shifted to much lower energy, as expected by the smaller cutoff of 50
meV/atom. This results are absolutely consistent with the MAE of 86 meV/atom
for the 2D model. The CGAT-hyperparameters are listed in the supplementary
material and the code can be found at https://github.com/hyllios/CGAT.git.

##  VII Data availability statement

All relevant data used in or resulting from this work is available at
Materials Cloud (https://doi.org/10.24435/materialscloud:sb-cy). The database,
including structures, distances to the hull, and other basic properties, can
be accessed at https://tddft.org/bmg/physics/2D/ through a simple web-based
interface.

##  VIII Acknowledgements

We acknowledge the computational resources awarded by XSEDE, a project
supported by National Science Foundation grant number ACI-1053575. The authors
also acknowledge the support from the Texas Advances Computer Center (with the
Stampede2 and Bridges supercomputers). We also acknowledge the Super Computing
System (Thorny Flat) at WVU, which is funded in part by the National Science
Foundation (NSF) Major Research Instrumentation Program (MRI) Award #1726534,
and West Virginia University. AHR and LW were funded in part, by the
Luxembourg National Research Fund (FNR), Inter Mobility 2DOPMA, grant
reference 15627293. AHR also recognizes the support of West Virginia Research
under the call research challenge grand program. JS and MALM gratefully
acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for
funding this project by providing computing time on the GCS Supercomputer
SUPERMUC-NG at Leibniz Supercomputing Centre (https://www.lrz.de) under the
project pn25co. We also thank Kristian Thygesen for kindly providing us with
full access to the C2DB database.

##  IX Author Contributions

AHR performed the systematic structure generation; MALM and HCW performed the
DFT high-throughput calculations; JS performed the training of the machines
and the machine learning predictions of the distance to the hull; AHR, LW, and
MALM directed the research; all authors contributed to the analysis of the
results and to the writing of the manuscript.

##  X Competing Interests

The authors declare that they have no competing interests.

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D.-e. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Novoselov  _et al._ (2005) K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. V. Morozov, and A. K. Geim, Proc. Natl. Acad. Sci. U.S.A. 102, 10451 (2005).
  * Mounet  _et al._ (2018) N. Mounet, M. Gibertini, P. Schwaller, D. Campi, A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Cepellotti, G. Pizzi, and N. Marzari, Nat. Nanotechnol. 13, 246 (2018).
  * Sorkun  _et al._ (2020) M. C. Sorkun, S. Astruc, J. Koelman, and S. Er, Npj Comput. Mater. 6, 1 (2020).
  * Zhou  _et al._ (2019) J. Zhou, L. Shen, M. D. Costa, K. A. Persson, S. P. Ong, P. Huck, Y. Lu, X. Ma, Y. Chen, H. Tang, and Y. P. Feng, Sci. Data 6, 1 (2019).
  * Haastrup  _et al._ (2018) S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, J. Gath, K. W. Jacobsen, J. J. Mortensen, T. Olsen, and K. S. Thygesen, 2D Mater. 5, 042002 (2018).
  * Gjerding  _et al._ (2021) M. N. Gjerding, A. Taghizadeh, A. Rasmussen, S. Ali, F. Bertoldo, T. Deilmann, N. R. Knøsgaard, M. Kruse, A. H. Larsen, S. Manti, T. G. Pedersen, U. Petralanda, T. Skovhus, M. K. Svendsen, J. J. Mortensen, T. Olsen, and K. S. Thygesen, 2D Mater. 8, 044002 (2021).
  * Lyngby and Thygesen (2022) P. Lyngby and K. S. Thygesen, Npj Comput. Mater. 8 (2022).
  * Song  _et al._ (2021) Y. Song, E. M. D. Siriwardane, Y. Zhao, and J. Hu, ACS Appl. Mater. Interfaces 13, 53303 (2021).
  * Jain  _et al._ (2013) A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, APL Mater. 1, 011002 (2013).
  * Gražulis  _et al._ (2009) S. Gražulis, D. Chateigner, R. T. Downs, A. Yokochi, M. Quirós, L. Lutterotti, E. Manakova, J. Butkus, P. Moeck, and A. Le Bail, J. Appl. Crystallogr. 42, 726 (2009).
  * Allen (2002) F. H. Allen, Acta Crystallogr., Sect. B: Struct. Sci. 58, 380 (2002).
  * Van Hove  _et al._ (2002) M. A. Van Hove, K. Hermann, and P. Watson, Acta Crystallogr., Sect. B: Struct. Sci. 58, 338 (2002).
  * Saal  _et al._ (2013) J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, and C. Wolverton, Jom 65, 1501 (2013).
  * Curtarolo  _et al._ (2012) S. Curtarolo, W. Setyawan, G. L. 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. Demchenko, and D. Morgan, Comput. Mater. Sci. 58, 218 (2012).
  * Talirz  _et al._ (2020) L. Talirz, S. Kumbhar, E. Passaro, A. V. Yakutovich, V. Granata, F. Gargiulo, M. Borelli, M. Uhrin, S. P. Huber, S. Zoupanos, C. S. Adorf, C. W. Andersen, O. Schütt, C. A. Pignedoli, D. Passerone, J. VandeVondele, T. C. Schulthess, B. Smit, G. Pizzi, and N. Marzari, Sci. Data 7, 1 (2020).
  * Bergerhoff and Brown (1987) G. Bergerhoff and I. Brown,  _Crystallographic Databases_ , edited by F. Allen, G. Bergerhoff, and R. Sievers, Vol. 360 (International Union of Crystallography, Chester, 1987) pp. 77–95.
  * Zagorac  _et al._ (2019) D. Zagorac, H. Müller, S. Ruehl, J. Zagorac, and S. Rehme, J. Appl. Crystallogr. 52, 918 (2019).
  * Belsky  _et al._ (2002) A. Belsky, M. Hellenbrandt, V. L. Karen, and P. Luksch, Acta Crystallogr., Sect. B: Struct. Sci. 58, 364 (2002).
  * Schmidt  _et al._ (2022a) J. Schmidt, N. Hoffmann, H.-C. Wang, P. Borlido, P. J. M. A. Carriço, T. F. T. Cerqueira, S. Botti, and M. A. L. Marques, “Large-scale machine-learning-assisted exploration of the whole materials space,”  (2022a).
  * Park and Wolverton (2020) C. W. Park and C. Wolverton, Phys. Rev. Mater. 4, 063801 (2020).
  * Schmidt  _et al._ (2022b) J. Schmidt, N. Hoffmann, H.-C. Wang, P. Borlido, P. J. M. A. Carriço, T. F. T. Cerqueira, S. Botti, and M. A. L. Marques, “Large-scale machine-learning-assisted exploration of the whole materials space,” (2022b), arXiv:2210.00579 .
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Litvin and Kopsky (2002) D. Litvin and V. Kopsky, “International tables for crystallography. volume d: Subperiodic groups,” (2002).
  * Aroyo  _et al._ (2011) M. I. Aroyo, J. M. Perez-Mato, D. Orobengoa, E. Tasci, G. de la Flor, and A. Kirov, Bulg. Chem. Commun 43, 183 (2011).
  * Aroyo  _et al._ (2006) M. I. Aroyo, A. Kirov, C. Capillas, J. Perez-Mato, and H. Wondratschek, Acta Crystallogr., Sect. A: Found. Crystallogr. 62, 115 (2006).
  * Fredericks  _et al._ (2021) S. Fredericks, K. Parrish, D. Sayre, and Q. Zhu, Comput. Phys. Commun. 261, 107810 (2021).
  * Davies  _et al._ (2016) D. W. Davies, K. T. Butler, A. J. Jackson, A. Morris, J. M. Frost, J. M. Skelton, and A. Walsh, Chem 1, 617 (2016).
  * Chen and Ong (2022) C. Chen and S. P. Ong, arXiv preprint arXiv:2202.02450 (2022).
  * Restrepo  _et al._ (2004) G. Restrepo, H. Mesa, E. J. Llanos, and J. L. Villaveces, J. Chem. Inf. Model. 44, 68 (2004).
  * Huran  _et al._ (2021a) A. W. Huran, H.-C. Wang, A. San-Miguel, and M. A. L. Marques, J. Phys. Chem. Lett. 12, 4972 (2021a).
  * Gao  _et al._ (2021) L. Gao, Y.-F. Zhang, and S. Du, Nano Res. 14, 2826 (2021).
  * Huran  _et al._ (2021b) A. W. Huran, H.-C. Wang, and M. A. L. Marques, 2D Mater. 8, 045002 (2021b).
  * Förster  _et al._ (2013) S. Förster, K. Meinel, R. Hammer, M. Trautmann, and W. Widdra, Nature 502, 215 (2013).
  * Förster  _et al._ (2020) S. Förster, S. Schenk, E. M. Zollner, O. Krahn, C.-T. Chiang, F. O. Schumann, A. Bayat, K.-M. Schindler, M. Trautmann, R. Hammer, K. Meinel, W. A. Adeagbo, W. Hergert, J. I. Flege, J. Falta, M. Ellguth, C. Tusche, M. DeBoissieu, M. Muntwiler, T. Greber, and W. Widdra, Phys. Status Solidi B 257, 1900624 (2020).
  * Grundmann  _et al._ (2013) M. Grundmann, F.-L. Schein, M. Lorenz, T. Böntgen, J. Lenzner, and H. von Wenckstern, Phys. Status Solidi A 210, 1671 (2013).
  * Kurdyumova and Baranova (1961) R. Kurdyumova and R. Baranova, Kristallografiya 6, 402 (1961).
  * Sakuma (1988) T. Sakuma, J. Phys. Soc. Jpn. 57, 565 (1988).
  * Villars and Cenzual (2009) P. Villars and K. Cenzual, eds.,  _Landolt-Börnstein - Group III Condensed Matter: Crystal Structures of Inorganic Compounds_ (Springer-Verlag Berlin Heidelberg, Heidelberg, 2009).
  * Akopyan  _et al._ (2010) I. K. Akopyan, V. V. Golubkov, O. A. Dyatlova, A. N. Mamaev, B. V. Novikov, and A. N. Tsagan-Mandzhiev, Phys. Solid State 52, 805 (2010).
  * Schmidt  _et al._ (2019) J. Schmidt, M. R. G. Marques, S. Botti, and M. A. L. Marques, Npj Comput. Mater. 5, (2019).
  * Kulik  _et al._ (2022) H. J. Kulik, T. Hammerschmidt, J. Schmidt, S. Botti, M. A. L. Marques, M. Boley, M. Scheffler, M. Todorović, P. Rinke, C. Oses, A. Smolyanyuk, S. Curtarolo, A. Tkatchenko, A. P. Bartók, S. Manzhos, M. Ihara, T. Carrington, J. Behler, O. Isayev, M. Veit, A. Grisafi, J. Nigam, M. Ceriotti, K. T. Schütt, J. Westermayr, M. Gastegger, R. J. Maurer, B. Kalita, K. Burke, R. Nagai, R. Akashi, O. Sugino, J. Hermann, F. Noé, S. Pilati, C. Draxl, M. Kuban, S. Rigamonti, M. Scheidgen, M. Esters, D. Hicks, C. Toher, P. V. Balachandran, I. Tamblyn, S. Whitelam, C. Bellinger, and L. M. Ghiringhelli, Electron. Struct. 4, 023004 (2022).
  * Vogt  _et al._ (2012) P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Phys. Rev. Lett. 108, 155501 (2012).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Ong  _et al._ (2013) S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, Comput. Mater. Sci. 68, 314 (2013).
  * Schmidt  _et al._ (2021) J. Schmidt, L. Pettersson, C. Verdozzi, S. Botti, and M. A. L. Marques, Sci. Adv. 7, (2021).


11institutetext:  School of Physics and Information Technology, Shaanxi Normal
University, Xi’an 710062, China
School of Physics and Electronic-Electrical Engineering, Ningxia University,
Yinchuan 750021, P. R. China
Beijing Advanced Innovation Center for Big Data and Brain Computing, Beihang
University, Beijing 100191, P. R. China
Robert Koch-Institute, Nordufer 20, 13353 Berlin, Germany

# Emergent route towards cooperation in interacting games: the dynamical
reciprocity

Qin-Qin Wang 11 Ri-Zhou Liang 11 Ji-Qiang Zhang 2233 Guo-Zhong Zheng 11 Lin Ma
11 Li Chen<EMAIL_ADDRESS>
###### Abstract

The success of modern civilization is built upon widespread cooperation in
human society, deciphering the mechanisms behind has being a major goal for
centuries. A crucial fact is, however, largely missing in most prior studies
that games in the real world are typically played simultaneously and
interactively rather than separately as assumed. Here we introduce the idea of
interacting games that different games coevolve and influence each other’s
decision-making. We show that as the game-game interaction becomes important,
the cooperation phase transition dramatically improves, a fairly high level of
cooperation is reached for all involved games when interaction goes to be
strong. A mean-field theory indicates that a new mechanism — _the dynamical
reciprocity_ , as a counterpart to the well-known network reciprocity, is at
work to foster cooperation, which is confirmed by the detailed analysis. This
revealed reciprocity is robust against variations in the game type, the
population structure, and the updating rules etc, and more games generally
yield a higher level of cooperation. Our findings point out the great
potential towards high cooperation for many issues are interwoven with each
other in the real world, and also the possibility of sustaining decent
cooperation even in extremely adverse circumstances.

##  1 Introduction

Recent withdrawals of the United States from a couple of “groups” like WHO,
Paris Agreement, UNESCO etc. signifies a degraded cooperation at the global
scale. Any solution to this sort of problems requires an understanding of what
processes drive and maintain human cooperation and what measures or
institutions could be implemented for its promotion. The key question to be
addressed is: why entities help each other who could potentially be in
competition and incur a cost to themselves? As the paradigm of _homo
economicus_ shows, people always try to maximize their earnings and avoid
irrational investments, which inevitably leads to the tragedy of the commons
[1].

Important progresses have been made with the help of evolutionary game theory
[2] by analysing the stylized social dilemmas such as prisoner’s dilemma and
the public goods game. Several mechanisms are proposed [3] in the past several
decades, such as reward and punishment [4], social diversity [5], direct [6]
or indirect reciprocity [7], kin [8] or group selection [9, 10], spatial or
network reciprocity [11]. In particular, theoretically accounting for the fact
that human populations are highly organized and individuals interact
repeatedly with their immediate neighbors can support cooperation [11]. The
rationale behind is that a structured neighborhood facilitates the formation
of cooperator clusters, which effectively resist the invasion of defectors, as
opposed to the well-mixed scenario. The ensuing years have witnessed a wealth
of theoretical studies that further confirm this so-called network reciprocity
for various population structures [12]. However, recent human behavioral
experiments show that structured populations do not promote cooperation in
general [13, 14], at least some conditions combining the game parameters and
the population structure must be met for cooperation to thrive [15]. One
explanation is that the complexities of human psychology make humans switch
strategies frequently that the assortment fails in the static networks [16].
But dynamic networks indeed offer an escape because the players are allowed to
adjust social ties and they are more like to cooperate under this peer
pressure [17, 18]. This unsatisfactory situation implies that some essential
elements could be missing in current game-theoretic models and the experiment-
driven modeling approach is needed. Note that, in most of these studies a
single game is considered, and they focus on the factors of interest like the
underlying structures of population, or the impact of some dynamical
processes, e.g. the punishment or reward; the conclusions drawn are supposed
to be applicable in more general circumstances.

Games, however, may not be unfolded in isolation but often in parallel. For
instance, we humans are engaged in different activities, works, sports, and
recreations; colleagues in a company could work on a couple of concurrent
projects; and countries have to deal with a whole range of conflicts such as
trade war, security issues, diplomatic crisis etc. Only when the evolution of
these games is independent from each other, the modeling efforts based on a
single game are then reasonable as most of existing work assumed. Observations
in aforementioned contexts, however, suggest that the decision-making of
entities in one game is often conditioned by what happened in another. Similar
observations are also made in biological games, like chimpanzees are more
likely to groom their fellows if they are skillful in hunting, and vice versa,
and actually these two behaviors together with sharing food, joint patrol the
borders, support one another in conflicts etc are all correlated in the
chimpanzee’s social life [19].

A closely related research line is the multigame dynamics, the existing work
shows that dynamical inconsistencies are already possible when two or more
non-repeated games are coupled [20, 21, 22, 23], meaning that the eventual
fate of games cannot be inferred from the single game dynamics. A more recent
work starts to study the repeated scenario and an evolutionary framework of
the so-called multichannel games is proposed [24], where they find that the
fixed game linkage is able to enhance cooperation in all games engaged in
general. Still, fundamental questions remain: _what typical evolutionary
dynamics are expected when more games are engaged, to what extent would such
game-game interaction alter the classic cooperation mode of single game, and
any new cooperation scenario arises therein?_

In this Letter, we mainly study two symmetrically interacting games, where
they have a stake in each other, and focus on clarifying the impact of game-
game interaction. We reveal a new type of reciprocity rooted in the game-game
interaction that is able to maintain high levels of cooperation. In
particular, fairly high cooperation is expected when the interaction goes to
the extreme that the decision-making of a given game is completely conditioned
by the other and vice versa. The mechanism behind lies in the new types of
interactions that lead to a persistent advantage of cooperators. Furthermore,
the uncovered reciprocity is found to be quite robust and more games generally
lead to be more cooperative for all games. For a detailed account of this
study, see [26].

Figure 1:  The cooperation evolution of two symmetrically interacting PD on
the 2d square lattice. (a) Phase transitions of cooperation prevalence
regarding game $G_{1}$ ($f^{G_{1}}_{C}=f_{CC}+f_{CD}$) versus the temptation
$b$ is shown for interaction strengths $\theta=0,0.25,0.5,0.75$ and 1. Due to
the symmetry, $f^{G_{2}}_{C}\approx f^{G_{1}}_{C}$ (data not shown for visual
clarity). (b) Typical time series are shown with fixed $b=1.1$. Parameters:
$L=1024$ and $K=0.1$, the random initial condition for both games, data is
obtained over 50 ensemble averages after transient in (a).

##  2 Modeling two interacting games

Suppose that two games $\mathbb{G}\\!=\\!\\{G_{1},G_{2}\\}$ are played
simultaneously in a population composed of $N$ players, where they are located
on an $L\times L$ square lattice with a periodic boundary condition. They can
adopt one of the two strategies for each game: cooperation ($C$) or defection
($D$), i.e. $\mathbb{S}_{1}=\\{C,D\\}$. Therefore, there are four possible
states $\mathbb{S}_{2}=\\{XY|CC,CD,DC,DD\\}$ in the two interacting games,
where $X,Y$ represent the state regarding game $G_{1,2}$ respectively. For
simplicity, we resort to the pairwise game defined as follows: mutual
cooperation brings both a reward $R$, mutual defection leads to a punishment
$P$ for each, and mixed encounter yields the cooperator a sucker’s payoff $S$
yet a temptation $T$ for the defector. Their ranking determines the game type.
Here, we follow the common practice for a weak prisoner’s dilemma (PD) with
$R=1$, $P=S=0$, $T=b>1$ for both games if not stated otherwise.

Following the standard Monte Carlo (MC) simulation procedure, firstly a game
$g\in\mathbb{G}$ is chosen at random to play in an elementary step, a player
$i$ is then randomly chosen and accumulates its payoff $\Pi_{i}$. Next, one of
$i$’s neighbors $j$ is picked randomly, and acquires its payoff $\Pi_{j}$ as
well. Lastly, player $i$ adopts $j$’s strategy regarding game $g$ with a
probability according to the Fermi rule [25]

$W^{g}_{j\rightarrow
i}=\frac{1}{1+\exp[(\widehat{\Pi}_{i}^{g}-\widehat{\Pi}_{j}^{g})/K]},$ (1)
$\widehat{\Pi}_{i,j}^{G_{1,2}}=(1-\theta)\Pi_{i,j}^{G_{1,2}}+\theta\Pi_{i,j}^{G_{2,1}},$
(2)

where $\widehat{\Pi}_{i,j}^{g}$ is _the effective payoffs_ , which captures
the reality that to imitate, players compare the overall payoff profiles in
all games rather than simply the one under play. Therefore the decision is
made based upon a combination of both payoffs, a simple case is as shown in
Eq. (2). We interpret the weight $\theta\in[0,1]$ as the game interaction
strength, a larger value means a stronger impact of the other game; two
extreme cases $\theta=0,1$ correspond to the two independent games and the
cross-playing scenario, respectively. $K$ is a temperature-like parameter,
measuring the uncertainties in the imitation process, its inverse can be
interpreted as the selection intensity in biology or the bounded rationality
in economical contexts. A full MC step consists of
$2\\!\times\\!L\\!\times\\!L$ elementary steps, where every player is updated
once for each game on average. Simulations are carried out for $L=1024$, and
the data for the cooperator fractions are averaged over $10^{6}$ MC steps
after a transient period of $10^{6}$ steps.

##  3 Results

Varying the game interaction strength $\theta$, we observe a continuing
promotion in the cooperation prevalence $f_{c}$ as a function of the
temptation $b$ [see Fig. 1(a)]. For the independent game case where
$\theta=0$, a second-order phase transition (PT) for cooperation is seen but
the cooperation region is rather small with the critical temptation
$b_{c}\approx 1.038$, beyond which the cooperators become extinct. As $\theta$
is increased, $b_{c}$ is shifted to the right, the prevalence also becomes
higher, the cooperation is lifted. Finally, as $\theta\rightarrow 1$, this
promotion is maximal, where the PTs become absent and nearly full cooperation
is seen across the whole parameter region $1\leq b\leq 2$ for both games. This
is quite unexpected since in the cross-playing scenario ($\theta=1$), the
decision-making of a game is entirely blind to its own payoff. This
observation of promotion is strengthened by the corresponding time series by
fixing $b=1.1$ shown in Fig. 1(b), where the initial decrease in $f_{c}$ is
also inhibited when $\theta$ becomes large.

Figure 2:  The evolution of cooperation patterns for $\theta=0$ (top row), 0.5
(middle row), and 1(bottom row). The system is prepared with full cooperators
CC within upper half domain versus full defectors DD within another half (the
leftmost panel). The characteristic snapshots are taken at
$t=0,50,100,200,300$ for $\theta=0$, whereas $t=0,100,1000,5000,10000$ for the
others. Parameters: the domain is of size $128\times 128$, $K=0.1$, and
$b=1.1$.

To gain some intuition of how the game-game interaction affects cooperation,
we first look at how the spatiotemporal evolution is influenced. Fig. 2 shows
the cases of $\theta=0,0.5,1$, but starting from a bulk initial condition
because it is more intuitive without altering $f_{c}$ compared to the random
initial condition cases. When $\theta=0$, defectors dominate in both games,
$DD$ players invade the $CC$’s domain, and cooperators quickly go extinct. At
the intermediate strength $\theta=0.5$, this advantage disappears where all
four species coexist. In the other extreme $\theta=1$, a reversed invasion is
seen where $CC$ players dominate and take over the whole domain in the end.
This suggests that a reversed advantage is expected between cooperators and
defectors as the game-game interaction is engaged.

##  4 A mean-field theory

To understand the rationale behind, we develop a mean-field theory based on
the replicator equation [27, 28], where the evolution of the four fractions
with respect to each game depends on their relative fitness measured by the
payoffs that can be formally described as
$\dot{f}_{s}=f_{s}(\widehat{\Pi}_{s}-\bar{\Pi})$, where
$g\\!\in\\!\mathbb{G}$, $s\\!\in\\!\mathbb{S}_{2}$, and
$\bar{\Pi}=\sum_{s}f_{s}\widehat{\Pi}_{s}$ is the average fitness. With some
algebra (see [26]), we obtain the ordinary differential equations of
cooperator fraction for game $G_{1}$ (i.e. $f^{G_{1}}_{C}=f_{CC}+f_{CD}$, the
exchange of 1 and 2 applies for game $G_{2}$) as

$\dot{f}^{G_{1}}_{C}=f^{G_{1}}_{C}f^{G_{1}}_{D}(\Pi^{G_{1}}_{C}-\Pi^{G_{1}}_{D})+(f_{CC}f_{DD}-f_{CD}f_{DC})(\Pi^{G_{2}}_{C}-\Pi^{G_{2}}_{D}),$
(3)

where $\Pi^{G_{1,2}}_{C,D}$ are the fitness in game $G_{1}$ or $G_{2}$
respectively for the cooperators and defectors. The first term in the rhs. is
well-known in the single game scenario [29] that comes from the game under
play, meaning that the fitness advantage in cooperators
$\Pi^{G_{1}}_{C}>\Pi^{G_{1}}_{D}$ converts the defectors into cooperators when
they meet up. The second term is new that captures the game-game interaction.
Specifically, the impact of the other game is through two interacting pairs:
i) when $CC$ players come across $DD$, the advantage of cooperators in game
$G_{2}$ ($\Pi^{G_{2}}_{C}>\Pi^{G_{2}}_{D}$) also facilitates the proliferation
of cooperators in game $G_{1}$ due to the game-game correlation; ii)
unexpectedly, in the opposite case when $\Pi^{G_{2}}_{D}>\Pi^{G_{2}}_{C}$, the
advantage of defectors in $G_{2}$ also helps the growth of cooperators in
$G_{1}$ when $CD$ encounters $DC$ players. Therefore, the above analysis shows
that potentially there are now new dynamical routes at work towards
cooperation in addition to the one in the single game case.

|  Individual scenario |  Bulk scenario |
---|---|---|---
Invasion |  | CC + DD $\xrightarrow{G_{1}/G_{2}}$ DC/CD + DD
---
| CC + DD $\xrightarrow{G_{1}/G_{2}}$ CC + CD/DC
---
Neutral |  | CC + DC $\xrightarrow{G_{1}}$ 2CC or 2DC
---
CC + CD $\xrightarrow{G_{2}}$ 2CC or 2CD
DD + CD $\xrightarrow{G_{1}}$ 2DD or 2CD
DD + DC $\xrightarrow{G_{2}}$ 2DD or 2CD
| CC + DC $\xrightarrow{G_{1}}$ 2CC or 2DC
---
CC + CD $\xrightarrow{G_{2}}$ 2CC or 2CD
DD + CD $\xrightarrow{G_{1}}$ 2DD or 2CD
DD + DC $\xrightarrow{G_{2}}$ 2DD or 2CD
Catalyzed |  | CD + DC $\xrightarrow{G_{1}/G_{2}}$ CD/DC + CC
---
| CD + DC $\xrightarrow{G_{1}/G_{2}}$ DC/CD + DD
---
Table 1: Classification of interactions in two interacting PD games for the
cross-playing scenario ($\theta=1$), where all six pairwise interactions can
be classified into three categories, in either individual or bulk scenario. In
invasion interactions, the payoff advantage of a given strategy in either game
is explicitly transformed into its reproduction. For neutral type, the state
change is purely random, no net conversion towards cooperation or defection is
expected. In the last category, the advantage of defectors/cooperators
regarding a given game cross-catalyzes the production of cooperators/defectors
in the other game. The last two categories are only possible in interacting
games. Figure 3:  Time evolution of all six interface proportions for two
interacting PD games starting with random initial conditions for $\theta=1$
(a) and $0.5$ (b), respectively. Inset show the relative fraction of the two
non-neutral types $P_{r}=P_{{}_{CD-DC}}/P_{{}_{CC-DD}}$. Parameters: $b=1.1$
and $L=1024$ for the 2d square lattice.

##  5 Mechanism analysis

To be more specific, it’s helpful to look into all six interactions in details
in our lattice system, as listed in Table 1. Here, we distinguish two
scenarios — individual and bulk scenarios. In the former, we only focus on the
evolution of the interaction pairs when without knowledge of their
surroundings such as the random state configuration. The bulk scenarios apply
for the circumstance when players of the same type are well-bulked, both
intra- and inter-bulk play are incorporated. The two scenarios are typically
present in the early phase of evolution and afterwards, respectively. For
simplicity, we consider the cross-playing case, where the six pairs of
interactions can be classified into three categories: invasion, neutral, and
catalyzed type for both scenarios (Table 1). While the neutral type of
interactions brings no net effect on cooperation, the other two categories
determine the evolution of cooperation prevalence, though they always have the
opposite effects either in individual or bulk scenario.

Typical evolution of all interactions is shown in Fig. 3(a) starting from
random initial conditions, where the evolution can be roughly divided into two
stages. (i) At the early stage $t\\!<\\!t_{c}$ ($t_{c}\\!\approx\\!10$ MC
steps) when no clear clusters are formed and thus the individual scenario
applies, the proportion $P_{{}_{CD-DC}}>P_{{}_{CC-DD}}$ is detected, meaning
that catalyzed interactions dominate over the invasion ones; and according to
the evolutionary dynamics in Table 1, a net production of cooperation is
expected. (ii) When $t>t_{c}$, clusters are gradually formed, where both the
size and compactness increase (see [26]), therefore the bulk scenario sets in.
Interestingly, a crossover is seen that $P_{{}_{CD-DC}}<P_{{}_{CC-DD}}$, the
reversed dominance again yields a net increase of cooperators since
$CC\\!-\\!DD$ pairs in bulk scenario favor the cooperators (Table 1).
Therefore, cooperation is preferred in the whole evolutionary processes. Back
to the mean-field equation Eq.(3), our analysis indicates that the second term
always brings a positive contribution to the cooperation evolution.
Thereinafter, we term this mechanism caused by game-game interactions as the
_dynamical reciprocity_. It also works for cases with $\theta<1$, as shown in
Fig. 3(b). However, the dynamical reciprocity only works in structured
population, no promotion is seen in the well-mixed population (see [26]).

##  6 Robustness

Within structured population, the revealed reciprocity is quite robust. In
[26], we show that when the interacting game is extended to be general
pairwise games (including snowdrift game and stage hunt etc), a similar
cooperation promotion is still observed irrespective of the game type, and
fairly high cooperation is expected for the whole parameter domain when games
are cross played. We also show it is also applicable to a multiplayer game
(the public goods game); Robust observations are made in model variants such
as asymmetrically interacting games, games with different updating rules (like
replicator rule, Moran rule, follow-the-best rule etc [28]), with different
time-scales, and even with two different games, i.e. a PD is coupled with a
snowdrift game.

Additional structural complexities [26] from underlying populations like
small-world networks and Erdős-Rényi random topologies also do not change the
working of the reciprocity. The structural heterogeneity neither alters the
promotion trend, as shown in the case of scale-free networks.

In particular, when the number of engaged games increases, a higher level of
cooperation is expected in general, see Fig. 4, where equal contribution for
each game is assumed when more than one game is played simultaneously. Since
potentially there are many issues interwoven with each other in reality, much
higher cooperation is expected than the case when only a single game is
unfolded.

Figure 4:  (a) Phase transitions of cooperation prevalence for one-, two-, and
three-game cases, where the each game is of equal contribution in the
effective payoffs for the later two cases (i.e. $\theta=1/2$ and $1/3$,
respectively). (b) Time series for fixed $b=1.1$. Parameters: $L=1024$ for the
2d square lattice and $K=0.1$.

##  7 Conclusions

In summary, the discussed game-game interaction is a natural ingredient that
may underpin a wealth of issues, from complex behaviors in animals, to inter-
personal activities in daily life, and even to international relationships at
the global scale. The potential for being highly cooperative, as revealed here
points out a promising route towards a cooperative world. It is worthwhile to
emphasize that contrary to the network reciprocity, where the underlying
structure of population plays the key role [11, 12], including the
interdependent network reciprocity [30, 31, 32]. The mechanism behind the
promotion here stems instead from the dynamical interaction among different
games. Our results suggest that the dynamical reciprocity could constitute a
new category of mechanisms behind the emergence of cooperation.

On the theoretic side, our finding of “more is different” [33] calls for more
systematic investigations in specific contexts, since the revealed mechanism
may offer valuable inspiration to avoid the cooperation crises. On the
experimental side, behavioral experiments are needed to justice the dynamical
reciprocity in realities and unveil other complexities that may arise in
interacting games.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
under Grants 61703257 and 12075144, and by the Fundamental Research Funds for
the Central Universities GK201903012. L. C. acknowledges the enlightening
discussions with Dirk Brockmann (HU and RKI) in the early phase of the project
and Ying-Cheng Lai (ASU) for helpful comments.

## References

  * [1] Hardin G. Science16219681243.
  * [2] Nowak M. A. Karl S. Science3032004793.
  * [3] Nowak M. A. Science31420061560.
  * [4] Sigmund K., Hauert C. Nowak M. A. Proc. Natl. Acad. Sci. U.S.A.98200110757.
  * [5] Santos F. C., Santos M. D. Pacheco J. M. Nature45420082131.
  * [6] Trivers R. L. Q Rev Biol46197135.
  * [7] Nowak M. A. Sigmund K. Nature3931998573.
  * [8] Hamilton W. D. J. Theor. Biol.7196417.
  * [9] Keller L. Levels of selection in evolution Princeton University Press, Princeton, New Jersey 1999\.
  * [10] Queller D. C. Nature20119641145.
  * [11] Nowak M. A. May R. M. Nature3591992826.
  * [12] Szabó G. Fath G. Phys. Rep.446200797.
  * [13] Carlos G.-L., Cuesta J. A., Sánchez A. Moreno Y. Sci. Rep.22012325.
  * [14] Carlos G.-L., Alfredo F., Gonzalo R., Alfonso T., Cuesta J. A., Angel S. Yamir M. Proc. Natl. Acad. Sci. U.S.A.109201212922.
  * [15] Rand D. G., Nowak M. A., Fowler J. H. Christakis N. A. Proc. Natl. Acad. Sci. U.S.A.111201417093.
  * [16] Traulsen A., Semmann D., Sommerfeld R. D., Krambeck H. J. Milinski M. Proc. Natl. Acad. Sci. U.S.A.10720102962.
  * [17] Rand D. G., Samuel A. and Christakis N. A. Proc. Natl. Acad. Sci. U.S.A.108201119193.
  * [18] Katrin F., Post D. J., Van D. Dirk S. Ecol. Lett.142011546.
  * [19] Hammerstein P. et al. Genetic and cultural evolution of cooperation MIT press, Cambridge,Massachusetts 2003\.
  * [20] Cressman R., Gaunersdorfer A. Wen J.-F. Int. J. Game Theory2200067.
  * [21] Chamberland M. Cressman R. Games Econ Behav302000319.
  * [22] Hashimoto K. J. Theor. Biol.2412006669.
  * [23] Venkateswaran V. R. Gokhale C. S. Proc. Royal Soc. B286201920190900.
  * [24] Donahue K., Hauser O. P., Nowak M. A. Hilbe C. Nat. Commun.1120201.
  * [25] Szabó G. Tőke C. Phys. Rev. E58199869.
  * [26] Liang, Rizhou et al. arXiv:2102.003602021.
  * [27] Taylor P. D. Jonker L. B. Math Biosci401978145.
  * [28] Roca C. P., Cuesta J. A. Sánchez A. Phys. Life Rev.62009208.
  * [29] Smith J. M. Evolution and the Theory of games Cambridge University Press, Cambridge,England 1982\.
  * [30] Wang Z., Szolnoki A. Matjaž P. Scientific Reports320131183.
  * [31] Jin Q., Wang L, Xia C. Y. Wang Z. Scientific Reports420144095.
  * [32] Xia C. Y., Meng X. Kun Wang Z PLOS ONE10(6)2015e0129542.
  * [33] Anderson P. W. Science1771972393.