# Surface Rearrangement and Evaporation Kinetics of Supported Gold
Nanoparticle Catalysts

James P. Horwath  University of Pennsylvania, Department of Materials Science
and Engineering, Philadelphia, PA 19104  Colin Lehman-Chong  University of
Pennsylvania, Department of Chemical and Biomolecular Engineering,
Philadelphia, PA 19104  Aleksandra Vojvodic  University of Pennsylvania,
Department of Chemical and Biomolecular Engineering, Philadelphia, PA 19104
<EMAIL_ADDRESS>Eric A. Stach  University of Pennsylvania, Department of
Materials Science and Engineering, Philadelphia, PA 19104
<EMAIL_ADDRESS>
###### Abstract

Heterogeneous catalysts consisting of supported metallic nanoparticles
typically derive exceptional catalytic activity from their large proportion of
under-coordinated surface sites which promote adsorption of reactant
molecules. Simultaneously, these high energy surface configurations are
unstable, leading to nanoparticle growth or degradation, and eventually a loss
of catalytic activity. Surface morphology of catalytic nanoparticles is
paramount to catalytic activity, selectivity, as well as degradation rates,
however, it is well-known that harsh reaction conditions can cause the surface
structure to change. Still, limited research has focused on understanding the
link between nanoparticle surface facets and degradation rates or mechanisms.
Here, we study a model Au supported catalyst system over a range of
temperatures using a combination of in situ Transmission Electron Microscopy,
kinetic Monte Carlo simulations, and density functional theory calculations to
establish an atomistic picture of how variations in surface structures and
atomic coordination environments lead to shifting evolution mechanisms as a
function of temperature. By combining experimental results, which yield direct
observation of dynamic shape changes and particle evaporation rates, with
computational techniques, which enable understanding the fundamental
thermodynamics and kinetics of nanoparticle evolution, we illustrate a two-
step evolution mechanism in which mobile adatoms form through desorption from
low-coordination facets and subsequently evaporate off the particle surface.
By understanding the role of temperature in the competition between surface
diffusion and evaporation, we are able to show how individual atomic movements
lead to particle-scale morphological changes, and rationalize why evaporation
rates vary between particles in a system of nearly identical nanoparticles.

Supported metal nanoparticles (NPs) are commonly used in industrial catalysis.
As heterogeneous catalysts, under-coordinated active sites on the surfaces of
nanoparticles promote adsorption of reaction gas molecules and facilitate
chemical reactions to enhance both catalytic activity and selectivity.
However, the high surface free energy associated with these under-coordinated
surfaces and instability of under-coordinated active sites simultaneously
drive nanoparticle evolution to reduce the total free energy of the system by
reducing the catalyst surface area. Maintaining catalytic activity for
important commercial reactions therefore commonly relies on preventing this
coarsening, and keeping nanoparticles small1. Research focusing on coarsening
and degradation mechanisms in supported nanoparticles has focused on the use
of mean-field theories for decades. While this approach is helpful for
developing an analytical, intuitive model which can be compared against simple
experiments, several key assumptions limit the use of these models for
specific applications and real systems with large complexity. Mean-field
theories which describe diffusional interactions between particles and their
surrounding field rely on the assumption of infinite separation between
particles to avoid accounting for direct interactions between particles, and
disregard any surface energy anisotropy by studying the behavior of a smooth,
(truncated) spherical particle.2, 3. However, these assumptions are not valid
in real systems with complex, heterogeneous and anisotropic particles which
are the key features of these systems responsible for their unique chemical
and catalytic behaviour.

In reality, it is well known that nanoscale particles are commonly faceted and
that the surface energy of each facet is not the same and can change with
environmental conditions e.g. temperature, oxidizing or reducing conditions,
reacting gases4. Therefore, understanding the competition between degradation
mechanisms in model catalysts requires moving beyond mean-field theories and
focusing on local changes which cause deviations from the expected mean-field
behavior. In our previous work, Au NPs heated to 900 ∘C were found to evolve
by a competition between evaporation out of the system and adatom diffusion
along the substrate5. While high temperatures are required to increase the
vapor pressure of solid gold such that evaporation is favorable, we expect
that diffusion of adatoms along the particle surface is still possible at
lower temperatures, and that this diffusion could lead to important
morphological changes. Several recent studies have focused on using atomic-
resolution scanning transmission electron microscopy (STEM) to quantify the
morphology of faceted nanoparticles. Using these methods, researchers have
observed beam-induced adatom diffusion on faceted surfaces, motion of entire
layers leading to stochastic rearrangement of particle morphology, and the
influence of edge atoms on the concentration of diffusing atoms6, 7. Others
have used time-resolved STEM to observe facet changes as a nanoparticle is
subjected to oxidizing and reducing environments and showed reaction-induced
surface roughening8. Catalysis studies have also shown enhanced reaction rates
on specific crystal surfaces and surface-dependent coarsening rates in
supported nanoparticles, further demonstrating the importance of understanding
how catalyst performance and degradation rates depend on nanoscale
morphology9, 10, 11.

In this manuscript, the goal is to understand how the dominant evolution
mechanisms in gold nanoparticles change as a function of temperature, and to
elucidate the role of surface structures on evolution kinetics. We use time-
resolved in situ transmission electron microscopy (TEM) to understand the
mechanisms which contribute to nanoparticle degradation as a function of
temperature. Motivated by our experimental findings, which suggest
temperatures above 750 ∘C are required to activate the evaporation process on
Au NP, we performed kinetic Monte Carlo (KMC) simulations to achieve atomic
scale understanding of the mechanisms which lead to particle-scale shape
changes. By investigating how fundamental particle properties, e.g. the
coordination number of atoms on different surfaces, lead to stochastic
evolution processes we are able to clearly observe the growth of new surface
layers, the formation of non-equilibrium facet orientations, and the
relationship between the particle structure and the evaporation rate. Further,
density functional theory (DFT) calculations are applied to study the low
energy reaction pathways which lead to nanoparticle degradation. This
combination of time-resolved experimental data, atomistic kinetic simulations,
and first principles thermodynamic calculations allow us to build a
comprehensive picture of the relationship between particle surface properties,
environmental conditions, and evaporation rates.

##  1 Results and Discussion

###  1.1 Temperature-dependent Evolution Behavior

Our model catalyst consisted of 5 nm Au nanoparticles supported on a silicon
nitride in situ heating chip. The system was heated to high temperature within
the TEM, and images were collected as a function of time (see Methods section
for further details). Nanoparticle sizes were measured from TEM images, and
particle size as a function of time was tracked for each particle.
Experimental data show that below 750 ∘C particle size changes very little
throughout the experiment. At 600 – 700 ∘C a few nanoparticles evaporate and
there are some examples of particle coalescence (seen as discontinuous
growth/disappearance), however, most particles do not change size
significantly (Figure 1). At 750 ∘C many of the particles in view begin to
rapidly evaporate, and the process accelerates as the temperature continues to
increase. At temperatures of 750 ∘C and above there is no clear evidence for
particle coalescence, and we find that evaporation dominates.

Figure 1: Size (R) vs. time profiles for all particles tracked in experimental
TEM image series. Each plot shows the real time evolution of nanoparticles as
a function of time and temperature.

By directly observing changes to the shapes and sizes of individual
nanoparticles, we can identify patterns in the evolution behaviors as
nanoparticles shrink. Figure 2 shows the evolution of an evaporating
nanoparticle which clearly suggests preferential evaporation from specific
crystallographic facets, as throughout the evaporation processes flat facets
and sharp angular notches form at the particle surface.

Figure 2: Evolution of an evaporating particle (circled in red at t = 0 s)
through time. Faceted notches form at the surface of the particle and trigger
fast evaporation. Scale bar represents 5 nm.

To better understand the kinetics and temperature-dependence of the system,
the evaporation rate was plotted as a function of temperature. Assuming that
evaporation is the only process which removes Au mass from the model catalyst,
the total volume of gold in the system (assuming spherical particles) was
calculated in each frame of the TEM dataset, fit to a line as a function of
time, and the evaporation rate was taken as the slope of the line (Figure 3
A). For all temperatures, a linear relationship between total Au volume as a
function of time showed clear agreement with experimental data. A dramatic
increase in evaporation rate is observed as the temperature approaches and
exceeds 750 ∘C (Figure 3A). We find an agreement with the familiar Arrhenius
form (Figure 3B) suggesting that the evaporation of Au from nanoparticle
surfaces is kinetically activated and not a dominant reaction at all
temperatures. The extracted activation energy for evaporation of Au was found
to be 1.2 eV. We performed DFT calculations showing that this activation
energy value is on the same order as, but still significantly less than, the
calculated enthalpy of sublimation from any low energy Au surface (see table
of desorption energies in Supplemental Table 1 and Supplemental Figure 1).
This combination of experimental observations and first-principles
calculations demonstrate that a more complicated mechanism than simple
desorption is responsible for nanoparticle evolution.

Figure 3: A.) The evaporation rate is determined by fitting the slope of the
the volume of Au in the system as a function of time. Very little mass is lost
below 700 ∘C. B.) The evaporation rates fit and Arrhenius form, such that an
activation barrier for evaporation can be extracted.

###  1.2 Atomic-scale Morphological Changes

Motivated by our experimental observations of evaporation as an activated
process and clear transitions between stable surface structures in coalescing
and evaporating nanoparticles, we carried out kinetic Monte Carlo (KMC)
simulations to stochastically model the atomic processes that lead to
morphological changes in the nanoparticles. Our simulation starts with model
particle corresponding to an equilibrium-shape Au nanoparticle with atoms
fixed on a perfect face-centered cubic (FCC) lattice and considers diffusion
of surface atoms to unoccupied neighboring lattice sites (diffusion into
vacuum, or lattice sites with coordination number zero, represent
evaporation). The nanoparticle is modeled as an without considering the
support even though the experimental system consists of supported particles.
Our previous research demonstrated that Au nanoparticle evolution (on an
identical support) is dominated by evaporation from free surfaces with only a
minimal contribution from diffusion along the support surface, and, due to our
particle deposition method, the particle-support contact area is very small5.
Directly modeling metal-support interactions would be more important in
systems where NP wet the support - here, the morphology of the entire particle
could change as the contact area changes. Activation energies are based on a
simple broken bond model described in the Methods section. Figure 4 shows the
model particle at several points in the KMC simulation. Using three
dimensional animations to track step-by-step changes in the structure, we
identified three unique processes by which the particle morphology changes:
the formation of adatom clusters on large facets (Figure 4A), the appearance
of facets not present in the initial structure (Figure 4B), and cluster
pinning at defect sites which leads to the growth of new layers (Figure 4C).
Due to the relatively large number of broken bonds required for evaporation
compared to diffusion and the exponential decrease in event probability as a
function of coordination number, most KMC events at all temperatures
correspond to adatom diffusion on existing crystal surfaces. We find that when
three or more adatoms meet on a surface they tend to form a mobile cluster
that migrates on large (1 1 1) surfaces. In some cases, these clusters are
transient and only exist for short time periods, but in others they stay
together while their configuration rapidly fluctuates. As more adatoms join
the cluster, it becomes more stable in shape and less mobile on the surface.
This is because cluster motion would require concerted movements of many atoms
which is unlikely in our simple model. However, we do see that small clusters
become pinned when they encounter facet edges or surface vacancies. The
pinning of these clusters at high energy defect sites allows them to remain
stationary and continue to grow into a low energy configuration through the
addition of more adatoms, enabling the growth of new layers/terraces on top of
the initial structure. In agreement with the well-known Burton-Cabrera-Frank
model for step-flow kinetics, we see that under-coordinated edge and corner
atoms are removed first, on average, which triggers removal of the remaining
edge atoms12, 13. The collective removal of the edges between two facets forms
new facets with different crystallographic orientation. For example, removal
of the edge between {1 1 1} and {1 0 0} facets creates a {3 1 1} plane, while
the removal of an edge between two {1 1 1} facets creates a {1 1 0} plane.

Figure 4: KMC results show the evolution of an evaporating nanoparticle. In
A.), a mobile cluster of adatoms is shown on a {1 1 1} surface. B.) shows the
development of a {1 1 0} facet through the removal of edge atoms. C.) shows
the growth of a new layer on top of the original particle surface.

In order to more clearly understand the changes induced by the formation of
new layers and facets, and to understand how this relates to the behavior of
evaporating nanoparticles we carried out diffraction simulations using the
XrDebye module of the ASE python library. The Debye equation calculates
diffraction patterns based on the wave interference caused by all unique pairs
of atoms in the system; summing over all interference patterns builds a
diffraction pattern which agrees with experimental results and calculations
based on other methods, and does not require a mathematical description for
shape and structure factors14. Simulated diffraction patterns are displayed as
a function of time in Supplemental Figure 2A. Using the assumption that, for a
monometallic particle, changes in diffraction peak intensities corresponds to
changes in the volume fraction of diffracting surfaces, peak intensity ratios
with respect to the dominant {1 1 1} diffraction intensity as a function of
time were used to qualitatively describe morphological changes in the
nanoparticle during evaporation15. Our results show increases in the relative
intensities of both (2 0 0) and (3 1 1) peaks. While the (2 0 0) intensity
increases linearly through the course of the simulation, the (3 1 1) peak
increases linearly with a shallow slope initially, before rapid increase at
later times. The rapid increase in (3 1 1) intensity coincides very closely
with the time step at which the particle evaporation rate begins to rapidly
increase (Supplemental Figure 2B and C) suggesting that the formation, or
exposure, of {3 1 1} is correlated with an increased evaporation rate.

The {$n$ 1 1} surfaces (for odd $n$) are vicinal surfaces for the {1 0 0}
family. This means that they are slightly misoriented from the low index {1 0
0} planes, and that increasing $n$ corresponds to increasing the width of {1 0
0} terraces. The surface science literature shows both experimental and
theoretical verification for the formation of specific, stable vicinal
surfaces during the roughening process to help stabilize high energy or highly
misoriented surfaces. For example, Bartolini demonstrated that {3 1 1}, {5 1
1} and {11 1 1} are “magic” vicinals for FCC crystals which are exceptionally
stable in evolving surfaces16. It is suggested that these orientations are
most stable because they are able to accommodate small local distortions which
produce a “corrugated surface”, rather than a surface with many high energy
step edges. Our KMC model shows two mechanisms by which {$n$ 1 1} planes can
develop during the diffusion/evaporation evolution process. First, the removal
of atoms at the edges between {1 0 0} and {1 1 1} faces forms {3 1 1} planes
directly (Figure 5). Further removal of edge atoms would increase $n$ and
increase the width of {1 0 0} terraces. Secondly, the pinning of mobile adatom
clusters at edge defects on the particle surface, and subsequent growth to
form a new layer, can form this type of surface. In agreement with this
mechanistic picture of {3 1 1} formation via fast evaporation from the removal
of edges, formation of steps, and subsequent layer evaporation, the angle
between the evaporating surfaces in our experimental images of evaporation
corresponds very closely to the angle between {1 0 0} and{1 1 1} faces (Figure
2, t = 180 s), where edge removal would form {$n$ 1 1} vicinal surfaces.
Moreover, the side length of the notch which initially forms corresponds with
the lattice parameter of the Au FCC structure suggesting a step-wise evolution
pattern.

Figure 5: Schematic view of how {$n$ 1 1} facets can be formed either by
removing atoms from the particle edge, or the growth of new layers. The purple
atoms represent the atoms which are added/removed to make the purple {3 1 1}
plane in the right-most structure.

###  1.3 Desorption Mechanism and Kinetic Equilibrium

DFT calculations were utilized to provide an atomistic explanation for
nanoparticle evaporation via preferential desorption from specific crystal
facets. The model consisted of calculated surface defect formation energies
($\Delta G_{defect}$) and gold adatom adsorption energies ($\Delta
G_{adsorption}$). We considered low Miller index (1 0 0), (1 1 0), (1 1 1), (2
1 1) and (3 1 1) facets which would be expected to dominate nanoparticles in
the size and temperature ranges considered experimentally. Since initial
surface defect formation energies and adsorption energies of adatoms on the (1
1 1) surface sites were significantly higher than the experimental activation
barrier (Supplemental Table 1), we hypothesized that the evaporation process
under these conditions is driven by an adatom formation step followed by a
rate-limiting adatom desorption step and that these processes are entropically
driven at elevated temperatures. Guided by the experimental system design in
which nanoparticles are under continuously pumped at high vacuum conditions
leading to non-conservative conditions, the DFT model suggests an equilibrium
between the existence of surface facets and mobile (1 1 1) adatoms (Figure 6)
and further assumes that evaporated gold atoms cannot redeposit on the
surfaces of the nanoparticles.

Figure 6: Schematic illustrating the evaporation mechanism whereby surface
facets are in equilibrium with (111) adatoms that undergo rate-limiting
desorption into the vacuum where they are swept away before equilibrium can be
established. Figure 7: Free energy changes and equilibrium constants
associated with the formation of (111) adatoms from the surface atoms of the
(111), (110), (100), (211), and (311) surface facets. At the top is also the
calculated free energy change for the desorption of a (111) adatom into the
gas phase at 750 C.

Hence, our model of the evaporation process is as follows: upon heating, the
nanoparticle reaches a point at which the equilibrium constant for the adatom
forming reaction step becomes sizeable enough for substantial adatom formation
and corresponding surface atom removal. These adatoms form from the low-
coordination surfaces (steps, kinks, etc.) as evidenced by the relatively low
formation energies for the (3 1 1) and (1 1 0) surfaces shown in Figure 7.
This process is ultimately entropically driven; the adatoms and gas phase gold
atoms have significant translational entropy contributions and at elevated
temperatures they become more stable relative to the higher index surface
facets. Formation energies in Figure 7 are calculated using $\Delta
G_{750}=\Delta G_{defect}+\Delta G_{adsorption}$, where $\Delta G_{750}$ is
the energy required to take a gold atom from its lattice position in the
surface and place it as an adatom adsorbed on the (1 1 1) surface at 750 ∘C.
$\Delta G_{defect}$ is the free energy change associated with creating a
defect in the surface by removing an atom to the gas phase, and $\Delta
G_{adsorption}$ is the adsorption energy of a gold adatom on the (1 1 1)
surface relative to the gas phase. Equilibrium constants are calculated as
$K_{750}=e^{\frac{-\Delta G_{750}}{k_{B}T}}$. Calculated adsorption energies
for gold adatoms on the three-fold and bridge sites of the (1 1 1) surface
confirm that the activation barrier for the surface diffusion of these adatoms
is very small (0.01 eV), implying that they have two essentially unhindered
translational degrees of freedom and substantial entropy at the temperatures
relevant to the experiments (Supplemental Table 2). The adatoms then can
undergo desorption with $\Delta G\approx 1$ eV at 750 ∘C, in direct agreement
with the experimentally determined activation barrier, supporting this two-
step picture of adatom formation and desorption from {1 1 1} terraces (Figure
7). To confirm the DFT findings, we calculated the equilibrium concentration
of {1 1 1} adatoms on model nanoparticle surfaces directly from the KMC
simulations described above. Here, {1 1 1} adatoms are defined as atoms in a
KMC step with a coordination number of three. Supplemental Figure 3 shows a
plot of the adatom concentration (calculated as the number of {1 1 1} adatoms
divided by the total number of surface atoms), where we can clearly see that
the equilibrium KMC adatom concentration at all temperatures converges to
1$\%$ in agreement with equilibrium constant $K_{eq,750}=0.012$.

For temperatures of 750 ∘C and below, where the evaporation rate is low and
surface roughening occurs slowly, we see experimental evidence for
preferential evaporation from specific crystallographic surfaces. Our
experimental data show the evolution of an evaporating particle (Figure 2).
The evaporation process begins by the formation of a notch in the edge of the
nanoparticle at a nearly 90∘ angle. The notch continues to grow and forms a
step edge while the surrounding particle surface becomes smoother.
Subsequently, the step edge rapidly disappears, and we see the formation of a
flat particle surface. The process continues in the same way until evaporation
is complete. It is clear that certain crystallographic faces of the Au
nanoparticles evaporate or evolve faster than others to cause these
discontinuous shrinkage patterns. Our KMC results give us an atomistic picture
of how these facet relationships form and evolve using the varied coordination
number of surface atoms as a proxy for anisotropic surface energy. The
simulations, in agreement with the mechanism proposed above, show that the
desorption process in which a surface atom is removed to form an adatom on the
(111) terrace (leading to evaporation or surface diffusion) is most likely
initiated from the edges and corners between facets since these sites have the
lowest coordination number. After removal from the edge site, adatoms are most
likely to diffuse along the surface rather than directly evaporate due to the
high activation barrier for evaporation; this leads to the presence of mobile
atoms on the crystal surfaces in equilibrium with the facets themselves. These
surface adatoms are even more under-coordinated after leaving the edge sites
and migrate in a random walk until they either find a low energy binding site
(near a surface defect, forming a cluster with other adatoms, or at step
edges) or evaporate17.

###  1.4 Nanoparticle Surface Roughening

Based on our results, which show experimentally that evaporation begins around
750 ∘C and, computationally that the evaporation process depends on a
morphological departure from a smooth surface structure, we expect that the
ability of Au atoms to evaporate is closely related to the roughening
transition. While the actual appearance of a discontinuity in surface
roughness to signify the roughening transition is debated in the literature,
our combination of experimental and simulated data show the presence of rough
crystal surfaces and rapid fluctuations in surface structure that cause
roughening18. Moreover, our measurements of evaporation rate as a function of
temperature from experimental data appear nearly discontinuous, with a clear
onset of evaporation near 750 ∘C. DFT confirms that the energetic barrier for
adatom creation is low for the high-index surfaces which form through the
evolution process. This reinforces our understanding that rough surfaces,
which form in part through the creation of {1 1 1} adatoms from surface edges,
are critical to the evaporation process throughout the temperature range we
considered. At higher temperatures the equilibrium will shift such that
adatoms on the (111) surface are actually more stable than the low-coordinated
surfaces themselves, eventually enabling rapid morphology change and
evaporation. With this in mind, though our KMC models assumed a perfect FCC
lattice and initial morpholgy determined byt the equiilibrium Wulff
construction, we suspect that the evolution behavior is not strongly dependent
on the initial configuration of the nanoparticle. In all cases, we observe
surface roughening and that evolution is driven by evaporation of adatoms on
rough surface, not by desportion from specific facets, edges, or corners.
Researchers have previously used ab initio modeling to map the morphological
phase diagram of Au nanoparticles as a function of size19. They showed, that
for particle sizes seen in our experiments, then onset of roughening occurs
around 600 ∘C and that by 850 ∘C the nanoparticle should be “quasi-molten”.
Based on this, our experiments essentially probe the entire temperature range
over which rough surfaces are expected for small gold nanoparticles. Our
experimental data do not show evidence for immediate and spontaneous
roughening or surface melting at any temperature, however, we do see rapid
changes in the surface structure which signal the progression between the
perfect faceted surface and the high-temperature rough configuration.

Roughening also decreases the mean coordination number of surface atoms
compared to the clean crystal surface. Considering the physics behind our KMC
model, where the probability of desorption/evaporation is inversely related to
the coordination number, we can see how the roughest surfaces at the highest
temperatures are exponentially more prone to evaporation due to thermal
activation and large concentrations of under-coordinated adatoms. Previous KMC
studies have shown the link between evaporation and surface roughening, and
described the acceleration of roughening and evaporation in systems where mass
is not conserved18. In a different system, it has been shown that the presence
of O2 gas reduces the difference in surface energy between varied Pt facets,
and leads to the formation of round particles with rough surfaces8. Similarly,
research has suggested that the desorption and surface diffusion PtOx species
from rough nanoparticle surfaces are responsible for rapid coarsening via
Ostwald Ripening for alumina-supported Pt nanoparticles20. As gold is more
noble than Pt, even at the nanoscale, the rough nanoparticle surface is
subjected to a driving force for evaporation which comes from the pressure
difference between the ambient vacuum and the high Gibbs-Thomson pressure at
the particle surface, perhaps explaining why we see evaporation rather than
any particle growth in our experiments.

##  2 Conclusions

Our use of a combination of experimental and computational techniques has
allowed us to understand the mechanisms and rates of supported nanoparticle
evaporation as a function of temperature and particle surface characteristics.
Using time-resolved in situ TEM, we were able to observe the real behavior of
systems containing hundreds of supported gold nanoparticles at various
processing temperatures, and track their evolution in time. Motivated by
experimental results which signaled that an activated mechanism drives
nanoparticle evaporation, KMC simulations allowed us to suggest atomistic
pathways that enable evaporation alongside morphological changes observed
experimentally. Finally, DFT calculations were used to reinforce the link
between surface morphology and propensity for degradation. These calculations
show that evaporation occurs through a two-step mechanism where, first, high-
energy edge atoms form adatoms on larger {1 1 1} terraces, and subsequently
the high temperature in the system enables evaporation with an activation
barrier confirmed in our experimental analysis. While our experimental and
computational results focus specifically on the evolution of Au nanoparticles,
our approach - utilizing common experimental techniques and simulations which
include minimal material-specific assumptions - is generalizable to other
systems.

With our model of the process by which high system temperatures induce surface
roughening and subsequent evaporation, and the knowledge that this process
initiates at the edges of high-index surfaces, we can move towards
understanding subtleties of supported catalyst degradation without relying on
simplifications and assumptions about nanoparticle shape and surroundings. In
this way, findings like ours, which suggest important parameters for catalytic
particles, can inform research in materials synthesis to create novel catalyst
materials which are less susceptible to degradation under harsh reaction
conditions.

##  3 Methods

###  3.1 In situ TEM Experiments

Au nanoparticles coated with dendrimer ligands to enable self-assembly at
specified particle spacings were synthesized according to previous study21.
The nanoparticles used in this experiment had a monodisperse particle size
distribution centered at R = $2.1\pm 0.3$ nm. After synthesis, nanoparticles
were drop-cast directly on to Hummingbird Scientific in situ heating chips22.
After drop-casting, heating chips were plasma cleaned under O2 (forward target
power = 25 W) for 15 seconds to remove the ligand coating before heating.
Finally, the chip was loaded into a Hummingbird Scientific heating holder and
inserted directly into a JEOL F200 S/TEM operating with accelerating voltage
of 200 kV. Chip windows were scanned to find regions containing large areas of
self-assembled nanoparticles, and the samples were subsequently heated to
temperature. Images were collected every five seconds for at least one hour
using a Gatan OneView electron detector with 4k pixel resolution. Experimental
temperatures range from 650 – 850 ∘C in 50 ∘C increments.

###  3.2 Image Segmentation and Analysis

After each experiment, an unsupervised CNN pipeline was used to segment images
and extract particle sizes and positions as a function of time; background
information, machine learening Python code, and an interactive image
segmentation tutorial are provided in the references23, 24. A small set of
experimental images ( 5 images) were randomly chosen from each dataset to
build a training set for an unsupervised autoencoder model and each 4k image
was cropped into non-overlapping 128 x 128-pixel sections. Each subsection was
then rotated by 90°, 180°, and 270° to augment the dataset . Since CNN are
translation-invariant and recognize image features regardless of their
position, each of these image subsections and rotations represents a unique
training example. The unique UNet-type unsupervised CNN was retrained for each
experimental dataset, and then applied to raw images using the PyTorch
framework25. Particle sizes and positions were extracted from binary images
using the scikit-image connected components algorithm26. After processing all
experimental datasets, a comprehensive Pandas DataFrame was used to store all
experimental data as a function of particle size, time, and temperature. All
further analysis and visualization were performed in Python, and made use of
numpy, matplotlib, and scipy libraries27, 28, 29. Jupyter notebooks
demonstrating the analysis and production of figures are available on github
at: https://github.com/jhorwath/FacetEvaporation.git.

###  3.3 Kinetic Monte Carlo Simulations

A homemade script was used for kinetic Monte Carlo (KMC) simulations
(simulation scripts are available on the github repository listed above). To
begin, a NP model was made by sampling the experimental particle size
distribution, and forming the equilibrium FCC structure for the particle based
on the Wulff construction applied in the Atomic Simulation Environment (ASE)
python package to use as simulation input30. KMC event probabilities were
calculated according to Equation 1:

$P_{event}=v_{0}\exp\left[\frac{-\epsilon N}{k_{B}T}\right]$ (1)

where $v_{0}$ is an attempt frequency chosen to match the experimental
physics, $\epsilon$ is an energy per bond (taken as 0.1 eV / bond) multiplied
by coordination number $N$31. $k_{B}$ and $T$ are Boltzmann’s constant and
temperature, respectively.

In each KMC step, the event probability was calculated for each atom in the
particle to determine which atoms will move/evaporate32, 33, 34, 35. In order
to accelerate the simulation, only surface atoms ($N<12$) were considered in
KMC time steps. For each event, all available lattice sites surrounding an
atom were considered as possible diffusion sites, and evaporation was
considered as diffusion into the vacuum. To determine which move to carry out,
a rate catalog based on the activation barrier for each candidate was built,
and then sampled using a random number generator. The diffusion event rate,
$r$, was determined by Equation 2:

$r=v_{0}\exp\left[\frac{-\Delta E}{k_{B}T}\right]$ (2)

where $\Delta E_{diffusion}=-\epsilon(N_{1}-N_{0})$ for coordination numbers
$N_{1}$ of the candidate site and $N_{0}$ at the current site.
Correspondingly, $\Delta E_{evaporation}=\epsilon N_{0}$ since evaporation
requires breaking all surface bonds. By this model, diffusion to a more
coordinated lattice site is favorable, and evaporation occurs with low
probability but is still possible for atoms with low coordination number at
high temperatures. Simulations were run until all atoms in the system
evaporated. To ensure that KMC results matched the expected physical behavior,
KMC evaporation rates were plotted to show that
$\frac{dR}{dt}\propto\frac{1}{R}$, in agreement with our previous research on
nanoparticle evaporation (Supplemental Figure 4)5. 3D visualizations of KMC
trajectories were made using the Plotly python library36. Diffraction
simulations of KMC nanoparticles were carried out using the Debye scattering
equation and were performed using the XrDebye code available in the ASE python
package.

###  3.4 Density Functional Theory

For the DFT work, all systems are built using the Atomic Simulation
Environment (ASE) 37 and all calculations are performed using Quantum
Espresso’s 38 implementation of generalized gradient approximation (GGA) DFT
using the revised Perdew-Burke-Ernzerhof (RPBE) exchange-correlation
functional 39 and a plane wave cutoff of 400 eV. A 5x5x1 Monkhorst-Pack grid
40 is used for k-point sampling on all relaxations. Periodic boundary
conditions are implemented with a dipole correction 41. All calculations are
converged to 10-5 eV and forces are converged to 0.03 eV/Å. Core electrons are
treated using ultrasoft pseudopotentials 42. Fermi-Dirac smearing is used with
a smearing width of 0.1 eV. Calculations are not spin-polarized.

Entropy corrections are calculated using the ASE Thermochemistry package 30.
Adatoms on the (111) surface are treated using the Hindered Translator method
43 with a diffusion barrier of 0.01 eV (equal to the difference in adsorption
energies for a gold adatom on the threefold and bridge sites of the (111)
surface) and only translational entropy is included. Surface lattice atoms are
treated as harmonic oscillators with only vibrational entropy contributions.
Gas phase atoms are treated as ideal gas with an entropy of 0.002 eV/K.

DFT calculated ground state energies are then corrected with these entropies
to get the free energy of the system according to Equation 3:

$G=E_{DFT}-TS$ (3)

##  4 Acknowledgments

J.P.H. and E.A.S acknowledge support through the National Science Foundation,
Division of Materials Research, Metals and Metallic Nanostructures Program
under Grant 1809398. We thank Katherine Elbert, Shengsong Yang, and
Christopher Murray from the University of Pennsylvania Department of Chemistry
for help with nanoparticle and dendrimer ligand synthesis. Additionally, we
acknowledge the support of Peter Voorhees, of the Northwestern University
Department of Materials Science and Engineering, for guidance while developing
KMC models and discussion of surface evolution kinetics. C.L.C. would like to
thank the Vagelos Institute for Energy Science and Technology (University of
Pennsylvania) for support through a graduate fellowship. C.L.C. and A.V.
acknowledge the use of the computer time allocation at the Extreme Science and
Engineering Discovery Environment (XSEDE) supported through National Science
Foundation Energy under Award Number DMR180108.

## References

  * Liu and Corma 2018 Liu, L.; Corma, A. Metal Catalysts for Heterogeneous Catalysis: From Single Atoms to Nanoclusters and Nanoparticles. _Chem. Rev._ 2018, 99
  * Balluffi et al. 2005 Balluffi, R. W.; Allen, S. M.; Carter, W. C.; Kemper, R. A. _Kinetics of materials_ ; J. Wiley & Sons: Hoboken, N.J, 2005
  * Voorhees 1985 Voorhees, P. W. The theory of Ostwald ripening. _Journal of Statistical Physics_ 1985, _38_ , 231–252
  * Williams and Bartelt 1989 Williams, E. D.; Bartelt, N. C. Surface faceting and the equilibrium crystal shape. _Ultramicroscopy_ 1989, _31_ , 36–48
  * Horwath et al. 2021 Horwath, J. P.; Voorhees, P. W.; Stach, E. A. Quantifying Competitive Degradation Processes in Supported Nanocatalyst Systems. _Nano Letters_ 2021, _21_ , 5324–5329
  * De wael et al. 2020 De wael, A.; De Backer, A.; Jones, L.; Varambhia, A.; Nellist, P. D.; Van Aert, S. Measuring Dynamic Structural Changes of Nanoparticles at the Atomic Scale Using Scanning Transmission Electron Microscopy. _Physical Review Letters_ 2020, _124_ , 106105
  * Liu et al. 2021 Liu, P.; Arslan Irmak, E.; De Backer, A.; De wael, A.; Lobato, I.; Béché, A.; Van Aert, S.; Bals, S. Three-dimensional atomic structure of supported Au nanoparticles at high temperature. _Nanoscale_ 2021, _13_ , 1770–1776
  * Altantzis et al. 2019 Altantzis, T.; Lobato, I.; De Backer, A.; Béché, A.; Zhang, Y.; Basak, S.; Porcu, M.; Xu, Q.; Sánchez-Iglesias, A.; Liz-Marzán, L. M.; Van Tendeloo, G.; Van Aert, S.; Bals, S. Three-Dimensional Quantification of the Facet Evolution of Pt Nanoparticles in a Variable Gaseous Environment. _Nano Letters_ 2019, _19_ , 477–481
  * Yuan et al. 2018 Yuan, W.; Zhang, D.; Ou, Y.; Fang, K.; Zhu, B.; Yang, H.; Hansen, T. W.; Wagner, J. B.; Zhang, Z.; Gao, Y.; Wang, Y. Direct In Situ TEM Visualization and Insight into the Facet‐Dependent Sintering Behaviors of Gold on TiO ${}_{\textrm{2}}$. _Angewandte Chemie International Edition_ 2018, _57_ , 16827–16831
  * Bi et al. 2011 Bi, Y.; Ouyang, S.; Umezawa, N.; Cao, J.; Ye, J. Facet Effect of Single-Crystalline Ag ${}_{\textrm{3}}$ PO ${}_{\textrm{4}}$ Sub-microcrystals on Photocatalytic Properties. _Journal of the American Chemical Society_ 2011, _133_ , 6490–6492
  * Chmielewski et al. 2019 Chmielewski, A.; Meng, J.; Zhu, B.; Gao, Y.; Guesmi, H.; Prunier, H.; Alloyeau, D.; Wang, G.; Louis, C.; Delannoy, L.; Afanasiev, P.; Ricolleau, C.; Nelayah, J. Reshaping Dynamics of Gold Nanoparticles under H ${}_{\textrm{2}}$ and O ${}_{\textrm{2}}$ at Atmospheric Pressure. _ACS Nano_ 2019, acsnano.8b08530
  * Burton et al. 1951 Burton, W. K.; Cabrera, N.; Frank, F. C.; Mott, N. F. The growth of crystals and the equilibrium structure of their surfaces. _Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences_ 1951, _243_ , 299–358, Publisher: Royal Society
  * Burton and Cabrera 1949 Burton, W. K.; Cabrera, N. Crystal growth and surface structure. Part I. _Discussions of the Faraday Society_ 1949, _5_ , 33
  * Derlet et al. 2005 Derlet, P. M.; Van Petegem, S.; Van Swygenhoven, H. Calculation of x-ray spectra for nanocrystalline materials. _Physical Review B_ 2005, _71_ , 024114
  * Zhang et al. 2020 Zhang, L.; Gonçalves, A. A. S.; Jaroniec, M. Identification of preferentially exposed crystal facets by X-ray diffraction. _RSC Advances_ 2020, _10_ , 5585–5589
  * Bartolini et al. 1989 Bartolini, A.; Ercolessi, F.; Tosatti, E. “Magic” vicinal surfaces stabilized by reconstruction. _Physical Review Letters_ 1989, _63_ , 872–875
  * Pimpinelli and Villain 1994 Pimpinelli, A.; Villain, J. What does an evaporating surface look like? _Physica A: Statistical Mechanics and its Applications_ 1994, _204_ , 521–542
  * Zhdanov and Kasemo 1998 Zhdanov, V. P.; Kasemo, B. Surface roughening: Kinetics, adsorbate-induced effects, and manifestation in catalytic reactions. _The Journal of Chemical Physics_ 1998, _108_ , 4582–4590
  * Barnard et al. 2009 Barnard, A. S.; Young, N. P.; Kirkland, A. I. Nanogold: A Quantitative Phase Map. 2009, _3_ , 6
  * Simonsen et al. 2010 Simonsen, S. B.; Chorkendorff, I.; Dahl, S.; Skoglundh, M.; Sehested, J.; Helveg, S. Direct Observations of Oxygen-induced Platinum Nanoparticle Ripening Studied by In Situ TEM. _Journal of the American Chemical Society_ 2010, _132_ , 7968–7975
  * Elbert et al. 2017 Elbert, K. C.; Jishkariani, D.; Wu, Y.; Lee, J. D.; Donnio, B.; Murray, C. B. Design, Self-Assembly, and Switchable Wettability in Hydrophobic, Hydrophilic, and Janus Dendritic Ligand–Gold Nanoparticle Hybrid Materials. _Chem. Mater._ 2017, 10
  * 22 MEMS Heating and Biasing. https://hummingbirdscientific.com/products/heating-biasing/
  * Horwath et al. 2020 Horwath, J. P.; Zakharov, D. N.; Mégret, R.; Stach, E. A. Understanding important features of deep learning models for segmentation of high-resolution transmission electron microscopy images. _npj Computational Materials_ 2020, _6_ , 108
  * 24 Vyas, L.; Horwath, J. P.; Stach, E. A. Tutorial on Unsupervised Image Segmentation for Electron Microscopy. https://mlforem.github.io/
  * 25 Paszke, A.; Gross, S.; Chintala, S.; Chanan, G.; Yang, E.; DeVito, Z.; Lin, Z.; Desmaison, A.; Antiga, L.; Lerer, A. Automatic differentiation in PyTorch. 4
  * van der Walt et al. 2014 van der Walt, S.; Schönberger, J. L.; Nunez-Iglesias, J.; Boulogne, F.; Warner, J. D.; Yager, N.; Gouillart, E.; Yu, T. scikit-image: image processing in Python. _PeerJ_ 2014, _2_ , e453
  * Harris et al. 2020 Harris, C. R.; Millman, K. J.; van der Walt, S. J.; Gommers, R.; Virtanen, P.; Cournapeau, D.; Wieser, E.; Taylor, J.; Berg, S.; Smith, N. J.; Kern, R.; Picus, M.; Hoyer, S.; van Kerkwijk, M. H.; Brett, M.; Haldane, A.; del Río, J. F.; Wiebe, M.; Peterson, P.; Gérard-Marchant, P. et al. Array programming with NumPy. _Nature_ 2020, _585_ , 357–362
  * Hunter 2007 Hunter, J. D. Matplotlib: A 2D Graphics Environment. _Computing in Science & Engineering_ 2007, _9_ , 90–95
  * Virtanen et al. 2020 Virtanen, P.; Gommers, R.; Oliphant, T. E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; van der Walt, S. J.; Brett, M.; Wilson, J.; Millman, K. J.; Mayorov, N.; Nelson, A. R. J.; Jones, E.; Kern, R.; Larson, E.; Carey, C. J. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. _Nature Methods_ 2020, _17_ , 261–272
  * Larsen et al. 2017 Larsen, A.; Mortensen, J.; Blomqvist, J.; Jacobsen, K. K. W.; Hjorth Larsen, A.; Jørgen Mortensen, J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Bjerre Jensen, P.; Kermode, J.; Kitchin, J. R.; Leonhard Kolsbjerg, E. et al. The atomic simulation environment—a Python library for working with atoms. _Journal of Physics: Condensed Matter_ 2017, _29_ , 273002
  * Pyykko 2004 Pyykko, P. Theoretical Chemistry of Gold. _Angewandte Chemie International Edition_ 2004, _43_ , 4412–4456
  * He et al. 2016 He, X.; Cheng, F.; Chen, Z.-X. The Lattice Kinetic Monte Carlo Simulation of Atomic Diffusion and Structural Transition for Gold. _Scientific Reports_ 2016, _6_ , 33128
  * Voter 2007 Voter, A. F. In _Radiation Effects in Solids_ ; Sickafus, K. E., Kotomin, E. A., Uberuaga, B. P., Eds.; Springer Netherlands: Dordrecht, 2007; Vol. 235; pp 1–23, Series Title: NATO Science Series
  * Jansen 2012 Jansen, A. P. J. _An introduction to kinetic Monte Carlo simulations of surface reactions_ ; Lecture notes in physics 856; Springer: Heidelberg, 2012; OCLC: 797324973
  * 35 Landau, D. P.; Binder, K. A Guide to Monte Carlo Simulations in Statistical Physics, Third Edition. 489
  * Inc 2015 Inc, P. T. Collaborative data science. 2015; https://plot.ly, Place: Montreal, QC Publisher: Plotly Technologies Inc.
  * Larsen et al. 2017 Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Du\lak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Jensen, P. B.; Kermode, J.; Kitchin, J. R.; Kolsbjerg, E. L.; Kubal, J.; Kaasbjerg, K.; Lysgaard, S.; Maronsson, J. B. et al. The atomic simulation environment—a Python library for working with atoms. _Journal of Physics: Condensed Matter_ 2017, _29_ , 273002, Publisher: IOP Publishing
  * Giannozzi et al. 2009 Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Corso, A. D.; Gironcoli, S. d.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. _Journal of Physics: Condensed Matter_ 2009, _21_ , 395502, Publisher: IOP Publishing
  * Hammer et al. 1999 Hammer, B.; Hansen, L. B.; Nørskov, J. K. Improved adsorption energetics within density-functional theory using revised Perdew-Burke-Ernzerhof functionals. _Physical Review B - Condensed Matter and Materials Physics_ 1999, _59_ , 7413–7421
  * Pack and Monkhorst 1977 Pack, J. D.; Monkhorst, H. J. ”Special points for Brillouin-zone integrations”—a reply. _Physical Review B_ 1977, _16_ , 1748–1749
  * Bengtsson 1999 Bengtsson, L. Dipole correction for surface supercell calculations. _Physical Review B - Condensed Matter and Materials Physics_ 1999, _59_ , 12301–12304
  * Vanderbilt 1990 Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. _Physical Review B_ 1990, _41_ , 7892–7895
  * Sprowl et al. 2016 Sprowl, L. H.; Campbell, C. T.; Árnadóttir, L. Hindered translator and hindered rotor models for adsorbates: Partition functions and entropies. _Journal of Physical Chemistry C_ 2016, _120_ , 9719–9731


Supp_GKA_arxiv.pdf ††footnotetext: † These authors contributed equally to this
work. ††footnotetext: ∗ Email<EMAIL_ADDRESS>
# Goodenough-Kanamori-Anderson high-temperature ferromagnetism in tetragonal
transition-metal xenes

U. Yorulmaz† Department of Physics, Eskisehir Osmangazi University, Eskisehir,
Turkiye  Department of Physics & NANOlab Center of Excellence, University of
Antwerp, Groenenborgerlaan 171, B-2020 Antwerp, Belgium  D. Šabani† Department
of Physics & NANOlab Center of Excellence, University of Antwerp,
Groenenborgerlaan 171, B-2020 Antwerp, Belgium  C. Sevik  Department of
Mechanical Engineering, Eskisehir Technical University, Eskisehir, Turkiye
Department of Physics & NANOlab Center of Excellence, University of Antwerp,
Groenenborgerlaan 171, B-2020 Antwerp, Belgium  M. V. Milošević∗ Department of
Physics & NANOlab Center of Excellence, University of Antwerp,
Groenenborgerlaan 171, B-2020 Antwerp, Belgium

###### Abstract

Seminal Goodenough-Kanamori-Anderson (GKA) rules provide the inceptive
understanding of the superexchange interaction of two magnetic metal ions
bridged with an anion, and suggest fostered ferromagnetic interaction for
orthogonal bridging bonds. However, there are no examples of two-dimensional
(2D) materials with structure that optimizes the GKA arguments towards
enhanced ferromagnetism and its critical temperature. Here we reveal that an
ideally planar GKA ferromagnetism is indeed stable in selected tetragonal
transition-metal xenes (tTMXs), with Curie temperature above 300 K found in
CrC and MnC. We provide the general orbitally-resolved analysis of magnetic
interactions that supports the claims and sheds light at the mechanisms
dominating the magnetic exchange process in these structures. With recent
advent of epitaxially-grown tetragonal 2D materials, our findings earmark
tTMXs for facilitated spintronic and magnonic applications, or as a desirable
magnetic constituent of functional 2D heterostructures.

††preprint: APS/123-QED

##  I Introduction

The experimental discovery of the premiere magnetic two-dimensional materials
(M2DMs) - CrI3 [1] and CrGeTe3 [2, 3] \- opened the floodgates to many
emergent 2D materials of this class. Numerous theoretical [4, 5, 6, 7] and
experimental studies [8, 9, 10, 11] followed, explaining the origins and
possible manipulations of the long-range magnetic order in the monolayer
limit. It is needless to emphasize that intrinsically room-temperature M2DMs
would be highly beneficial for applications in sensing, spintronics, and
otherwise, and bear promise towards high tunability by diverse mechanical,
chemical, and electronic means. However, it quickly became clear that critical
temperatures (T${}_{\textrm{c}}$) for magnetic order to vanish are by rule
always smaller in 2D materials compared to their bulk counterparts [1, 12,
13]. Namely, in order to host sizable regions with magnetic order (and
circumvent limitations imposed by the Mermin-Wagner theorem [14]), M2DMs
require anisotropy in magnetic exchange. That needed anisotropy is known to
originate from the spin-orbit coupling (SOC), which is in general much weaker
compared to e.g. Coulomb attraction or repulsion of charged particles, hence
attains very low magnitudes (of typically 0.01-0.1 meV). Another reason for
T${}_{\textrm{c}}$ to decrease as one thins the magnetic material from bulk to
a monolayer stems from the correspondingly diminishing magnetic exchange along
the third, out-of-plane direction. One may therefore expect that magnetic
order in 2D materials is strictly limited to the very low T${}_{\textrm{c}}$,
but that is not necessarily the case. For example, Fe3GeTe2 hosts the
ferromagnetic (FM) order up to 130 K in the monolayer (ML) limit [15]. In
addition, the FM order and T${}_{\textrm{c}}$ of 213 K were measured in few-
layer thick 1T-CrTe2 [16], with the unusual trend that T${}_{\textrm{c}}$
increases as one goes from bulk to few-layers. Furthermore, in a few-layer
FePS3, the antiferromagnetic (AFM) order was observed up to the
T${}_{\textrm{c}}$ of 120 K [17]. Finally, FM order was measured even at room
temperature in thicker films ($\approx 10$ nm, see e.g. [18, 19]), but also in
monolayer MnSex [20, 21], VSe2 [11, 22], and Cr3Te4 [23]. The latter samples
were deposited epitaxially, which in general involves structural defects [22]
and strong interfacial effects with the substrate [22, 23] into the origins of
the observed robust magnetic interactions, which complicates theoretical
interpretations. Otherwise, the isotropic magnetic interactions in crystalline
monolayers are relatively straightforward to extract theoretically in all
available first-principles codes. Such studies, on predominantly in silico
created 2D materials, have yielded many predictions of high- or even room-
temperature intrinsic magnetism [24, 25, 26, 27, 28, 21]. However, those
predictions typically failed to quantify the underlying microscopic mechanisms
for such a large predicted magnetic exchange.

Figure 1: Crystal structure of a monolayer tTMX. TM (cyan) atoms are sketched
with one $d$ orbital, and X (yellow) atoms with one $p$ orbital, since these
orbitals and their hybridization are essential for magnetic interactions in
this system.

In this paper we therefore take a step back, and explore the route towards
room-temperature 2D ferromagnetism starting from the well-established set of
Goodenough-Kanamori-Anderson (GKA) empirical theoretical rules from the late
1950s [29, 30, 31, 32]. We accordingly aim at monolayers with 90∘ between the
nearest magnetic atoms (A) connected by a ligand (X), thus ideally a Lieb-
lattice material of A2X type. However, a magnetic 2D material of such specific
planar structure has not been reported to date, although some Lieb-lattice 2D
materials have been considered computationally for other purposes (see e.g.
Ref. [33]). As a best available choice, for not only geometry but also sizable
SOC, we instead focus on the family of tetragonal transition-metal xenes
(tTMXs; see Fig. 1), seeking a square-lattice planar material among them -
still with 90∘ TM-X-TM nearest-neighbor bonds. Using advanced methodology on
top of the standard first-principles approaches based on Density Functional
Theory (DFT), we computationally validate the structural stability and strong
intrinsic magnetic interactions in these materials, detail the microscopic
(orbital-resolved) origin of enhanced magnetic exchange, and identify CrC and
MnC as premiere square-lattice monolayer ferromagnets with a particularly high
T${}_{\textrm{c}}$.

##  II Results

We commence our analysis with a computationally crude throughput screening of
dynamical stability and magnetic interactions in tTMXs (where TM = V, Cr, Mn,
Fe, Co, Ni, Cu, Zn, and X = C, N, O, S, Se, Te). For each stable material we
perform total energy mapping between the density functional theory corrected
with on-site Coulomb repulsion (DFT+U) and the Heisenberg model Hamiltonian,
for six particular magnetic orders - namely FM and AFM orders along three
Cartesian directions (see Appendix B in † Supplementary Information), in order
to extract the governing magnetic interactions in the system. For the sake of
screening, Hubbard parameter U in the calculations is taken from the online
database, based on bulk oxides of transition metals [34]. In order to decrease
computational cost, we consider only the first nearest-neighbor (NN) magnetic
interactions and the single-ion anisotropy (SIA) in the Heisenberg
Hamiltonian.

The considered primitive unit cell of these materials consists of two TM, and
two X atoms - such that each TM atom has four X atoms as the nearest
neighbors, and vice versa (see Fig. 1), where one can apply symmetry rules for
the exchange matrix (see Appendix A in † Supplementary Information) and SIA
(cf. Ref. [35]). In Fig. 1 each TM atom is sketched with one $d$ orbital, and
each X atom with one $p$ orbital, since these atomic orbitals, and their
mutual hybridization, are essential for the interactions between magnetic
moments on TM atoms. Though the positions of the atoms in the structure are
uniquely determined with respect to the in-plane primitive lattice vectors due
to the symmetry of tetragonal structures, both atomic species are allowed to
relax out-of-plane (along $\vec{a}_{3}\equiv z$ axis).

Table 1: Magnetic properties of the stable monolayer tTMX structures. FM, AFM, and NM stand for ferromagnetic, antiferromagnetic, and non-magnetic order, respectively. Jxx and Jyy mark in-plane, and Jzz out-of-plane exchange interactions. Out-of-plane exchange anisotropy $\Delta$ stands for the difference Jxx-Jzz. A${{}_{ii}^{zz}}$ is the single-ion anisotropy (SIA). If signs of SIA and exchange interactions are the same, SIA favors out-of-plane anisotropy, otherwise the in-plane one. Jxx, Jyy, Jzz, $\Delta$, and A${{}_{ii}^{zz}}$ are all given in meV. T${}_{\textrm{c}}$ stands for the critical temperature of the magnetic order, Curie for FM and Néel temperature for AFM monolayers.  | Magnetic |  |  |  |  |
---|---|---|---|---|---|---
| order |  Jxx=Jyy |  Jzz | $\Delta~{}~{}$ |  A${{}_{ii}^{zz}}$ |  T${}_{\textrm{c}}$(K)
CrC | FM | -52.68 | -52.66 | 0.03 | -0.35 | 515.6
MnC | FM | -105.98 | -106.01 | -0.03 | -0.67 | 1065.8
VN | FM | -3.73 | -10.54 | -6.81 | -13.94 | 152.4
CoN | FM | -12.57 | -12.94 | -0.38 | 1.28 | 132.2
NiTe | FM | -7.42 | -5.15 | 2.27 | 4.53 | 31.2
FeC | AFM | 53.00 | 53.28 | 0.28 | 0.21 | 5.3
FeO | AFM | 5.47 | 5.46 | -0.01 | 2.48 | 20.1
MnS | AFM | 55.07 | 55.95 | 0.88 | 1.69 | 15.2
FeS | AFM | 16.13 | 16.38 | 0.25 | -0.56 | 5.7
MnSe | AFM | 50.66 | 50.65 | -0.01 | -0.20 | 10.3
FeSe | AFM | 15.30 | 15.24 | -0.06 | -1.20 | 5.8
FeTe | AFM | 31.80 | -28.89 | -60.69 | 119.63 | 35.6
CuC | NM | - | - | - | - | -
NiN | NM | - | - | - | - | -
CuN | NM | - | - | - | - | -
NiS | NM | - | - | - | - | -
CuS | NM | - | - | - | - | -
ZnS | NM | - | - | - | - | -
CuSe | NM | - | - | - | - | -
ZnSe | NM | - | - | - | - | -
ZnTe | NM | - | - | - | - | -

Our throughput computational screening revealed that out of a 48 materials in
total, only the 21 listed in Table 1 possess dynamical stability. Out of those
21, we identified five materials with FM order as lowest in energy (out of six
possibilities considered), seven materials with AFM order as the lowest-energy
one, and nine materials not exhibiting magnetic order, i.e. with magnetic
moments on individual atoms below 0.5$\mu_{B}$ in every of the six considered
magnetic configurations, and the configuration with magnetic moments exactly 0
as lowest in energy. Out of five FM materials, two of them appear to have
completely flat, square lattice when on-site Coulomb repulsion between
electrons (U) is properly included - CrC and MnC. Such a structure ensures
that the first-nearest-neighbor TM-X-TM bond comprises 90 degrees angle, and
second-nearest-neighbor 180 degrees angle, making these structures optimal -
according to GKA rules - for “maximization” of magnetic exchange in the
system. This fact served as motivation to thoroughly analyze the structural,
electronic, and magnetic properties of these two materials, and discuss the
subtle interplay of these properties that optimizes the ferromagnetic order.

The first important finding concerns the structure. Namely, as indicated
above, the relative atomic arrangement along the $z$ axis strongly depends on
the on-site Coulomb repulsion (U) between electrons on TM atoms. In Fig. 2(a),
we present the relaxed structures for CrC and MnC with no U included, and for
the U value calculated self-consistently - using the linear response theory as
introduced by Timorov et al. [36]. The latter U value is found as 3.2949 eV
for CrC, and 3.8137 eV for MnC. In case of U=0, the structure buckles, and two
X atoms relax above and below the plane of TM atoms. Furthermore, the X atoms
create the tetrahedral structure around each TM atom. On the contrary, when
using a realistic value for U, both CrC and MnC exhibit an ideally planar
structure, i.e. a 2D square lattice as desired in the GKA argumentation
towards the enhanced ferromagnetic order. For certainty, we tested the
dynamical stability of both these planar structures - and the phonon
dispersions shown in Fig. 2(b) exhibited no imaginary phonon frequencies.

Figure 2: Stability of CrC and MnC ferromagnetic monolayers. (a) The effect of
Hubbard parameter U on monolayer structures of CrC and MnC, and (b) the phonon
dispersions of two materials for optimal U shown in (a), proving dynamical
stability of the flat structures.

After determining the on-site Coulomb repulsion and the planar structural
stability, we move on to the magnetic properties of CrC and MnC. The
necessary, yet insufficient condition for long-range magnetism is the non-zero
magnetic moment per unit space. The TM atoms are expected to provide the
latter, due to the localized unpaired electrons in their 3$d$ shell, each
electron carrying spin $\frac{1}{2}$, and spin magnetic moment of
$1~{}\mu_{B}$. In most of the structures based on TMs, the contribution of the
orbital magnetic moment is negligible compared to the spin magnetic moment,
due to the quenching of orbital momentum, hence the magnetic moment can be
assumed to originate purely from the spin of the electron. Consequently, the
total magnetic moment on each TM atom is $N\times 1~{}\mu_{B}$, where $N$ is
the total number of unpaired electrons in the 3$d$ shell of each TM atom.

The basic ionic theory suggests $4+$ oxidation state of Cr and Mn cations in
our monolayers, since C atom receives 4 electrons to reach stable octet
configuration. Furthermore, Cr4+ ion has 20 electrons and the
$1s^{2}2s^{2}2p^{6}3s^{2}3p^{6}3d^{2}$ electronic configuration, while Mn4+
ion 21 electrons and the $1s^{2}2s^{2}2p^{6}3s^{2}3p^{6}3d^{3}$ electronic
configuration. This means that one expects $2~{}\mu_{B}$ per Cr atom and
$4~{}\mu_{B}$ per primitive unit cell in CrC, and $3~{}\mu_{B}$ per Mn atom
and $6~{}\mu_{B}$ per primitive unit cell in MnC. Ionic theory predicts no
moment per C atom in either cases, due to the mentioned stable octet
configuration.

For a more precise account of magnetic moments per atom we resort to DFT
calculations, and find that: (1) in case of CrC, the magnetization per Cr atom
is $2.8588~{}\mu_{B}$ and per C atom $-0.7236~{}\mu_{B}$, resulting in
$4.2703~{}\mu_{B}$ per primitive unit cell; (2) in case of MnC, the
magnetization per Mn atom is $3.8825~{}\mu_{B}$ and per C atom
$-0.7718~{}\mu_{B}$, resulting in $6.2214~{}\mu_{B}$ per primitive unit cell.
The DFT results do corroborate the crude ionic theory regarding the total
magnetization of the unit cell, but also reveal the significant hybridization
between ($d$) orbitals of TM atoms and ($p$) orbitals of C atoms - causing a
rather significant magnetization on otherwise non-magnetic C atoms.

However, in order for a system to host measurable magnetic order, next to the
non-zero magnetic moments on TM ions, it must also host significant
interaction between them. Unlike the initial estimates using the method based
on mapping between total energies of the DFT and the Heisenberg Hamiltonians,
we now calculate the magnetic exchange by mapping the energy variations due to
the infinitesimal rotation of the magnetic moment on TM atoms from the
reference FM state, between the DFT Hamiltonian in the localized-orbital basis
set, and the Heisenberg Hamiltonian
$H=\sum_{<i,j>}\textbf{S}{{}_{i}}\textbf{J}{{}_{i,j}}\textbf{S}{{}_{j}}$, as
implemented in the TB2J code [37]. In the latter Hamiltonian, $\textbf{S}_{i}$
denotes the unit 3D vector of the magnetic moment on the $i^{th}$ TM atom;
$3\times 3$ matrix $\textbf{J}_{i,j}$ stands for interaction between magnetic
moments on $i$-th and $j$-th TM atoms; and $<i,j>$ denotes $i\neq j$, with
avoided double counting. The SIA matrix cannot be calculated using this
formalism (for details see [37]), therefore it is not explicitly written in
the used Heisenberg Hamiltonian. However, SIA is not negligible, and
correspondingly must be taken into account during e.g. the calculation of
T${}_{\textrm{c}}$ for the FM order. Therefore, we combine SIA reported in
Table 1, together with the $\textbf{J}_{i,j}$ calculated using TB2J to
construct the total model Hamiltonian for eventual 2nd-principles
calculations:

$H=\sum_{<i,j>}\textbf{S}{{}_{i}}\textbf{J}{{}_{i,j}}\textbf{S}{{}_{j}}+\sum_{i}\textbf{S}{{}_{i}}\textbf{A}{{}_{i,i}}\textbf{S}{{}_{i}}.$
(1)

The main advantage of TB2J and the Green’s-function-based methodology over the
total energy mapping is the ability to orbitally disentangle the origins of
magnetic interactions [38], and also the ability to calculate interactions
between all different neighbors within a large supercell upon a single DFT
calculation on the primitive unit cell. In particular, we have calculated the
matrices $\textbf{J}_{i,j}$ up to the 284th NN for both materials.

Figure 3: The isotropic magnetic exchange, per neighboring pair. The isotropic
exchange interactions in CrC (blue) and MnC (red), calculated using Green’s
function method as implemented in TB2J. The inset depicts the numerical
labeling of the nearest-neighbor sites within the structure.

Our results obtained using TB2J obey the symmetry-imposed constraints - i.e.
all off-diagonal elements of all $\textbf{J}_{i,j}$ matrices are exactly 0.
The diagonal part of each matrix can further be split - as presented in Table
1 \- into the isotropic (J${}_{i,j}=$ J${}^{xx}_{i,j}=$ J${}^{yy}_{i,j}$) and
the anisotropic part ($\Delta=$J${}^{zz}_{i,j}-$J${}^{xx}_{i,j}$). The
anisotropic part of the exchange, $\Delta$, in either system does not exceed
few (1-6) $\mu$eV, and is comparable with the rounding error in our
calculations. Therefore, we consider $\Delta$ as effectively 0, and we
subscribe the stabilization of the magnetic order in these 2D materials to
just Ji,j and SIA ††This is in qualitative agreement with results reported in
Table 1, where $\Delta$ was found to be much smaller than SIA and Ji,j..

In Fig. 3 we present our TB2J results for Ji,j up to the 6th NN. One notices
in Fig. 3 that 1st and 2nd NN interactions are strong (few tens of meV), an
order of magnitude larger compared to usually encountered isotropic exchange
values (few meV) in 2D materials. Another observation is that even pairs over
1 nm distance have small but non-zero exchange interaction, of few hundreds of
$\mu$eV. However, since these are much smaller than the 1st and 2nd NN
interactions, in what follows, we focus on the first two NN pairs with giant
Ji,j, being essential for the high-T${}_{\textrm{c}}$ GKA ferromagnetism. We
use exchange interaction parameters calculated with TB2J, and SIA from Table 1
to build the Heisenberg Hamiltonian as explained above, to then employ Monte
Carlo simulations to evaluate the stability of the FM order with respect to
temperature. In Fig. 4, we present the thereby obtained evolution of
(normalized) magnetization (Mst/Ms), magnetic susceptibility ($\chi$), and
specific heat (Cv), as a function of temperature. One clearly sees that
estimated critical temperatures of the FM order in both CrC and MnC exceed the
room temperature, being 307 K and 428 K respectively. Even though anisotropy
is in general required to stabilize magnetic order in 2D above 0 K - in our
case that is SIA - the main reason for such a large T${}_{\textrm{c}}$ lies in
the particularly large isotropic exchange between the 1st as well as the 2nd
NN pairs of TM atoms. Table 1 provides Curie and Neel temperature values for
t-TMX structures. Notably, VN, CoN, and NiTe structures exhibit ferromagnetic
properties; however, their respective Curie temperatures are significantly
lower compared to those of CrC and MnC. This phenomenon can be primarily
attributed to two key factors: weakened exchange interactions and strong
Single Ion Anisotropy (SIA). This weakened exchange interaction and strong SIA
inherently limits their ability to maintain ferromagnetic order at elevated
temperatures, thus resulting in lower Curie temperatures when contrasted with
CrC and MnC. Conversely, antiferromagnetic materials exhibit Neel temperatures
in proximity to absolute zero (0 K) due to the significant impact of SIA. This
low Neel temperature signifies the point at which antiferromagnetic materials
undergo a transition to a non-magnetic state. The details about calculations
of Curie and Neel temperatures can be found in † Supplementary Information.

Figure 4: Thermal stability of the long-range magnetic order. Magnetization,
specific heat and magnetic susceptibility of (a) CrC and (b) MnC, as a
function of temperature.

##  III Discussion

Figure 5: Orbitally-decomposed magnetic interactions. The magnetic exchange
interactions in CrC and MnC, for the three nearest-neighbor pairs. As
illustrated in the inset, for every pair the contributions are discerned per
orbital $a$ of TM1 ($a$ = $3s,4s,3p,3d,4p$), each interacting with just one
orbital ($b$ = 3$d$) of TM2.

Having presented the core results, we next detail the origin of the magnetic
interactions behind the observed high Curie temperature in monolayer tTMXs,
bearing in mind the original assumptions following from the GKA rules. To shed
the light on the source of the large exchange Ji,j, we look into the
orbitally-resolved contributions.

Our initial assumption was that the gross of the exchange interactions
originates from the exchange between the $d$ orbitals of the interacting TM
atoms ($d$-$d$ exchange). In order to properly describe the physics behind our
observations, we treat 3$s$, 3$p$, 3$d$, 4$s$, and 4$p$ as valence orbitals on
TM atoms. Further, we quantify the contributions of each orbital-to-orbital
interaction between the neighboring TM atoms (e.g. 4$s$ on TM1 and 3$d$ on
TM2), to the total exchange between them. For facilitated interpretation, we
consider the contributed interactions between each type of orbital on TM1 ($a$
= 3$s$, 4$s$, 3$p$, 3$d$, and 4$p$) and only 3$d$ orbitals on TM2 ($b$ =
3$d$), for the first three NN magnetic interactions $J^{a,b}_{\\#NN}$, as
shown in Fig. 5. It is rather obvious from Fig. 5 that TM1(3$d$)-TM2(3$d$)
interactions dominate (being several tens of meV strong), and determine the
magnetic order in the system - in this case the FM one. Interactions between
other types of orbitals on one TM and 3$d$ orbitals on the other are generally
at least an order of magnitude smaller. Therefore, after $d$-$d$ exchange is
proven to be crucial for the large exchange and room-temperature magnetism in
tTMXs, we next decompose it into the exchanges between individual $d$-orbitals
of the interacting pair of TM atoms, to disentangle the key contributors. The
results per $d$-orbital and per material are shown in Fig. 6.

Figure 6: Sub-orbitally-decomposed $3d$-$3d$ interactions. The magnetic
exchange between five different 3d ($d_{xy}$, $d_{yz}$, $d_{z^{2}}$, $d_{xz}$
and $d_{x^{2}-y^{2}}$) orbitals on first three NN pairs, for both CrC and MnC
monolayers. CrC and MnC are indicated using dark and light shading,
respectively.

###  III.1 The nearest-neighbor interaction

Owing to the square lattice symmetry of these monolayers, our original premise
holds, and the empirical GKA rules are validated in the case of the 1st NN
exchange - the ideally 90∘ TM-X-TM bond arrangement fosters a particularly
strong FM interaction, stemming mainly from the $d$-$d$ exchange.

As seen in Fig. 6, the largest contribution to the 1st NN magnetic interaction
in both considered materials comes from the interaction between $d_{xz}$ on
one TM atom, and $d_{yz}$ on the other, together with its symmetric twin -
TM1($d_{yz}$)-TM2($d_{xz}$). We prescribe this exactly to the square-lattice
geometry of the structure and the fact that the dumbbells of the $d_{xz}$ on
TM1 and $d_{yz}$ on TM2 point along the bonds to the (X1) ligand atom between
them in the structure (analogously for $d_{yz}$ on TM1 and $d_{xz}$ on TM2,
interacting via the adjacent ligand X2). Moreover, after having a closer look
at the DFT Hamiltonian in the localized basis set, we noticed that both
$d_{xz}$ on TM1 and $d_{yz}$ on TM2 interact only with $p_{z}$ on X1, while
the other hopping matrix elements (with $p_{x}$ and $p_{y}$ on X1, and
$p_{x}$, $p_{y}$, and $p_{z}$ on X2) are 0, because of the symmetry of the
system.

Another significant contribution to the total exchange between the 1st NN pair
originates from TM1($d_{xy}$)-TM2($d_{xy}$) interactions. In that case, the
$d_{xy}$ orbitals on TM1 and TM2 point towards each other directly and have
significant direct overlap, however, there is also significant hopping between
$d_{xy}$ on both TM1 and TM2 and $p_{x}$ and $p_{y}$ orbitals on both nearest
X atoms (X1 and X2). This increased complexity of the physical picture causes
differences in the sign and the strength of those particular contributions in
the two considered compounds: in case of CrC this contribution is AFM, while
in case of MnC it is FM. The main reason for the difference in this orbital
contribution lies in the different atomic environment \- the ordering and
occupation of the atomic $d$($p$) orbitals on TM(X) atoms - and different
behavior of direct exchange in the particular environment. In case of the
atomic environment in CrC, the superexchange solely determines the orbital
contribution to the total magnetic exchange and it is AFM, while in case of
MnC, direct- and superexchange compete in such atomic environment that the
result is the FM coupling. Details of these findings are made available in the
† Supplementary Information.

Both compounds host several other non-zero interactions between different $d$
orbitals, e.g. $d_{xy}$ on TM1 and $d_{z^{2}}$ on TM2. Even though these terms
are much smaller than the dominating ones discussed above, they are still
sizable (several meV) and they do affect the total exchange, albeit on
tertiary level. In particular TM1($d_{xy}$)-TM2($d_{z^{2}}$) contribution is
very sensitive to the alteration of the superexchange hopping
TM($d_{xy}$)-X($p_{x/y}$), hence we note that superexchange mechanism plays an
important role here. For more discussion on their origin and behavior, we
refer the reader to the † Supplementary Information.

###  III.2 The next-nearest and further-neighbor interactions

Although our initial premise of strong FM interactions between the first
nearest neighbors was validated, we point out at this stage that plain GKA
rules are not sufficient to interpret the magnetic behavior of a 2D material,
even in the case of an ideally planar square-lattice structure. Namely, the
2nd NN interaction is expected to be AFM according to the GKA rules, due to
the 180∘ TM-X-TM bond alignment - but we have observed the (strong) opposite
in both materials of interest.

In what follows, we present the results for the TM1-X-TM2 bonds being aligned
with the global $x$ coordinate (as the case of TM1-X-TM2 along $y$ is
completely analogous). As was the case with the 1st NN interaction, the
orbitals aligned with the direction of TM-X-TM bonds are mainly responsible
for the large total exchange of the 2nd NN pair as well - i.e.
TM1($d_{xz}$)-TM2($d_{xz}$) interaction is the dominant one in this case. The
dumbbells of these two $d$ orbitals both point towards the common X atom, and
only interact with its $p_{z}$ orbital.

The fact that the dominating contributions to the 1st and 2nd NN interactions
are originating from the same physical process - i.e. from two $d$ orbitals on
two TM atoms that are oriented towards the common ligand atom, and interact
only with its $p_{z}$ orbital - leads towards the conclusion that they should
be of the same sign and comparable strength, as they indeed are in our results
(strongly FM). That said, the GKA rules assume the dominant contribution to
the AFM superexchange in case of 180∘ TM-X-TM bonds to be via the $p$ orbital
of X, whose dumbbell is aligned with TM-X-TM direction [40, 30] \- which would
be the $p_{x}$ orbital in the above discussion. However, even though we find
these contributions to be AFM as GKA rules would suggest, we also find that
they are an order of magnitude smaller than the dominant contributions,
involving the $p_{z}$ orbital. Therefore, the disagreement between our results
and GKA rules originates in the fact that mechanism considered dominant by GKA
(superexchange involving the $p_{x}$, or $p_{\sigma}$) is secondary in our
case, and vice versa - the mechanism considered secondary by GKA
(superexchange involving the $p_{z}$, or $p_{\pi}$) appears to be dominant in
the two materials of our interest. For interested readers, this is discussed
in more detail within the † Supplementary Information.

In case of the 3rd NN exchange, the geometry is very similar to the 1st NN
exchange, hence one expects to have the same dominant contributions. We find
that the three main contributions in both materials are indeed the same ones
as in the case of the 1st NN, i.e. $d_{xz}$-$d_{yz}$, $d_{yz}$-$d_{xz}$, and
$d_{xy}$-$d_{xy}$. However, their observed behavior is more complicated. In
case of CrC the TM1($d_{xz}$)-TM2($d_{yz}$) and TM1($d_{yz}$)-TM2($d_{xz}$)
are contributing the most, however to AFM order (positive exchange parameter),
while the second largest contribution comes from TM1($d_{xy}$)-TM2($d_{xy}$)
and it is FM. In case of MnC the main contribution comes from
TM1($d_{xy}$)-TM2($d_{xy}$) and it is AFM. The interaction between
TM1($d_{xz}$)-TM2($d_{yz}$) and TM1($d_{yz}$)-TM2($d_{xz}$) is again FM,
however their magnitude affects the total 3rd NN exchange significantly less.
Since the 3rd NN exchange is negligible compared to the 1st and 2nd NN
exchange, we will not detail these interactions. We see however that sign and
strength of the orbital contributions to the 3rd NN exchange are mainly
determined by the orbital ordering and occupation of Cr and Mn atoms - i.e.
their atomic environment - and the fact that in different environment, the
different mechanisms may be dominating. In case of the
TM1($d_{xy}$)-TM2($d_{xy}$) contribution, our results suggest that in case of
CrC there is the competition between direct and superexchange, while in MnC it
is the usual superexchange through the X ligand that dominates. By comparing
the results for the first- and the third-nearest-neighbor exchange in two
materials, one could argue that in case that atomic environment stimulates the
superexchange alone, the result will be AFM interaction (the first-nearest-
neighbor in CrC, and third-nearest-neighbor in MnC). On the other hand, when
atomic environment stimulates the competition between different exchange
mechanisms (the first-nearest-neighbor in MnC, and third-nearest-neighbor in
CrC) resulting interaction between TM1($d_{xy}$)-TM2($d_{xy}$) will be FM. For
interested readers, we provide brief additional discussion of the effect of
atomic environment to the third-nearest-neighbor magnetic exchange in the †
Supplementary Information.

##  IV Conclusions

In summary, after a computational throughput screening of the whole family of
tTMX materials, we have identified two dynamically stable and ideally 2D flat
Lieb-like magnetic crystals, CrC and MnC. According to the seminal Goodenough-
Kanamori-Anderson rules, materials of such symmetry are prone to host
pronouncedly high ferromagnetic exchange interaction. Our detailed analysis of
the magnetic properties has shown that both Cr and Mn ions indeed have sizable
magnetic moments, that materials host large isotropic magnetic interactions
(order of 10 meV), that exchange anisotropy is negligible in these systems
(order of 1 $\mu$eV), and that stabilization of the long-range magnetic order
in these materials should be prescribed to single-ion anisotropy (order of 0.1
meV). As a result of this large isotropic exchange and non-zero single-ion
anisotropy, the Curie temperature for ferromagnetic transition of these
materials exceeds the room temperature - 307 K in CrC and 428 K in MnC.

In detailed analysis, we showed how the symmetry/geometry of the system
selects the dominant contributions to the total magnetic exchange. In case of
the 1st NN and 90 degrees TM-C-TM lattice direction, the interaction between
$d_{xz}$ on TM1 and $d_{yz}$ on TM2 - TM1($d_{xz}$)-TM2($d_{yz}$) - dominates
together with its twin counterpart TM1($d_{yz}$)-TM2($d_{xz}$). We presented
the hypothesis that this interaction is mediated by the $p_{z}$ orbital of the
intervening ligand C. In case of the 2nd NN and 180 degrees TM-C-TM lattice
direction, the dominant contributions are again geometry-selected, being
either TM1($d_{xz}$)-TM2($d_{xz}$) or TM1($d_{yz}$)-TM2($d_{yz}$), depending
on whether TM-C-TM direction is along global Cartesian $x$ or $y$ coordinate.
In this case we also hypothesized that interaction is mediated by the $p_{z}$
orbital of the C atom. In case of the 3rd NN, which was of secondary
importance in our systems, we could not establish a criterion for single
dominant contribution based on geometry, however, even there two contributions
were identified as leading ones. Finally we have shown that interactions
beyond the 3rd NN are far smaller and may be neglected in any prediction of
experimentally measurable physical properties of these systems.

In conclusion, we provided convincing evidence that CrC and MnC in their two-
dimensional tetragonal phase are room-temperature Lieb-like-lattice 2D
ferromagnets. Furthermore, we highlighted the effect of geometry on the
magnetic properties - both the magnetic moments and the exchange interactions
between them. That said, we expect any emergent material with flat square-
lattice structure and magnetic atom from the first half of the first row of
transition metals (V, Cr, Mn, Fe) to: (i) have sizable magnetic moments
(mainly from unpaired $d_{xz}$ and $d_{yz}$ orbitals), and (ii) host large
geometry-selected magnetic interactions between these two types of orbitals.
Consequently, in any such system with sufficient anisotropy in magnetic
interactions, the latter features will foster the large isotropic exchange and
high critical temperature of the long-range magnetic order in two dimensions,
as already shown for 2D metallic MnB monolayer.[41]

##  V Acknowledgement

This work was supported by the Research Foundation-Flanders (FWO) and the
Technological Research Council of Turkey (TUBITAK) under Contract No. 118F512.
D.Š. is a doctoral fellow of FWO under Contract No. 11J4322N. The
computational resources and services for this work were provided by the VSC
(Flemish Supercomputer Center), funded by the FWO and the Flemish Government –
department EWI.

## References

  * Huang _et al._ [2017] B. Huang, G. Clark and E. e. a. Navarro-Moratalla, _Nature_ , 2017, 546, 270–273.
  * Gong _et al._ [2017] C. Gong, L. Li and Z. e. a. Li, _Nature_ , 2017, 546, 265–269.
  * Xing _et al._ [2017] W. Xing, Y. Chen, P. M. Odenthal, X. Zhang, W. Yuan, T. Su, Q. Song, T. Wang, J. Zhong, S. Jia, X. C. Xie, Y. Li and W. Han, _2D Materials_ , 2017, 4, 024009.
  * Lado and Fernández-Rossier [2017] J. L. Lado and J. Fernández-Rossier, _2D Materials_ , 2017, 4, 035002.
  * Webster and Yan [2018] L. Webster and J.-A. Yan, _Phys. Rev. B_ , 2018, 98, 144411.
  * Xu _et al._ [2018] C. Xu, J. Feng and H. e. a. Xiang, _npj Comput Mater_ , 2018, 4, 0.
  * Kashin _et al._ [2020] I. V. Kashin, V. V. Mazurenko, M. I. Katsnelson and A. N. Rudenko, _2D Materials_ , 2020, 7, 025036.
  * Deng _et al._ [2018] Y. Deng, Y. Yu and Y. e. a. Song, _Nature_ , 2018, 563, 94–99.
  * Jiang _et al._ [2018] S. Jiang, L. Li and Z. e. a. Wang, _Nature Nanotech_ , 2018, 13, 549–553.
  * Huang _et al._ [2018] B. Huang, G. Clark and D. e. a. Klein, _Nature Nanotech_ , 2018, 13, 544–548.
  * Bonilla _et al._ [2018] M. Bonilla, S. Kolekar and Y. e. a. Ma, _Nature Nanotech_ , 2018, 13, 289–293.
  * Roemer _et al._ [2020] R. Roemer, C. Liu and K. Zou, _npj 2D Mater Appl_ , 2020, 4, 33.
  * Song _et al._ [2022] Q. Song, C. Occhialini and E. e. a. Ergeçen, _Nature_ , 2022, 602, 601–605.
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner, _Phys. Rev. Lett._ , 1966, 17, 1133–1136.
  * Fei _et al._ [2018] Z. Fei, B. Huang and P. e. a. Malinowski, _Nature Materials_ , 2018, 17, 778–782.
  * Meng _et al._ [2021] L. Meng, Z. Zhou and M. e. a. Xu, _Nature Communications_ , 2021, 12, 809.
  * Ghosh _et al._ [2021] A. Ghosh, M. Palit, S. Maity, V. Dwij, S. Rana and S. Datta, _Phys. Rev. B_ , 2021, 103, 064431.
  * Zhang _et al._ [2022] G. Zhang, F. Guo, H. Wu, X. Wen, L. Yang, W. Jin, W. Zhang and H. Chang, _Nature Communications_ , 2022, 13, 5067.
  * Wu _et al._ [2021] H. Wu, W. Zhang, L. Yang, J. Wang, J. Li, L. Li, Y. Gao, L. Zhang, J. Du, H. Shu and H. Chang, _Nature Communications_ , 2021, 12, 5688.
  * O’Hara _et al._ [2018] D. J. O’Hara, T. Zhu, A. H. Trout, A. S. Ahmed, Y. K. Luo, C. H. Lee, M. R. Brenner, S. Rajan, J. A. Gupta, D. W. McComb and R. K. Kawakami, _Nano Letters_ , 2018, 18, 3125–3131.
  * Kan _et al._ [2014] M. Kan, S. Adhikari and Q. Sun, _Phys. Chem. Chem. Phys._ , 2014, 16, 4990–4994.
  * Boukhvalov and Politano [2020] D. W. Boukhvalov and A. Politano, _Nanoscale_ , 2020, 12, 20875–20882.
  * Chua _et al._ [2021] R. Chua, J. Zhou, X. Yu, W. Yu, J. Gou, R. Zhu, L. Zhang, M. Liu, M. B. H. Breese, W. Chen, K. P. Loh, Y. P. Feng, M. Yang, Y. L. Huang and A. T. S. Wee, _Advanced Materials_ , 2021, 33, 2103360.
  * Sun _et al._ [2017] Y. Sun, Z. Zhuo, X. Wu and J. Yang, _Nano Lett._ , 2017, 17, 2771–2777.
  * Jiang _et al._ [2018] Z. Jiang, P. Wang, J. Xing, X. Jiang and J. Zhao, _ACS Appl. Mater. Interfaces_ , 2018, 10, 39032–39039.
  * Zhang _et al._ [2019] C. Zhang, Y. Nie, S. Sanvito and A. Du, _Nano Lett._ , 2019, 19, 1366–1370.
  * Liu _et al._ [2020] L. Liu, S. Chen, Z. Lin and X. Zhang, _J. Phys. Chem. Lett._ , 2020, 11, 7893–7900.
  * Wu _et al._ [2021] D. Wu, Z. Zhuo, H. Lv and X. Wu, _J. Phys. Chem. Lett._ , 2021, 12, 2905–2911.
  * Goodenough [1955] J. B. Goodenough, _Phys. Rev._ , 1955, 100, 564–573.
  * Kanamori [1959] J. Kanamori, _Journal of Physics and Chemistry of Solids_ , 1959, 10, 87–98.
  * Anderson [1959] P. W. Anderson, _Phys. Rev._ , 1959, 115, 2–13.
  * Anderson [1963] P. W. Anderson, _Theory of Magnetic Exchange Interactions:Exchange in Insulators and Semiconductors_ , Academic Press, 1963, vol. 14, pp. 99–214.
  * Gao _et al._ [2021] L. Gao, Y.-F. Zhang and S. Du, _Nano Research_ , 2021, 14, 2826–2830.
  * Ong [2021] S. P. Ong, _Materials Project Documentation_ , 2021.
  * Šabani _et al._ [2020] D. Šabani, C. Bacaksiz and M. V. Milošević, _Phys. Rev. B_ , 2020, 102, 014457.
  * Timrov _et al._ [2022] I. Timrov, N. Marzari and M. Cococcioni, _Computer Physics Communications_ , 2022, 279, 108455.
  * He _et al._ [2021] X. He, N. Helbig, M. J. Verstraete and E. Bousquet, _Computer Physics Communications_ , 2021, 264, 107938.
  * Kashin [2020] I. V. e. a. Kashin, _2D Materials_ , 2020, 7, 025036.
  * [39] This is in qualitative agreement with results reported in Table 1, where $\Delta$ was found to be much smaller than SIA and Ji,j.
  * van Vleck [1951] J. H. van Vleck, _J. Phys. Radium_ , 1951, 262–274.
  * Jiang _et al._ [2018] Z. Jiang, P. Wang, X. Jiang and J. Zhao, _Nanoscale Horiz._ , 2018, 3, 335–341.

See pages 1, of Supp_GKA_arxiv.pdf See pages 0, of Supp_GKA_arxiv.pdf