Experimental work suggests that biological soil crusts, dominant primary producers in drylands and tundra, are particularly vulnerable to disturbances that cause reverse ecological succession. To model successional transitions in biocrust communities, we propose a resource-firing game that captures succession dynamics without specifying detailed function forms. The game is evaluated in idealized terrestrial ecosystems, where disturbances are modeled as a reduction in available resources that triggers inter-species competition. The resource-firing game is executed on a finite graph with nodes representing species in the community and a sink node that becomes active when every species is depleted of resources. First, we discuss the theoretical basis of the resource-firing game, evaluate it in the light of existing literature, and consider the characteristics of a biocrust community that has evolved to equilibrium. We then examine the dependence of resource-firing and game stability on species richness, showing that high species richness increases the probability of very short and long avalanches, but not those of intermediate length. Indeed, this result suggests that the community’s response to disturbance is both directional and episodic, proceeding towards reverse succession in bursts of variable length. Finally, we incorporate the spatial structure of the biocrust community into a Cayley Tree and derive a formula for the probability that a disturbance, modeled as a random attack, initiates a large species-death event.

§ INTRODUCTION

Global change threatens to destabilize terrestrial ecosystems by triggering ecological state transitions, significant shifts in community structure and function. These shifts may occur in response to diminished resource availability or climate-induced disturbances, including increased frequency and intensity of extreme events, altered precipitation patterns, and rising temperatures. In some cases, changes in disturbance regimes can cause productive ecosystems to transition to less productive ones [1, 2, 3]. Understanding these dynamics is crucial to inform conservation strategies and accurately assess terrestrial ecosystem health. However, predicting state transitions can be difficult—communities of primary producers, which help drive ecosystem dynamics, often exhibit non-deterministic, non-linear, and non-monotonic transition behaviors [4, 5, 6].

Biological soil crusts are communities of cyanobacteria, lichen, and bryophytes that appear ubiquitously in terrestrial ecosystems. As dominant primary producers in drylands and tundra, biocrusts account for 7 percent of terrestrial primary productivity and over 50 percent of biological nitrogen fixation [7].

Biocrusts also form assemblages of extracellular secretions and biomass filaments that regulate local hydrological cycling and channel resources from nitrogen-enriched, hydrated biocrusts to depleted plants [8, 9, 10, 11]. These products rehabilitate damaged soil, forming the first stage of colonization in ecosystems recovering from disaster. Indeed, biocrusts are expected to regulate climate impacts on up to 70 percent of global dryland surface [12]. Given these ecosystem services and extent of cover, biocrust state transitions are of key consequence.

Ecological succession is the development of community structure over time, with each state characterized by different dominant organisms [13]. In the initial stage, bare soil is colonized by pioneer cyanobacteria, which bind filamentous sheaths together to stabilize soil microstructure [14]. Over time, early stage cyanobacteria are replaced by a more diverse community of late-stage biocrusts, including lichen, cyanolichen, and bryophytes in order of succession. Lichen form from a symbiotic relationship between a photosynthetic partner and fungal components, while mosses perform photosynthesis using stem leaves. Both penetrate soil surfaces through hyphal extension [15, 16]. Biocrusts perform different ecosystem services depending on their dominant succession stage. For example, moss-dominated biocrusts sequester carbon and dissimilate nitrogen more effectively, while lichen-dominated biocrusts are better suited to metabolize aromatic compounds and exhibit higher nitrogen uptake [17].

Previous research has established that biocrusts are vulnerable to a range of disturbances, including trampling, wildfire, flooding, vehicle traffic, and livestock grazing. Global change may also exert differential effects on biocrust communities based on succession stage, compromising vital ecosystem processes [18]. Several studies have observed the effects of these disturbances on biocrusts. Chronic physical disturbance and natural disaster have been shown to alter community structure for decades, notably accelerating erosion and reducing C and N concentrations in soil [19]. More recently, climate manipulation studies have observed reduced biocrust metabolic activity and diminished network connectivity in response to experimental warming [20]. Some work suggests that all types of disturbance trigger biocrust reverse ecological succession, a state transition from late to early stages. For instance, Ferrenberg et al. (2014) reported shifts from moss cover to cyanobacteria cover in response to both experimental warming and human trampling [12]. Other studies have shown that when subjected to climate manipulation and grazing, early stage cyanobacteria dominate the disturbance gradient, while late-stage cyanobacteria dominate the recovery gradient [21, 22]. Prior to reverse succession, disturbed biocrusts may also often undergo dramatic declines in taxonomic diversity [23]. Though disturbances may vary in the magnitude and rate of reverse succession they trigger, these results suggest that biocrusts exhibit a general, disturbance-agnostic response to suboptimal environmental conditions. The consequences of these induced reverse state transitions can be severe, including deficits in global carbon and nitrogen budgets.

In this study, we develop a time-independent resource-firing game that captures succession dynamics without specifying detailed function forms. Our goals were first, to propose a general model of biocrust succession transitions based on unified theory and second, to examine dynamics in idealized terrestrial ecosystems. For this purpose, we model disturbances as a reduction in available resources, triggering competition between species within the community. The resource-firing game is evaluated on (1) its ability to predict reverse succession in response to a generic disturbance, (2) the episodic and directional character of changes in community structure, and (3) the stability of equilibria.

First, we discuss the theoretical basis of the resource-firing game and show how the model may be implemented computationally. Then, we turn to the characteristics of a biocrust community that has evolved to an equilibrium state, the various dependences of the model, and attacks on the network.

§ MODEL

§.§ Overview and Theory

Shifts in succession often exhibit the same profile when subjected to very different environmental disturbances. We compared the effect of climate manipulation treatments, physical disturbance, and disaster events on the relative cover of four biotic groups that dominate biocrusts: early-stage cyanobacteria, late-stage cyanobacteria, lichen, and bryophytes. Relative cover and geographic data from the literature are summarized in Appendix Table <ref>. In 47 of the 50 studies in the dataset, biocrusts exhibit a reverse successional response to disturbance or reduced species diversity, irrespective of the type of disturbance. As a result, a key criterion of this model is its ability to predict disturbance-agnostic reverse succession. In each of the 4 studies that did not conform to this trend, moss cover increased in response to warming, physical disturbance, or resource limitation. Notably, this phenomenon appeared in temperate semiarid regions, including Southeast Spain and the Colorado Plateau [24, 25]. After a disturbance, the order of reverse succession in biocrusts follows a predictable pattern: moss, lichen, late-stage cyanobacteria, and early-stage cyanobacteria. However, when biocrusts recover from disturbance, the order of recolonization is often reversed [26, 27]. We suggest that while different disturbances may vary in the size and frequency of their effects on community structure, succession can arise from the internal dynamics of community-wide resource allocation, rather than the triggering stressor. The resource-firing game proposed in this paper captures these disturbance-agnostic internal dynamics.

Each study in Appendix Table <ref> observes biocrust cover both prior to disturbance and after disturbance, where $t_{obs}$ is the elapsed observation period. However, few have captured intermediate succession dynamics. The available literature is sufficient to supply initial conditions and suggest equilibria, but cannot robustly constrain model trajectories at $t < t_{obs}$. Therefore, we preliminarily assume that state transition behaviours proceed monotonically and call for field work to addresses this gap in existing literature. To build our model assumptions from this dataset, we do not specify function forms of relative cover change, instead focusing on macroscopic changes in succession stage. The strategy of biocrust ecological succession is estimated as the attempt to optimally distribute finite resources throughout an interconnected biomass system. This objective is achieved by considering the entire biomass as a nexus of species and attempting to allocating resources based on their fitness [28].

§.§ Order Theory

Biocrust ecological succession is the development of its community structure over time, with each succession stage characterized by different dominant organisms. In the initial stage, bare soil is colonized by pioneer cyanobacteria. Over time, early stage cyanobacteria are replaced by late-stage biocrusts, including lichen, cyanolichen, and bryophytes in order of succession. We formalize this successional sequence as the partially ordered set $S$. Let $S$ be the chain

\begin{align}
\centering
\text{\textit{cyanobacteria}} \longrightarrow \text{\textit{lichen}} \longrightarrow \text{\textit{bryophytes}}
\end{align}

Pioneer early succession populations require fewer resources and thrive in resource-limited conditions, whereas late succession populations outcompete the former in resource-abundant conditions [26]. Let $\phi$ be a function that maps dominant succession stage to the resources it requires to sustain itself. The notion of resources, while highly abstract, might be captured by applying the species-area relationship, generalized fertility indices, or net primary productivity, the details of which are beyond the scope of this study.

\begin{equation}
\phi(s) =
\begin{cases}
r_{f} & \text{if \textit{s} = fungi} \\
r_{c} & \text{if \textit{s} = cyanobacteria} \\
r_{l} & \text{if \textit{s} = lichen} \\
r_{b} & \text{if\textit{ s} = bryophytes}
\end{cases}
\end{equation}

such that $r_f \leq r_c \leq r_l \leq r_b$. We now turn to the order structure of our game. Explicit $r_f$, $r_c$, $r_l$, and $r_b$ values depend on (1) total resources $R$ available to a community and (2) the number of species of each succession stage $s_f, s_c, s_l, s_b$. The equation for resource allocation is given by

\begin{equation} \label{e:resources}
\left[\begin{array}{c}
s_f \\
s_c \\
s_l \\
\end{array}\right]^T
\left[\begin{array}{c}
r_f \\
r_c \\
r_l \\
\end{array}\right]
= R
\end{equation}

Formally, two partially ordered sets $S$ and $Q$ are order-isomorphic ($S \cong Q$) if there exists a map $\phi$ from $S$ onto $Q$ such that $x \leq y$ in $S$ $\iff$ $\phi(x) \leq \phi(y)$ in Q [29].

That is, $\phi$ must faithfully mirror the order structure of $S$. Let $Q = \phi(s)$ be a partially ordered set of allocated resources for $\phi$ defined in <ref>. Then, $S$ and $Q$ are order-isomorphic, allowing us to directly map the successional composition of a community to its resource allocation.

\begin{align}
\centering
\phi(\text{\textit{cyanobacteria}}) \longrightarrow \phi(\text{\textit{lichen}}) \longrightarrow \phi(\text{\textit{bryophytes}})
\end{align}

§.§ Representing Biocrust Networks

A biocrust community is modeled as an interconnected network that distributes resources to several populations of species competing to maximize ecological fitness. The strategy of biocrust ecological succession is estimated as the attempt to optimally distribute finite resources throughout an interconnected biomass [28].

A complete graph is a simple undirected graph $G$ in which every pair of nodes is connected by a unique edge.

We represent a generic biocrust network as a complete graph. A complete graph with $n$ nodes has $\frac{n(n-1)}{2}$ edges. Each node represents a species in the community, and an edge represents competition for shared resources. In other words, each species in the biocrust community is assumed to be in competition with every other species. Each node has two fields. The first field is its designated succession stage, $s$: fungus, early cyanobacteria, late cyanobacteria, lichen, or bryophyte. The second field is the concomitant number of resources, $\phi(s)$. The succession stage with a plurality of designated nodes is said to dominate the community. Figure <ref> shows two complete graphs, representing a fungi- and late cyanobacteria-dominated community with 8 species and a late cyanobacteria-dominated community with 12 species. The following section describes the rules of the resource-firing game on such graphs.

[scale = 0.7]
[label=bare soil,draw,fill=white] (node1) ;
[label=[name=l] fungi,draw,fill=brown] (node2) at ([xshift=2cm]node1.east);
[label=[name=l] early cyan.,draw,fill=yellow] (node3) at ([xshift=2cm]node2.east);
[label=[name=l] late cyan.,draw,fill=lime] (node4) at ([xshift=2cm]node3.east);
[label=[name=l] lichen,draw,fill=green] (node5) at ([xshift=2cm]node4.east);
[label=[name=l] moss,draw,fill=blue] (node6) at ([xshift=2cm]node5.east);

[scale = 0.7]
[typeset=(), math nodes, nodes=draw, circle, n=8, radius=2cm, clockwise]
1[as=$ $,fill=white,opacity = 1];
2[as=$ $,fill=brown,opacity = 1];
3[as=$ $,fill=brown,opacity = 1];
4[as=$ $,fill=lime,opacity = 1];
5[as=$ $,fill=blue,opacity = 1];
6[as=$ $,fill=lime,opacity = 1];
7[as=$ $,fill=green,opacity = 1];
8[as=$ $,fill=white,opacity = 1];
subgraph K_n ;

[scale = 0.7]
[typeset=(), math nodes, nodes=draw, circle, n=12, radius=2cm, clockwise]
1[as=$ $,fill=white,opacity = 1];
2[as=$ $,fill=lime,opacity = 1];
3[as=$ $,fill=yellow,opacity = 1];
4[as=$ $,fill=lime,opacity = 1];
5[as=$ $,fill=blue,opacity = 1];
6[as=$ $,fill=lime,opacity = 1];
7[as=$ $,fill=green,opacity = 1];
8[as=$ $,fill=brown,opacity = 1];
9[as=$ $,fill=white,opacity = 1];
10[as=$ $,fill=yellow,opacity = 1];
11[as=$ $,fill=lime,opacity = 1];
12[as=$ $,fill=green,opacity = 1];
subgraph K_n ;

Biocrust networks of varying scale. Both graphs are complete, with 8 and 12 species respectively (Color Required)

§.§ Resource-Firing Game

The chip-firing game models the loss of energy in networks of individuals whose states depend on the state of their neighbors. In this paper, we define the resource-firing game, a portmanteau of chip-firing game for resource allocation. The game is executed on a finite graph $G = (V,E)$ with $|V| = n$ nodes, in addition to a sink node, which becomes active when every species in the community is depleted of resources [30].

Given any biocrust network, we can play the resource allocation game as follows. First, we add a sink node to $G$. For $G$ to remain complete, the sink node must have edges to every species in the network. Thus, every community $G$ has $n$ nodes, 1 sink node and $n-1$ species, retaining the property of completeness. Next, a non-negative number of resources are distributed to each node, with the exception of the sink node. Recall that the distribution of resources to each species as a function of its succession stage is given by Equation <ref>. We write the resource configuration as follows:

\begin{equation}
\mathbf{r} = \{r_1,r_2,...,r_k\} \in \mathds{N}_0
\end{equation}

where $r_k$ is the number of resources on node $n_k$. The basic step of this game is firing initiated by a generic disturbance, the dynamics of which are detailed in Algorithm <ref>. Firing represents a transfer of resources from one stressed species to its non-stressed neighbors. Accordingly, the node with highest $r_n$, equivalently lowest fitness, fires.

A node $n_k$ can be fired if $r_k \geq deg(n_k)$, where $deg(n_k)$ is equivalent to the number of neighbors $n_k$ has. For example, if $n_1$ has edges to $n_2, n_4,$ and $n_5$, those nodes are considered its neighbors and $deg(n_1) = 3$. The set of neighbor pairs is denoted $E(G)$. When node $n_k$ fires, it sends 1 chip to each of its neighbors, losing $deg(n_k)$ resources. Biologically, this indicates that the species occupying $n_k$ no longer has sufficient resources to sustain itself. The released resources now become available to competing species. If $n_1$ is selected to fire and initially has $5$ resources, it sends $1$ resource to each its $3$ neighbors and retains $2$ resources. As a result, the succession stage of the species occupying $n_1$ regresses to one that can be supported by 2 resources, a process which we will call reverse succession. Dually, $n_2, n_4,$ and $n_5$, $n_1$'s neighbors, gain resources and can support species with advanced succession stages, a process which we will call forward succession. The mechanisms of forward and reverse succession are beyond the scope of the present manuscript, and we refer readers to the work of others. Formally, the piece-wise function $r'(i)$ denotes the successor configuration of of any node $i$ in the graph after a given node $j$ fires:

\begin{equation}
\mathbf{r'(i)} =
\begin{cases}
r_j - deg(n_j) & \text{if } i = j\\
r_i + 1 & \text{if }(n_i,n_j) \in E(G)\\
r_i & \text{if } (n_i,n_j) \notin E(G)
\end{cases}
\end{equation}

Redistribution of resources after firing may enable new non-sink nodes to fire. The game ends when the community arrives at a stable configuration in which no nodes are eligible to fire. At this point, the sink vertex releases its chips to each of its neighbors and the game restarts. The length of the game is defined as the number of firings until a stable configuration is reached.

Because chip firing exhibits terminating behavior, the stable configuration is always unique. Chip-firing also displays local confluence, meaning nodes can be fired in any order without affecting the length or final configuration of the game [31].

Resource-Firing Algorithm

$s$ is sink node, $V \neq \emptyset$ and contains all nodes except $s$
$t_{limit}$ iterations of resource-firing on a graph with $n$ nodes and $r$ resources

$t < t_{limit}$
$n_{max} = \{n_{i} \in V$ $|$ $\forall n_{j} \in V \quad$ $r_{j} \leq r_{i}\}$

$r_{max} \geq deg(n_{max})$
$r_{max} = r_{max} - deg(n_{max})$
$ \forall (n_i, n_j) \in E(G) \quad$ $r_{j} = r_{j} + 1$
$r_{sink} = 0$
$\forall n_{j} \in V \quad$  $r_j = r_j + 1$

§ DYNAMICS IN TERRESTRIAL ECOSYSTEMS

§.§ Example: $n = 5$

Figure <ref> shows a simple firing sequence on a biocrust community with $n = 5$ nodes (4 species) and $r = 10$ resources. As the sink node accumulates resources and the number of resources available to the community decreases, it undergoes reverse succession. Conversely, the release of resources from $s$ triggers forward succession, yielding the initial configuration.

[scale = 0.7]
[label=bare soil,draw,fill=white] (node1) ;
[label=[name=l] fungi,draw,fill=brown] (node2) at ([xshift=2cm]node1.east);
[label=[name=l] early cyan.,draw,fill=yellow] (node3) at ([xshift=2cm]node2.east);
[label=[name=l] late cyan.,draw,fill=lime] (node4) at ([xshift=2cm]node3.east);
[label=[name=l] lichen,draw,fill=green] (node5) at ([xshift=2cm]node4.east);
[label=[name=l] moss,draw,fill=blue] (node6) at ([xshift=2cm]node5.east);

[node distance=4cm]
(a) at (0,0) [scale = 0.6]
[typeset=(), math nodes, nodes=draw, circle, n=5, radius=1cm, clockwise]
1[as=$s_0$,fill=white,opacity = 1,accepting]; 2[as=$1$,fill=brown,opacity = 1]; 3[as=$4$,fill=blue,opacity = 1]; 4[as=$2$, fill=lime,opacity = 1]; 5[as=$3$,fill=green,opacity = 1]; subgraph K_n ;
(b) at (3,0) [scale = 0.6]
[math nodes, nodes=draw, circle, n=5, radius=1cm, clockwise]
1[as=$s_1$,fill=white,opacity = 1,accepting]; 2[as=$2$,fill=lime,opacity = 1]; 3[as=$0$,fill=white,opacity = 1]; 4[as=$3$,fill=green,opacity = 1]; 5[as=$4$,fill=blue,opacity = 1]; subgraph K_n ;
(c) at (6,0) [scale = 0.6]
[math nodes, nodes=draw, circle, n=5, radius=1cm, clockwise]
1[as=$s_2$,fill=white,opacity = 1,accepting]; 2[as=$3$,fill=green,opacity = 1]; 3[as=$1$,fill=brown,opacity = 1]; 4[as=$4$,fill=blue,opacity = 1]; 5[as=$0$,fill=white,opacity = 1]; subgraph K_n ;
(d) at (9,0) [scale = 0.6]
[math nodes, nodes=draw, circle, n=5, radius=1cm, clockwise]
1[as=$s_3$,fill=white,opacity = 1,accepting]; 2[as=$4$,fill=blue,opacity = 1]; 3[as=$2$,fill=lime,opacity = 1]; 4[as=$0$,fill=white,opacity = 1]; 5[as=$1$,fill=brown,opacity = 1]; subgraph K_n ;
(e) at (4.5,-3.0) [scale = 0.6]
[math nodes, nodes=draw, circle, n=5, radius=1cm, clockwise]
1[as=$s_4$,fill=white,opacity = 1,accepting]; 2[as=$0$,fill=white,opacity = 1]; 3[as=$3$,fill=green,opacity = 1]; 4[as=$1$,fill=brown,opacity = 1]; 5[as=$2$,fill=lime,opacity = 1]; subgraph K_n ;

[->] (a) – node[above]  (b);
[->] (b) – node[above]  (c);
[->] (c) – node[above]  (d);
[ ->] (d.south) to [out=270,in=0] (e.east);
[ ->] (e.west) to [out=180,in=270] (a.south);

Sequence of chip firings on a biocrust community with $n=5$ node and 10 resources. The sink nodes in each configuration are denoted $s_0,s_1, s_2, s_3, s_4$, respectively. Graph 5 represents the stable configuration

The resulting stable configuration appears on the bottom. Mathematically, the stable configuration refers to a configuration of resources in which each node has at most as many chips as its degree. Furthermore, there is no subset of nodes that can fire. The resource-firing game will always proceed to a stable configuration and remain there forever, irrespective of how nodes fire. In other words, the stable configuration represents an equilibrium or fixed point.

Now we turn to the stability of the configurations in our $n=5$ example. Let $c_0$ represent the fixed point of $f'(r'(i)) = (r'(1),r'(2),r'(3),r'(4))$, where i is the node index. Given that $f^{(n)} \rightarrow c_0$, we conclude that $c_0$ is a stable fixed point. For the given initial configuration, $f$ has a four-cycle,

\begin{equation}
f^{(4)}(c_0) = c_0, f'(c_0) = f''(c_0) = f^{(3)}(c_0) \neq c_0.
\end{equation}

We might also attempt to determine the stability of $c_0$ by using the Tutte polynomial.

For an undirected graph $G=(V,E)$, the Tutte polynomial is $T_{G}(x,y)=\sum \nolimits _{{A\subseteq E}}(x-1)^{{k(A)-k(E)}}(y-1)^{{k(A)+|A|-|V|}}$, where $k(A)$ denotes the number of connected components of the graph $(V,A)$ [32].

A minimum spanning tree is a tree that includes every node of a given graph and the minimal set of edges needed to connect each node without cycles [33].

Evaluating the Tutte polynomial of a graph $T_G$ at $T_G(1,1)$ yields the number of spanning trees, which provides information about stable configurations. [34]. The Tutte polynomial of $G_5$ in Figure <ref> is given in Figure <ref>. Note that the coefficients are non-negative.

major grid style=dotted,
y domain=-3:3,
zmin=-10, zmax=30,
title=$x^3 + y^3 + 3x^2 + 4xy + 3y^2 + 2x + 2y$,
3[mesh] x^3 + y^3 + 3*x^2 + 4*x*y + 3*y^2 + 2*x + 2*y;

Tutte Polynomial of the $n=5$ game shown in Figure <ref>

§.§ Species Richness Dependence

In the following two sections, we consider additive factors to the basic model introduced in Section <ref>. We begin with a first approximation of the species richness of global terrestrial ecosystems compiled from available literature data in Appendix Table <ref>.

Desert     $0<n<30$

Steppe     $30<n<50$

Temperate Forest     $50<n<100$

First, we examine how the dependence of the resource-firing on the number of species ($n$) affects game stability. A common approach derived from topological graph theory to measure chip-firing stability is to compute the genus of the graph.

The genus of a graph is the minimum integer $k$ such that the graph can be topologically embedded onto a sphere with $k$ handles without crossing itself [35].

If the graph embedding of $G$ is topologically simple, $G$ has a low genus and the game may likely requires fewer resources to be available at the outset. Dually, if $G$ has a high genus, the game may require a larger number of resources [35, 36]. In Figure <ref>, we compute the genus of resource-firing games for $0 < n < 100$ for idealized desert, steppe, and temperate forest biocrust communities. One consequence of this analysis is that genus increases quadratically with respect to $n$. Ergo, the resource stability threshold of idealized forest communities may be significantly higher than those of desert and semi-desert communities.

[scale = 0.8]
only marks,
xmin=0, xmax=100,
ymin=0, ymax=5005,

legend pos=north west,
legend cell align=left,
[blue!20] (0,0) rectangle (32,5000);
[red!20] (32,0) rectangle (52,5000);
[green!20] (52,0) rectangle (102,5005);
[mark=o,blue] coordinates
[mark=+,red]  coordinates
[mark=triangle,teal] coordinates
Temperate Forest;

Genus by Ecosystem (Color Required)

Next, we consider how species richness impacts the avalanche distribution of the resource-firing game. When a node fires, a stressed species releases resources to non-stressed neighbors. This reallocation of resources may enable other species to release resources, resulting a percolation of firings across the graph. An avalanche is precisely this trickling sequence of firings. The length of the avalanche is the length of the sequence, that is, the number of firings that occur until firing must be induced for the game to proceed. The avalanche distribution is the probability distribution of the lengths of avalanches that occur, given some initial configuration of chips. For idealized desert, steppe, and temperate forest communities with $n=15$, $n=40$, and $n=75$ species respectively, the resource-firing game was initiated 10,000 times and the avalanche distribution in Figure <ref> generated.

[scale = 0.8]
only marks,
xmin=0, xmax=5,
ymin=0, ymax=10,

legend pos=north west,
legend cell align=left,
[mark=o,blue] coordinates
Desert ($n=15$);
[mark=+,red] coordinates
Steppe ($n=40$);
[mark=triangle,teal] coordinates
(1.75, 0)
Temperate Forest ($n=75$);

Avalanche Distribution by Ecosystem (Color Required)

The distribution suggests that reverse succession is both episodic and directional, proceeding towards reverse succession in bursts of variable length. In all ecosystems, avalanches are frequently short or long, but not often medium-sized. As $n$ increases, two effects are observed: the probability of medium-sized avalanches tends towards zero and the average length of avalanches increases.

§.§ Attacks on Biocrust Communities

The resource-firing game assumes a uniformly competitive interaction between species. In reality, species compete most aggressively with species in their proximity for local pools of resources. For a complete graph $G$ representing a biocrust community, we define an accompanying path graph—a Cayley tree—that approximates the spatial structure of the community. By Cayley's formula, for $G$ with $n$ nodes, $n^{n-2}$ unique spanning trees exist [37]. To build our Biocrust Cayley Tree (BCT), we replicate all $n$ nodes on $G$, then place edges on them such that a desired spanning tree is generated. The structure of the BCT is as follows. There is a singular node $v$ located at the root at the tree with edges to $k$ neighbors. These neighbors comprise the first shell of the tree. Each node in the first shell has $k-1$ neighbors, forming the second shell. Each node in every subsequent shell also has $k-1$ neighbors [38].

[scale = 0.4]
(0,0) node[disc=c-0-1];
/in 3/disc,3/disc,3/disc,3/disc
[white] (c-0-1) circle(pt);
÷in 1,...,
(c-0-1) ++(÷*(360/)-180/:pt) node[=c--÷];
(c--) – (c--÷);

Cayley Tree with $n=45$ and $k = 3$

A random attack on the BCT compromises community structure by arbitrarily removing one node at a time. Attacks represent critical disturbances, such as climate change, trampling, or altered precipitation. Three types of nodes are defined: alive, incapacitated, and dead. If a node is attacked, the species it represents dies, and as the attack percolates through the community, its descendants are critically incapacitated. The species represented by all other nodes are considered alive. Let $p_A = \frac{n_A}{n}$ be the fraction of alive species in the community, where $n$ is the total number of nodes in the BCT. We define $A$ as the event that the species selected for attack is alive and $D$ as the event that it is already dead. Given layers $L$ and the number of neighbors $k$, Equation <ref> approximates the conditional probability of $A$ given $\bar{D}$ [38]. Note that a very large $L$ is assumed.

\begin{equation}\label{eqn:PA}
P(A|\bar{D}) = \frac{(k-2)(1-p_A)^{L+1}}{(k-2) - (k-1)p_A}
\end{equation}

Based on their degree of severity, critical disturbances can lead to species death events of various sizes. We denote such events $E$. We estimate the probability of $E$ with a death toll of at least $k$ as

\begin{equation}
P(E) = \frac{L(k-1)}{k-2} \cdot P(L_D)
\end{equation}

where $P(L_D)$ is the probability that the disturbance is severe enough to attack the BCT, that is, cause species death [38].

§ CONCLUSIONS

We have presented the resource-firing game as a model of resource allocation in complex networks, with a particular focus on its application to biocrust communities. First, the model was motivated by experimental results. Then, we discussed its theoretical basis and the characteristics of a biocrust community that has evolved to an equilibrium state.

We note here that the initial formulation of the resource-firing game hinges on some simplifying assumptions. For example, it assumes species are equally likely to compete with one another. In real-world ecosystems, species interaction also depends on differences in behavior or spatial proximity. Additionally, the model assumes that interactions are either present or absent, though they may exist on a spectrum of strength or frequency. However, the simple rules of our formulation adequately predict reverse succession in response to a generic disturbance, the episodic and directional character of changes in community structure, and the stability of equilibria.

One key result of our analysis is that changes in species richness have nonlinear effects on network structure and the avalanche distribution of the resource-firing game. In particular, the resource stability threshold for idealized forest communities is significantly higher than that of desert and semi-desert communities, and high species richness increases the probability of very short and long avalanches, but not those of intermediate length. These results suggests that the community’s response to disturbance is both directional and episodic, proceeding towards reverse succession in bursts of variable length.

We also introduce the BCT as an approximation of the spatial structure of a biocrust community and derive probabilities of critical disturbances that can compromise community structure.

Reference     Ecosystem Type     Disturbance     Data Reported

Campbell (1979) [39]      Semi-desert     Sand burial     Other parameters

Johansen et al. (1984) <cit.>     Desert     Fire     Relative cover

Johansen and St. Clair (1986) [41]     Desert     Grazing     Relative cover

Harper and Marble (1988) [42]     –     Vascular Plants     General trends

Johansen et al. (1993) [43]     Steppe     Fire     Species richness

Kaltenecker et al. (1997) [44]     Steppe     Fire     Species richness

Eldridge and Koen (1998) [45]     Desert     Trampling, Fire     Relative cover

Belnap and Phillips (2001) [46]     Grassland     Grass invasion     Species richness

Belnap and Eldridge (2001) [26]     –     –     General trends

Belnap and Lange (2001) [47]     –      Frequent fires     General trends

Belnap et al. (2004) [48]     Desert     Rain     Other parameters

Myers and Davis (2003) [49]     Temperate forest      Fire     Species richness

Hilty et al. (2004) [50]     Steppe     Fire     Species richness

Bowker et al. (2005) [51]     Semi-desert     Micronutrients     Relative cover

Thompson et al. (2005) [52]     Desert     Soil fertility     Relative cover

Belnap et al. (2006) [53]     Semi-desert     Grass invasion     Relative cover

Callaway (2007) [54]     Desert     Micronutrients     Species richness

Williams et al. (2008) [55]     Desert     Grazing, drought     Relative cover

Budel et al. (2008) [56]     Desert     Trampling, Rain     Other parameters

Lazaro et al. (2008) [57]     Semi-desert     Micronutrients     Other parameters

Kidron et al. (2008) [58]     Desert     Grazing     Other parameters

Wang et al. (2009) [59]     Desert     Sand burial     Relative cover

Li et al. (2010) [60]     Desert     Geomorphology     Other parameters

Langhans et al. (2010) [61]     Temperate forest     Fine-disturbance     General trends

Kuske et al. (2011) [19]     Semi-desert     Trampling     Species richness

Li et al. (2011) [62]     Desert     Ant nests     Other parameters

Escolar et al. (2012) [25]     Desert     Warming, Rain     Relative cover

Reed et al. (2012) [63]     Semi-desert     Rain     Relative cover

Johnson et al. (2012) [64]     Semi-desert     Warming, Rain     Other parameters

Zelikova et al. (2012) [65]     Semi-desert     Warming, Rain     Relative cover

Yeager et al. (2012) [66]     Semi-desert     Warming     Other parameters

Maestre et al. (2013) [67]     Semi-desert     Warming     Relative cover

Concostrina-Zubiri et al. (2013) [68]     Desert     Micronutrients     Species richness

Concostrina-Zubiri et al. (2014) [69]     Desert     Grazing     Species richness

Zhang et al. (2014) [70]     Desert     Sand burial     Other parameters

Ferrenberg et al. (2015) [12]     Semi-desert     Trampling     Relative cover

Maestre et al. (2015) [71]     Desert     Warming     Relative cover

Steven et al. (2015)[72]     Semi-desert     Trampling     Species richness

Budel et al. (2016) [73]     –     –     General trends

Belnap et al. (2016) [26]     –     –     General trends

Rutherford et al. (2017) [74]     Semi-desert     Warming     Relative cover

Fernandes et al. (2017) [75]     Desert     Rain     Species richness

Li et al. (2017) [76]     Desert     Rain, Topsoil     Species richness

Zhang et al. (2017) [77]     Desert     Sand burial     Species richness

Eldridge et al. (2018) [78]     Desert     Warming     Species richness

Ladron de Guevara et al. (2018) [79]     Desert     Warming     Relative cover

Li et al. (2018) [80]     Desert     Warming     Relative cover

Aanderud et al. (2019) [81]     Desert     Fire     Species richness

Berdugo et al. (2021) [82]     Temperate forest     Warming     Other parameters

Li et al. (2021) [83]     Desert     Warming     Relative cover

Lan et al. (2021) [84]     Desert     Sand burial     Other parameters

Finger-Higgens et al. (2022) [85]     Semi-desert     Grass invasion     Relative cover

Ecological succession of biocrust covers from available literature data (Appendix)

[1]
Reyer, C. P., Bathgate, S., Blennow, K., Borges, J. G., Bugmann, H., Delzon, S., ... and Hanewinkel, M. (2017). Are forest disturbances amplifying or canceling out climate change-induced productivity changes in European forests?. Environmental Research Letters: ERL, 12(3), 034027. https://doi.org/10.1088/1748-9326/aa5ef1

[2]
Shaver, G. R., Canadell, J., Chapin, F. S., Gurevitch, J., Harte, J., Henry, G., ... and Rustad, L. (2000). Global Warming and Terrestrial Ecosystems: A Conceptual Framework for Analysis. Bioscience, 50(10), 871-882. https://doi.org/10.1641/0006-3568(2000)050[0871:GWATEA]2.0.CO;2

[3]
Hitz, S., and Smith, J. (2004). Estimating global impacts from climate change. Global Environmental Change, 14(3), 201-218. https://doi.org/10.1016/j.gloenvcha.2004.04.010

[4]
McClintock, B. T., Langrock, R., Gimenez, O., Cam, E., Borchers, D. L., Glennie, R., and Patterson, T. A. (2020). Uncovering ecological state dynamics with hidden Markov models. Ecology letters, 23(12), 1878-1903. https://doi.org/10.1111/ele.13610

[5]
Scheffer, M., Carpenter, S., Foley, J.A., Folke, C. and Walker, B. (2001). Catastrophic shifts in ecosystems. Nature, 413, 591. https://doi.org/10.1038/35098000

[6]
Fasiolo, M., Pya, N. and Wood, S.N. (2016). A comparison of inferential methods for highly nonlinear state space models in ecology and epidemiology. Stat. Sci., 31, 96–118.

[7]
Elbert, W., Weber, B., Burrows, S., Steinkamp, J., Büdel, B., Andreae, M. O., and Pöschl, U. (2012). Contribution of cryptogamic covers to the global cycles of carbon and nitrogen. Nature Geoscience, 5(7), 459-462. https://doi.org/10.1038/ngeo1486

[8]
Costa, O. Y., Raaijmakers, J. M., and Kuramae, E. E. (2018). Microbial extracellular polymeric substances: ecological function and impact on soil aggregation. Frontiers in Microbiology, 9, 1636. https://doi.org/10.3389/fmicb.2018.01636

[9]
Carvajal Janke, N., and Coe, K. K. (2021). Evidence for a fungal loop in shrublands. Journal of Ecology, 109(4), 1842-1857. https://doi.org/10.1111/1365-2745.13610

[10]
Dettweiler‐Robinson, E., Sinsabaugh, R. L., and Rudgers, J. A. (2020). Fungal connections between plants and biocrusts facilitate plants but have little effect on biocrusts. Journal of Ecology, 108(3), 894-907. https://doi.org/10.1111/1365-2745.13310

[11]
Rudgers, J. A., Dettweiler‐Robinson, E., Belnap, J., Green, L. E., Sinsabaugh, R. L., Young, K. E., and Darrouzet‐Nardi, A. (2018). Are fungal networks key to dryland primary production?. American Journal of Botany, 105(11), 1783-1787. https://doi.org/10.1002/ajb2.1184

[12]
Ferrenberg, S., Reed, S. C., and Belnap, J. (2015). Climate change and physical disturbance cause similar community shifts in biological soil crusts. Proceedings of the National Academy of Sciences, 112(39), 12116-12121. https://doi.org/10.1073/pnas.1509150112

[13]
Prach, K., and Walker, L. R. (2011). Four opportunities for studies of ecological succession. Trends in Ecology and Evolution, 26(3), 119-123. https://doi.org/10.1016/j.tree.2010.12.007

[14]
Belnap, J., and Gardner, J. S. (1993). Soil microstructure in soils of the Colorado Plateau: the role of the cyanobacterium Microcoleus vaginatus. The Great Basin Naturalist, 40-47.

[15]
Rosentreter, R., Eldridge, D. J., Westberg, M., Williams, L., and Grube, M. (2016). Structure, composition, and function of biocrust lichen communities. Biological soil crusts: An organizing principle in drylands, 121-138. https://doi.org/10.1007/978-3-319-30214-0_7

[16]
Niinemets, Ü., and Tobias, M. (2013). Scaling light harvesting from moss “leaves” to canopies. Photosynthesis in bryophytes and early land plants, 151-171. https://doi.org/10.1007/978-94-007-6988-5_9

[17]
Tian, C., Pang, J., Bu, C., Wu, S., Bai, H., Li, Y., ... and Siddique, K. H. (2022). The microbiomes in lichen and moss biocrust contribute differently to carbon and nitrogen cycles in arid ecosystems. Microbial Ecology, 1-12. https://doi.org/10.1007/s00248-022-02077-7

[18]
Zaady, E., Eldridge, D. J., and Bowker, M. A. (2016). Effects of local-scale disturbance on biocrusts. Biological soil crusts: An organizing principle in drylands, 429-449. https://doi.org/10.1007/978-3-319-30214-0_21

[19]
Kuske, C. R., Yeager, C. M., Johnson, S., Ticknor, L. O., and Belnap, J. (2012). Response and resilience of soil biocrust bacterial communities to chronic physical disturbance in arid shrublands. The ISME Journal, 6(4), 886-897. https://doi.org/10.1038/ismej.2011.153

[20]
Dou, W., Xiao, B., Yao, X., and Kidron, G. J. (2023). Asymmetric responses of biocrust respiration to precipitation manipulation under a changing semiarid climate. Geoderma, 430, 116318. https://doi.org/10.1016/j.geoderma.2022.116318

[21]
Sun, J., and Li, X. (2021). Role of shrubs in the community assembly of biocrusts: the biotic and abiotic influences along a biocrust succession gradient. Plant and Soil, 460, 163-176.

[22]
Darby BJ, Neher DA, and Belnap J (2007) Soil nematode communities are ecologically more mature beneath late- than early-successional stage biological soil crusts. Applied Soil Ecology 35:203–212. https://doi.org/10.1016/j.apsoil.2006.04.006

[23]
Antoninka, A., Bowker, M. A., Chuckran, P., Barger, N. N., Reed, S., and Belnap, J. (2018). Maximizing establishment and survivorship of field-collected and greenhouse-cultivated biocrusts in a semi-cold desert. Plant and Soil, 429, 213-225. https://doi.org/10.1007/s11104-017-3300-3

[24]
Maestre, F. T., Escolar, C., Bardgett, R. D., Dungait, J. A., Gozalo, B., and Ochoa, V. (2015). Warming reduces the cover and diversity of biocrust-forming mosses and lichens, and increases the physiological stress of soil microbial communities in a semi-arid Pinus halepensis plantation. Frontiers in microbiology, 6, 865. https://doi.org/10.3389/fmicb.2015.00865

[25]
Escolar, C., Martínez, I., Bowker, M. A., and Maestre, F. T. (2012). Warming reduces the growth and diversity of biological soil crusts in a semi-arid environment: implications for ecosystem structure and functioning. Philosophical Transactions of the Royal Society B: Biological Sciences, 367(1606), 3087-3099. https://doi.org/10.1098/rstb.2011.0344

[26]
Belnap, J., Weber, B., and Büdel, B. (2016). Biological soil crusts as an organizing principle in drylands (pp. 3-13). Springer International Publishing. https://doi.org/10.1007/978-3-319-30214-0_1

[27]
Colesie, C., Felde, V. J. M. N. L., and Büdel, B. (2016). Composition and macrostructure of biological soil crusts. Biological soil crusts: An organizing principle in drylands, 159-172. https://doi.org/10.1007/978-3-319-30214-0_9

[28]
Odum, E. P. (1969). The Strategy of Ecosystem Development: An understanding of ecological succession provides a basis for resolving man's conflict with nature. Science, 164(3877), 262-270. https://doi.org/10.1126/science.164.3877.262

[29]
Davey, B. A., and Priestley, H. A. (2002). Introduction to lattices and order. Cambridge University Press. https://doi.org/10.1017/cbo9780511809088.003

[30]
Klivans, C. J. (2018). The mathematics of chip-firing. CRC Press. https://doi.org/10.1201/9781315206899

[31]
Newman, M. H. A. (1942). On theories with a combinatorial definition of" equivalence". Annals of Mathematics, 223-243. https://doi.org/10.2307/1968867

[32]
Ellis-Monaghan, J. A., and Merino, C. (2011). Graph polynomials and their applications I: The Tutte polynomial. Structural analysis of Complex Networks, 219-255. https://doi.org/10.1007/978-0-8176-4789-6_9

[33]
Cheriton, D., and Tarjan, R. E. (1976). Finding minimum spanning trees. SIAM Journal on Computing, 5(4), 724-742.

[34]
Welsh, D. (1999). The tutte polynomial. Random Structures & Algorithms, 15(3‐4), 210-228. https://doi.org/10.1002/(SICI)1098-2418(199910/12)15:3/4<210::AID-RSA2>3.0.CO;2-R

[35]
Nordhaus, E. A., Stewart, B. M., and White, A. T. (1971). On the maximum genus of a graph. Journal of Combinatorial Theory, Series B, 11(3), 258-267.

[36]
McDonough, A. (2018). Genus from sandpile torsor algorithm. arXiv preprint, arXiv:1804.07807. https://doi.org/10.46298/dmtcs.6176

[37]
Harary, F. (Ed.). (2015). A seminar on graph theory. Courier Dover Publications.

[38]
Patron, A. (2022). Avalanche-size distribution of Cayley tree. arXiv preprint arXiv:2210.14078.

[39]
Campbell, S. (1979). Soil stabilization by a prokaryotic desert crust: implications for Precambrian land biota. Origins of Life, 9(4), 335-348. https://doi.org/10.1007/978-94-009-9085-2_8

[40]
Johansen, J. R., St. Clair, L. L., Webb, B. L., and Nebeker, G. T. (1984). Recovery patterns of cryptogamic soil crusts in desert rangelands following fire disturbance. Bryologist, 238-243. https://doi.org/10.2307/3242798

[41]
St. Clair, L. L., and Johansen, J. R. (1993). Introduction to the symposium on soil crust communities. The Great Basin Naturalist, 1-4. https://doi.org/10.5194/bg-2017-178

[42]
Harper, K. T., and Marble, J. R. (1988). A role for nonvascular plants in management of arid and semiarid rangelands. Vegetation science applications for rangeland analysis and management, 135-169. https://doi.org/10.1007/978-94-009-3085-8_7

[43]
Johansen, J. R., Ashley, J., and Rayburn, W. R. (1993). Effects of rangefire on soil algal crusts in semiarid shrub-steppe of the lower Columbia Basin and their subsequent recovery. The Great Basin Naturalist, 73-88.

[44]
Kaltenecker, J. H. (1997). The recovery of microbiotic crusts following post-fire rehabilitation on rangelands of the western Snake River Plain. Boise State University Theses and Dissertations. 1072.

[45]
Eldridge, D. J., and Koen, T. B. (1998). Cover and floristics of microphytic soil crusts in relation to indices of landscape health. Plant Ecology, 137, 101-114. https://doi.org/10.1023/A:1008036214140

[46]
Belnap, J., and Phillips, S. L. (2001). Soil biota in an ungrazed grassland: response to annual grass (Bromus tectorum) invasion. Ecological Applications, 11(5), 1261-1275.

[47]
Belnap, J. (2001). Biological soil crusts: structure, function, and management (Vol. 150, pp. 241-261). O. L. Lange (Ed.). Berlin: Springer. https://doi.org/10.1890/1051-0761(2001)011[1261:sbiaug]2.0.co;2

[48]
Belnap, J., Phillips, S. L., and Miller, M. E. (2004). Response of desert biological soil crusts to alterations in precipitation frequency. Oecologia, 141, 306-316. https://doi.org/10.1007/s00442-003-1438-6

[49]
Myers, P., and Edward Davis, J. (2003). Recolonization of soils by algae in a northcentral Florida pine forest after controlled fire and soil sterilization. Nova Hedwigia, 207-219. https://doi.org/10.1127/0029-5035/2003/0076-0207

[50]
Hilty, J. H., Eldridge, D. J., Rosentreter, R., Wicklow-Howard, M. C., and Pellant, M. (2004). Recovery of biological soil crusts following wildfire in Idaho. Journal of Range Management, 57(1), 89-96. https://doi.org/10.2458/azu_jrm_v57i1_hilty

[51]
Bowker, M. A., Belnap, J., Davidson, D. W., and Phillips, S. L. (2005). Evidence for micronutrient limitation of biological soil crusts: importance to arid‐lands restoration. Ecological Applications, 15(6), 1941-1951. https://doi.org/10.1890/04-1959

[52]
Thompson, D. B., Walker, L. R., Landau, F. H., and Stark, L. R. (2005). The influence of elevation, shrub species, and biological soil crust on fertile islands in the Mojave Desert, USA. Journal of Arid Environments, 61(4), 609-629. https://doi.org/10.1016/j.jaridenv.2004.09.013

[53]
Belnap, J., Phillips, S. L., and Troxler, T. (2006). Soil lichen and moss cover and species richness can be highly dynamic: the effects of invasion by the annual exotic grass Bromus tectorum, precipitation, and temperature on biological soil crusts in SE Utah. Applied Soil Ecology, 32(1), 63-76. https://doi.org/10.1016/j.apsoil.2004.12.010

[54]
Callaway, R. M. (2007). Positive interactions and interdependence in plant communities. Springer Science and Business Media. https://doi.org/10.1007/978-1-4020-6224-7

[55]
Williams, W. J., Eldridge, D. J., and Alchin, B. M. (2008). Grazing and drought reduce cyanobacterial soil crusts in an Australian Acacia woodland. Journal of Arid Environments, 72(6), 1064-1075. https://doi.org/10.1016/j.jaridenv.2007.11.017

[56]
Büdel, B., Darienko, T., Deutschewitz, K., Dojani, S., Friedl, T., Mohr, K. I., and Weber, B. (2009). Southern African biological soil crusts are ubiquitous and highly diverse in drylands, being restricted by rainfall frequency. Microbial ecology, 57, 229-247. https://doi.org/10.1007/s00248-008-9449-9

[57]
Lázaro, R., Cantón, Y., Solé-Benet, A., Bevan, J., Alexander, R., Sancho, L. G., and Puigdefábregas, J. (2008). The influence of competition between lichen colonization and erosion on the evolution of soil surfaces in the Tabernas badlands (SE Spain) and its landscape effects. Geomorphology, 102(2), 252-266. https://doi.org/10.1016/j.geomorph.2008.05.005

[58]
Kidron, G. J., Vonshak, A., and Abeliovich, A. (2008). Recovery rates of microbiotic crusts within a dune ecosystem in the Negev Desert. Geomorphology, 100(3-4), 444-452. https://doi.org/10.1016/j.geomorph.2008.01.012

[59]
Wang, W., Liu, Y., Li, D., Hu, C., and Rao, B. (2009). Feasibility of cyanobacterial inoculation for biological soil crusts formation in desert area. Soil Biology and Biochemistry, 41(5), 926-929. https://doi.org/10.1016/j.soilbio.2008.07.001

[60]
Li, X. R., He, M. Z., Zerbe, S., Li, X. J., and Liu, L. C. (2010). Micro‐geomorphology determines community structure of biological soil crusts at small scales. Earth Surface Processes and Landforms, 35(8), 932-940. https://doi.org/10.1002/esp.1963

[61]
Langhans, T. M., Storm, C., and Schwabe, A. (2010). Regeneration processes of biological soil crusts, macro-cryptogams and vascular plant species after fine-scale disturbance in a temperate region: recolonization or successional replacement?. Flora-Morphology, Distribution, Functional Ecology of Plants, 205(1), 46-60. https://doi.org/10.1016/j.flora.2008.12.001

[62]
Li, X. R., Jia, R. L., Chen, Y. W., Huang, L., and Zhang, P. (2011). Association of ant nests with successional stages of biological soil crusts in the Tengger Desert, Northern China. Applied Soil Ecology, 47(1), 59-66. https://doi.org/10.1016/j.apsoil.2010.10.010

[63]
Reed, S. C., Coe, K. K., Sparks, J. P., Housman, D. C., Zelikova, T. J., and Belnap, J. (2012). Changes to dryland rainfall result in rapid moss mortality and altered soil fertility. Nature Climate Change, 2(10), 752-755. https://doi.org/10.1038/nclimate1596

[64]
Johnson, S. L., Kuske, C. R., Carney, T. D., Housman, D. C., Gallegos‐Graves, L. V., and Belnap, J. (2012). Increased temperature and altered summer precipitation have differential effects on biological soil crusts in a dryland ecosystem. Global Change Biology, 18(8), 2583-2593. https://doi.org/10.1111/j.1365-2486.2012.02709.x

[65]
Zelikova, T. J., Housman, D. C., Grote, E. E., Neher, D. A., and Belnap, J. (2012). Warming and increased precipitation frequency on the Colorado Plateau: implications for biological soil crusts and soil processes. Plant and Soil, 355, 265-282. https://doi.org/10.1007/s11104-011-1097-z

[66]
Yeager, C. M., Kuske, C. R., Carney, T. D., Johnson, S. L., Ticknor, L. O., and Belnap, J. (2012). Response of biological soil crust diazotrophs to season, altered summer precipitation, and year-round increased temperature in an arid grassland of the Colorado Plateau, USA. Frontiers in Microbiology, 3, 358. https://doi.org/10.3389/fmicb.2012.00358

[67]
Maestre, F. T., Escolar, C., de Guevara, M. L., Quero, J. L., Lázaro, R., Delgado‐Baquerizo, M., and Gallardo, A. (2013). Changes in biocrust cover drive carbon cycle responses to climate change in drylands. Global Change Biology, 19(12), 3835-3847. https://doi.org/10.1111/gcb.12306

[68]
Concostrina‐Zubiri, L., Martínez, I., Rabasa, S. G., and Escudero, A. (2014). The influence of environmental factors on biological soil crust: from a community perspective to a species level approach. Journal of Vegetation Science, 25(2), 503-513. https://doi.org/10.1111/jvs.12084

[69]
Concostrina-Zubiri, L., Huber-Sannwald, E., Martínez, I., Flores, J. F., Reyes-Agüero, J. A., Escudero, A., and Belnap, J. (2014). Biological soil crusts across disturbance–recovery scenarios: effect of grazing regime on community dynamics. Ecological Applications, 24(7), 1863-1877. https://doi.org/10.1890/13-1416.1

[70]
Zhang, B., Zhou, X., and Zhang, Y. (2015). Responses of microbial activities and soil physical-chemical properties to the successional process of biological soil crusts in the Gurbantunggut Desert, Xinjiang. Journal of Arid Land, 7, 101-109. https://doi.org/10.1007/s40333-014-0035-3

[71]
Maestre, Fernando T., et al. Warming reduces the cover and diversity of biocrust-forming mosses and lichens, and increases the physiological stress of soil microbial communities in a semi-arid Pinus halepensis plantation. Frontiers in Microbiology 6 (2015): 865. https://doi.org/10.3389/fmicb.2015.00865

[72]
Steven, B., Kuske, C. R., Gallegos-Graves, L. V., Reed, S. C., and Belnap, J. (2015). Climate change and physical disturbance manipulations result in distinct biological soil crust communities. Applied and Environmental Microbiology, 81(21), 7448-7459. https://doi.org/10.1128/aem.01443-15

[73]
Büdel, B., Dulić, T., Darienko, T., Rybalka, N., and Friedl, T. (2016). Cyanobacteria and algae of biological soil crusts. Biological soil crusts: an organizing principle in drylands, 55-80. https://doi.org/10.1007/978-3-319-30214-0_4

[74]
Rutherford, W. A., Painter, T. H., Ferrenberg, S., Belnap, J., Okin, G. S., Flagg, C., and Reed, S. C. (2017). Albedo feedbacks to future climate via climate change impacts on dryland biocrusts. Scientific Reports, 7(1), 1-9. https://doi.org/10.1038/srep44188

[75]
Fernandes, V. M., Machado de Lima, N. M., Roush, D., Rudgers, J., Collins, S. L., and Garcia‐Pichel, F. (2018). Exposure to predicted precipitation patterns decreases population size and alters community structure of cyanobacteria in biological soil crusts from the Chihuahuan Desert. Environmental Microbiology, 20(1), 259-269. https://doi.org/10.1111/1462-2920.13983

[76]
Li, X. R., Song, G., Hui, R., and Wang, Z. R. (2017). Precipitation and topsoil attributes determine the species diversity and distribution patterns of crustal communities in desert ecosystems. Plant and Soil, 420, 163-175. https://doi.org/10.1007/s11104-017-3385-8

[77]
Zhang, B., Zhang, Y., Li, X., and Zhang, Y. (2018). Successional changes of fungal communities along the biocrust development stages. Biology and Fertility of Soils, 54, 285-294. https://doi.org/10.1007/s00374-017-1259-0

[78]
Eldridge, D. J., and Delgado‐Baquerizo, M. (2019). The influence of climatic legacies on the distribution of dryland biocrust communities. Global Change Biology, 25(1), 327-336. https://doi.org/10.1111/gcb.14506

[79]
Ladrón de Guevara, M., Gozalo, B., Raggio, J., Lafuente, A., Prieto, M., and Maestre, F. T. (2018). Warming reduces the cover, richness and evenness of lichen‐dominated biocrusts but promotes moss growth: insights from an 8 yr experiment. New Phytologist, 220(3), 811-823. https://doi.org/10.1111/nph.15000

[80]
Li, X. R., Jia, R. L., Zhang, Z. S., Zhang, P., and Hui, R. (2018). Hydrological response of biological soil crusts to global warming: A ten‐year simulative study. Global Change Biology, 24(10), 4960-4971.  https://doi.org/10.1111/gcb.14378

[81]
Aanderud, Z. T., Bahr, J., Robinson, D. M., Belnap, J., Campbell, T. P., Gill, R. A., and St. Clair, S. (2019). The burning of biocrusts facilitates the emergence of a bare soil community of poorly-connected chemoheterotrophic bacteria with depressed ecosystem services. Frontiers in Ecology and Evolution, 7, 467. https://doi.org/10.3389/fevo.2019.00467

[82]
Berdugo, M., Mendoza-Aguilar, D. O., Rey, A., Ochoa, V., Gozalo, B., García-Huss, L., and Maestre, F. T. (2021). Litter decomposition rates of biocrust-forming lichens are similar to those of vascular plants and are affected by warming. Ecosystems, 24(6), 1531-1544.

[83]
Li, X., Hui, R., Zhang, P., and Song, N. (2021). Divergent responses of moss-and lichen-dominated biocrusts to warming and increased drought in arid desert regions. Agricultural and Forest Meteorology, 303, 108387. https://doi.org/10.1016/j.agrformet.2021.108387

[84]
Lan, S., Thomas, A. D., Rakes, J. B., Garcia‐Pichel, F., Wu, L., and Hu, C. (2021). Cyanobacterial community composition and their functional shifts associated with biocrust succession in the Gurbantunggut Desert. Environmental Microbiology Reports, 13(6), 884-898. https://doi.org/10.1111/1758-2229.13011

[85]
Finger-Higgens, R., Duniway, M. C., Fick, S., Geiger, E. L., Hoover, D. L., Pfennigwerth, A. A., and Belnap, J. (2022). Decline in biological soil crust N-fixing lichens linked to increasing summertime temperatures. Proceedings of the National Academy of Sciences, 119(16), e2120975119. https://doi.org/10.1073/pnas.2120975119

§ DECLARATIONS

§.§ Acknowledgements

The author thanks Erwan Monier, who provided feedback on technical aspects of the model.

§.§ Funding

The author declares that no funds, grants, or other support were received during the preparation of this manuscript.

§.§ Competing Interests

The author declares no competing interests.


# Remarkably enhanced Curie temperature in monolayer CrI3 by hydrogen and
oxygen adsorption: A first-principles calculations

Maedeh Rassekh1,3<EMAIL_ADDRESS>Junjie He2 Saber Farjami Shayesteh3
Juan Jose Palacios 1,4<EMAIL_ADDRESS>1Departamento de Física de la
Materia Condensada, Universidad Autónoma de Madrid, 28049 Madrid, Spain.
2Bremen Center for Computational Materials Science, University of Bremen, Am
Fallturm 1, 28359 Bremen, Germany
3Department of Physics, University of Guilan, 41335-1914 Rasht, Iran
4 Instituto Nicolás Cabrera (INC) and Condensed Matter Physics Center
(IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain.

###### Abstract

Two-dimensional (2D) materials unique properties and their promising
applications in post-silicon microelectronics have attracted significant
attention in the past decade. Recently, ferromagnetic order with out-of-plane
easy axis in a monolayer of CrI3 has been observed and reported, with a Curie
temperature of 45 Kelvin. Here we study, using density functional theory (DFT)
calculations, how hydrogen and oxygen adsorption affects the structural,
electronic, and magnetic properties of a CrI3 monolayer. Our results show that
while the structure remains almost unchanged by the adsorption of hydrogen,
adsorption of oxygen completely distorts it. We have also found that both the
adsorption of hydrogen and oxygen atoms significantly influences the
electronic and magnetic properties of the monolayer. While hydrogen quenches
the magnetic moments of Cr atoms, oxygen introduces an impurity band in the
gap. Interestingly, we find a strong enhancement of the Curie temperature by
full hydrogenation, while the results are not conclusive for O. This result
suggests a simple and effective approach to manipulate the electronic and
magnetic properties of 2D magnets for spintronics applications.

††preprint: APS/123-QED

##  I Introduction

Since the discovery of graphene in 2004, there is a tremendous interest in
experimental and theoretical studies as well as development of two-dimensional
(2D) materials Novoselov  _et al._ (2004). Various 2D materials such as
hexagonal boron nitride Dean  _et al._ (2010), silicene Cahangirov  _et al._
(2009), and transition metal dichalcogenides Mak  _et al._ (2010) have been
identified as important candidates for next generation electronic and
optoelectronic devices due to their unique properties. Although almost all
types of electronic behavior such as metallic, semiconducting, insulating
(topological and trivial), and superconducting have been observed and studied,
magnetism has remained elusive until very recently. The discovery of the first
truly 2D magnet, CrI3, has completed the family of 2D materials. A single
atomic layer of CrI3 is magnetic and its Curie temperature is about 45∘K, what
can be considered a relatively high temperature Huang  _et al._ (2017). Still
far from room temperature, we can hope that reaching this goal will not be
impossible, which would be of unlimited potential in nanoscale magnetic memory
devices and spintronics where spin-based logic operations instead of charge-
based ones could guarantee low power consumption and extended device
applications Wang  _et al._ (2019). Technically, one typically needs two
magnetic metals or half-metal electrodes and a non-magnetic semiconductor in
between or, alternatively, a magnetic semiconductor between non-magnetic
electrodes to have a spintronic device. The recent reports of magnetic order
in different 2D crystals such as $Fe_{3}GeTe_{3}$ Fei  _et al._ (2018),
$VSe_{2}$ Bonilla  _et al._ (2018), $Cr_{2}Ge_{2}Te_{6}$ Gong  _et al._
(2017), $FePS_{3}$ Lee  _et al._ (2016); Wang  _et al._ (2016a), and MXenes He
_et al._ (2019a) provide us with an exciting new platform in research of 2D
materials and exploration of fundamental theory of magnetism.

Chromium trihalides, CrX3 (X = Cl, Br, I), are a class of van der Waals
bonded, layered semiconductors. These compounds are magnetic and have been
known for many decades McGuire  _et al._ (2015). Weak van der Waals (vdW)
interaction between stacking layers in this kind of materials offer the
possibility for exfoliation into single-layer nanosheets Wang  _et al._
(2016b) and, remarkably, their intrinsic ferromagnetism can persist down to
few layers and even single layer. Among this family, the exceptional
ferromagnetic behavior and the highest magnetic ordering temperature of CrI3
can be ascribed to strong covalent interaction between the Cr atom and its
nearest-neighbor I atoms. Owing to its electronic band gap, CrI3 monolayer is
a promising material for spintronics applications. This is remarkable because,
typically, ferromagnets tend to be quite conductive, while
antiferromagnets/ferrimagnets are mostly insulators or semiconductors Gao  _et
al._ (2019).

All this has prompted a number of studies Frisk  _et al._ (2018); Klein  _et
al._ (2018); Jiang  _et al._ (2018). For example, based on density functional
theory (DFT) calculations, Jiang et al. Jiang  _et al._ (2018) found that a
direct-to-indirect band gap transition can occur by manipulating spin
direction of CrI3 from out-of-plane to in-plane, which can induce a magnetic
field controlled photoluminescence. Jiang et al. Jiang  _et al._ (2019) have
studied CrI3 zigzag nanoribbons, which are strips of CrI3 monolayer. They
found that these are also ferromagnetic, that the edge states dominate the
band structure around the Fermi level, and that this can be tuned by the
different edge atomic structures. Moreover, because of the weak van der Waals
(vdW) interlayer coupling, designing heterostructures with interesting
properties using CrI3 and other vdW materials is possible without considering
lattice mismatch. As found by Zhang et al. Zhang  _et al._ (2018), reducing
the distance between graphene and CrI3 in a heterostructure to about 3.3 and
2.4 Å, a Chern insulating state can be acquired. A very strong gap opening of
about 150 meV is also found in graphene which is induced by the magnetic
proximity of CrI3. Song et al. Song  _et al._ (2018) have reported that
atomically thin chromium triiodide (CrI3) can act as a spin-filter tunnel
barrier sandwiched between graphene contacts. Based on their result,
increasing the number of CrI3 layers can enhance tunneling magnetoresistance.
The tunneling current ($I_{t}$) also has a strong magnetic field dependence
which implies a spin-dependent tunneling probability related to the field-
dependent magnetic structure of CrI3.

Strain engineering is another effective strategy for tuning and controlling 2D
materials as well as CrI3. The electronic band gap of CrI3 along with
characteristic ferromagnetism can be controlled by uniaxial and biaxial
compressive as well as tensile strains Mukherjee  _et al._ (2019). Wu et al.
Wu  _et al._ (2019) have investigated magnetic and electronic properties of
monolayer CrI3 under strain. They have found that biaxial strain can induce a
ferromagnetic-antiferromagnetic transition. By varying the amount of strain, a
transition in its electronic state occurs from magnetic-metal to half-metal
and to half-semiconductor. Zheng et al. Zheng  _et al._ (2018) results suggest
that the out of plane magnetic orientation can flip to in plane by a tensile
strain. Webster et al. Webster and Yan (2018) have shown that the Curie
temperature of CrI3 monolayer is sensitive to the strain, and its highest
amount is close to the equilibrium point. As the compressive strain increases
the Curie temperature decreases.

Figure 1: Schematic illustration of the single-layer CrI3 unit cell: bond
length of Cr-I (a), bond angles and vertical distance between the plane
containing I atoms and the plane containing Cr atoms (b), top view of CrI3
crystalline structure, showing the honeycomb arrangement of the Cr atoms (c),
side view of a single Cr site with an arrow representing its out-of-plane
magnetic moment and the splitting of d-levels into a higher energy $e_{g}$
doublet and a lower energy $t_{2g}$ triplet in an octagonal environment (d),
and the edge-sharing octahedral Cr3+ ion with six I- ions (e).

Tuning the electronic and magnetic properties of materials has always been a
central topic in condensed matter physics and material science. For instance,
in this regard, interactions between adatoms and CrI3 monolayer are important
because the electronic structure of the system can be altered, which could
lead to appealing new properties. However, the effect of adsorption of atoms
on the magnetism of intrinsic 2D magnets is not yet sufficiently understood.
In this work, motivated by previous investigations on other 2D materials, we
intend to reveal the effect of H and O adsorption on the structural,
electronic, and magnetic properties of CrI3 monolayer. Our paper is structured
in the following manner. Details of computational methodology are described in
Sec. 2. We discuss our results in Sec. 3. Finally, conclusions and a summary
of our results are given in Sec. 4.

##  II computational Methods and preliminary tests

We have performed DFT Hohenberg and Kohn (1964); Kohn and Sham (1965)
calculations using the OpenMX code Ope , which is based on norm-conserving
pseudo-potential method Bachelet (1982); Troullier and Martins (1991);
Vanderbilt (1990); Blöchl (1990); Morrison  _et al._ (1993) with a partial
core correction and a linear combinations of pseudo atomic orbitals (LCPAO) as
a basis functions Ozaki (2003); Ozaki and Kino (2004) here specified by
Cr$6.0-s^{2}p^{2}d^{1}$, I$7.0-s^{2}p^{2}d^{1}$, H$6.0-s^{2}p^{1}$, and
O$6.0-s^{2}p^{2}d^{1}$. For example, in the case of the Cr atom,
Cr$6.0-s^{2}p^{2}d^{1}$ means that the cutoff radius is 6.0 Bohr Ozaki (2003);
Ozaki and Kino (2004), and that two primitive orbitals for each of s and p
components and one primitive orbital for d components are used. The exchange-
correlation functional was the spin-polarized GGA-PBE Perdew  _et al._ (1996).
All calculations were performed until the change in total energy between two
successive iteration steps, converged to less than $10^{-6}$ Hartree. For the
purpose of the calculation of the density of states (DOS), a Gaussian smearing
of the energy levels was applied with standard deviation set to 0.15 eV. A
CrI3 unit cell with two Cr and six I atoms was initially constructed using
experimental values McGuire  _et al._ (2015). A vacuum spacing of 20 Å in the
z direction is used to prevent the interaction of the single-layer CrI3 with
its periodic images. We converged the total energy with respect to the cutoff
energy and kgrid for the unit cell. Based on these convergence tests, a cutoff
energy of 220 Ry and $11\times 11$ kgrid have been used in all presented
results. The unit cell structure was also geometrically optimized. The quasi-
Newton eigenvector following method was executed for structural relaxation of
all structures until the change in forces between two successive iteration
steps were less than $1\times 10^{-3}$ Hartree/Bohr.

After optimization we get a Cr-I covalent bond length of 2.73 Å with bond
angles of 92.88∘, 89.55∘, and 87.12∘. These values are consistent with
previous DFT results Wang  _et al._ (2016b); Mukherjee  _et al._ (2019); Ghosh
_et al._ (2019). As shown in (Fig. 1), the Cr layer is sandwiched in between
two I layers in the two-dimensional CrI3 sheet. The plane containing I atoms
is at a vertical distance of 1.58 Å from the plane containing Cr atoms (Fig.
1) which form a honeycomb network. The Cr3+ ions are surrounded by 6 nearest-
neighbor I- ions each of them bonded to two Cr atoms and arranged in an edge
sharing octahedral coordination (Fig. 1). In an ionic picture, Cr3+ has three
electrons, with an electronic configuration of $3d^{3}4s^{0}$. Owing to the
octahedral symmetry, the d-orbitals split into three low-lying $t_{2g}$ and
two high-lying $e_{g}$ orbitals Skomski (2008). Therefore the three electrons
of Cr3+ ions in this environment occupy the $t_{2g}$ with $S=3/2$. Our spin-
polarized DFT calculations reveal an energy difference $(E_{\rm AFM}-E_{\rm
FM})$ between ferromagnetic (FM) and antiferromagnetic (AFM) order of 15.81
meV/Cr atom, meaning that the FM state of CrI3 monolayer is energetically more
favorable than the AFM state and in good agreement with previous works Wang
_et al._ (2016b); Ghosh  _et al._ (2019); Guo  _et al._ (2018).

Hydrogen and O adsorption on CrI3 monolayer is modeled using single adatom in
a $2\times 2$ supercell containing 8 Cr and 24 I atoms. For this supercell we
re-converged in kgrid and, based on these convergence tests, a $7\times 7$
kgrid has been finally used throughout. We have considered adsorption of these
atoms at 6 different sites: hollow (H) site at the center of star, hollow (h)
site at the center of diamond, bridge (B) site at the midpoint of the Cr-I
bond, top ( TI) site directly above the upper I atom, valley ( VI) site
directly above the lower I atom, and top ( TCr) site directly above the Cr
atom (Fig. 2). Unless otherwise mentioned, for each adsorption site, the
adatom and CrI3 atoms are relaxed in all x, y and z directions. To estimate
the adsorption energy of adatoms, the calculations for the isolated adatom,
isolated CrI3 single layer and CrI3 with adatom system are performed in same-
sized CrI3 supercell.

Figure 2: Symbolic presentation of six adsorption sites on CrI3 lattice.

##  III Results and discussion

First we have confirmed that the electronic structure of CrI3 monolayer is in
good agreement with prior studies Wang  _et al._ (2016b); Lado and Fernández-
Rossier (2017); Zhang  _et al._ (2015). Our calculations show that the single-
layer of CrI3 is a ferromagnetic indirect semiconductor with a band gap of
about 1.1 eV. The magnetic moment is concentrated mostly at the Cr atoms. The
total magnetic moments are 6$\mu_{B}$ and 24$\mu_{B}$ for the CrI3 unit cell
and the $2\times 2$ supercell, respectively, which is about $3\mu_{B}$ per Cr
atom. The total density of states (TDOS) and the band structure for pure CrI3
monolayer are plotted in Fig. 3.

Figure 3: TDOS and band structure of pristine CrI3 monolayer (a), and PDOS
(b). The positive (negative) values correspond to spin-up (spin-down) states
(b). Figure 4: Effect of coordination matching on orbital projections of the
CrI3 octahedron. Schematic of rotation of CrI3 octahedron (a), and
corresponding projected DOS in CrI3 octahedron before (b) and after (c)
coordination matching. Before rotation, crystal structure aligned by lattice
symmetry: b along y axis, c along z axis, results in unclear projections.
After rotation, the Cr-I bonds are exactly aligned parallel to Cartesian axis,
allowing DOS orbitally resolved. Blue and purple balls represent Cr and I
atoms, respectively.

The imbalance between the spin-up and spin-down components is evidence of the
magnetic nature of the system. The valence band and conduction band edges are
formed by spin-up bands, showing a spin-polarized half-semiconductor
character. (The band gap for the spin-down channel is about 1.8 eV.) This is
also reflected in the projected density of states (PDOS), as shown in Fig.
3-b, where the top of the valence band and bottom of the conduction band are
formed mostly by spin-up p orbitals of I atoms and spin-up d orbitals of Cr
atoms. We can also see the hybridization between Cr-3d and I-5p states. From
the PDOS related to the d orbitals of the Cr atom, it is also clear that
occupied Cr-3d orbitals are mainly found in the spin-up direction, as we
described above using crystal field theory. For more details, we also
calculated the PDOS for individual Cr 3d orbitals (see Fig. 4-b).

Table 1:  Calculated values of various quantities for H and O adsorbed CrI3: The most stable site, adsorption energy $E_{a}$, the magnetic moment $\mu$ of the adatom/CrI3 system, the magnetic moment of an isolated adatom $\mu_{A}$, and the charge transfer $\Delta Q$ between adatom and the CrI3i monolayer, respectively. Coverage | Site |  $E_{a}$(eV) | $\mu(\mu_{B})$ | $\mu_{A}(\mu_{B})$ |  $\Delta Q$(e)
---|---|---|---|---|---
HCrI3 ($2\times 2$) | $T_{Cr}$ | 1.02 | 23 | 1 | +0.01
OCrI3 ($2\times 2$) | $V_{I}$ | 3.35 | 24 | 2 | +0.37
Figure 5: Optimized $2\times 2$ CrI3 supercell with a H atom adsorbed: top
view and nearest Cr-I bonds length (a) and H-Cr bond length (b). Optimized
$2\times 2$ CrI3 supercell with an O atom attached: top view and nearest O-Cr
bonds length (c) and side view and O-I bond length (d). Blue, purple, pink,
and red balls represent Cr, I, H and O atoms, respectively. Figure 6:
Schematic view of adsorption energy in six adsorption sites: HCrI3 (left), and
OCrI3 (right) systems.

It can be seen that the five Cr 3d states split into three groups
($a_{1}=d_{z^{2}}$, $e_{1}=d_{xy}$ and $d_{x^{2}-y^{2}}$, $e_{2}=d_{xz}$ and
$d_{yz}$), which is trigonal prismatic crystal field splitting instead of the
expected octahedral one. Following Hu et al. Hu  _et al._ (2017), this can be
explained by the mismatch between the Cartesian coordination with the internal
coordination of the octahedral crystal field, making the PDOS orbitally
unresolved Gao  _et al._ (2019). To enable an orbitally resolved PDOS that
distiguinshes $t_{2g}$ ($d_{xy}$, $d_{xz}$, and $d_{yz}$) from $e_{g}$
($d_{x^{2}-y^{2}}$, and $d_{z^{2}}$) orbitals, a manipulation is needed to
align the Cr-I bonds in the octahedrons along x, y, and z axis. For this we
first find the proper rotation matrix using the VESTA program VES . Then we
calculated the corresponding rotation angles and rotated the CrI3 unit cell
using Atomsk Ato . After this change, as we can see in Fig. 4-c, the PDOS is
orbitally resolved and the Cr 3d orbital splits into two groups, $e_{g}$ and
$t_{2g}$, which is the expected octahedral splitting. We can also refer this
situation to a distorted octahedral crystal field He  _et al._ (2017). For the
edge-sharing I atoms, their p orbitals split into two $p_{x}/p_{y}$ and
$p_{z}$ groups. As we can see in Fig. 4-c, the stronger $pd\sigma$
hybridization between in plane $p_{x}/p_{y}$ orbitals and the Cr
$d_{x^{2}-y^{2}}$ orbital results in lower energies compared to the weaker
$pd\pi$ hybridization between out of plane $p_{z}$ and Cr $d_{xz}/d_{yz}$
orbitals Wang  _et al._ (2016b). In an ionic picture, we can also see that the
I-3p orbitals now are almost completely filled.

Now we turn our attention to the modifications induced by H and O adatoms. The
results for the most stable adsorption site for H and O are given in (Table
1), in which HCrI3 and OCrI3 stand for H- and O-modified systems,
respectively. For each system, six adsorption sites (high symmetry sites) were
considered. The relative stabilities of these sites are quantified by their
adsorption energy $E_{a}$, defined as the the energy we need to remove the
adatom from the CrI3 monolayer: $E_{a}=(E_{\rm ad}+E_{\rm CrI_{3}})-E_{\rm
tot}$, where $E_{\rm ad}$, $E_{\rm CrI_{3}}$, and $E_{\rm tot}$ represent the
total energies of a single free adatom, the clean CrI3, and the XCrI3 system
(X = H,O), respectively. For the case with H, the $T_{Cr}$ site is found to be
the most stable structure. As can be seen in Fig. 6, all adsorption sites are
possible due to their positive absorption energies, except for hollow site.
For B , TI, and VI sites, we fixed x and y directions of H adatom to prevent
it from moving off these sites, and the remaining dimensions are relaxed. For
the OCrI3 system, the most stable configuration is the VI site. Optimized
$2\times 2$ CrI3 supercell with a H atom adsorbed is shown in Figs. 5-a and b.
The $d_{\rm HCr}$ stands for the distances between the H atom and its nearest
Cr atom. In this system, the nearest Cr atom is just below the adatom. By the
adsorption of the H atom the nearest Cr-I bond length stretches from $\sim
2.73$ Å to $\sim 2.76$ Å and 2.77 Å. While the H atom hardly influences the
original CrI3 structure, this is greatly distorted by the adsorbed O atom
(Fig. 5-c and d). In our definition, adsorption energies for all OCrI3
structures are positive (see Fig. 6), which indicates that attaching oxygen is
always energetically more favorable. From Table 1 it can be seen that the
adsorption sites are strongly dependent on the species of adatoms with an
adsorption energy for the OCrI3 system visibly higher than for HCrI3 system.

For all calculated TDOS and PDOS we have rotated the structures as we did for
the pure CrI3 monolayer. We also checked that the results (adsorption energy,
magnetic moment, bond length, …) of the optimized rotated and not-rotated
structure are all the same. Again, the only reason for this rotation is to get
the appropriate splitting. The TDOS of HCrI3 and OCrI3 systems are plotted in
Fig. 7-c along with that of the clean CrI3 for comparative purposes. As we can
see in this figure, H leaves the TDOS of CrI3 almost without change, while O
introduces visible changes. Most importantly, the latter introduces an
impurity state spin-up band that appears in the gap (see the band structure in
Fig. 7-b).

Figure 7: Charge density difference with the isosurface value of 0.005
$e/Bohr^{3}$. The purple and pink represent the charge accumulation and
depletion (a), from left to right, the band structures for pure CrI3, HCrI3,
and OCrI3 in $2\times 2$ systems, respectively. The spin-up and spin-down
bands are represented by black and red lines, respectively (b), and TDOS of
HCrI3 and OCrI3 systems. The positive (negative) values correspond to spin-up
(spin-down) states. Black lines show The TDOS of the pure CrI3 (c). Figure 8:
TDOS and PDOS of HCrI3 ($2\times 2$) system OCrI3 ($2\times 2$) system (a),
the positive (negative) values correspond to spin-up (spin-down) states. A
schematics of the Lewis structure before and after the H and O atoms
adsorption in CrI3 (b).

On the other hand, the Mulliken population analysis Mulliken (1955) yields an
excess charge of electrons on H and O adatoms (Fig. 7-a). The obtained results
of charge transfer are given in Table 1. The electronegativity of H, O, Cr and
I are 2.20, 3.44, 1.66 and 2.66 respectively. When a Cr atom is exposed to H
atom, due to the higher electronegativity of H, charge is transferred from Cr
to H. In the case of O, the amount of charge transfer is more, due to its even
higher electronegativity.

The DOS projected on different atomic species for the HCrI3 ($2\times 2$)
system is presented in Fig. 8-a. (The TDOS is for comparative purposes.)
Clearly, we can see the hybridization between the 1s state of H atom and the
3d states of the nearest Cr atom. For more details we have also calculated the
orbital projected PDOS for Cr 3d orbitals (Fig. 9-a).

Figure 9: DOS projected onto Cr-d and H-s orbitals for HCrI3 system (a), and
onto the Cr-d , I-p, and O-p orbitals for OCrI3 system (b).

The down-spin $t_{2g}$ d orbitals are now partially filled, but the octahedral
splitting is preserved. As we mentioned, the $t_{2g}$ levels of Cr in an
octahedral environment have 3 unpaired electrons, totalling $3\mu_{B}$ per
atom. As a result, when Cr binds to H, it shares one of its d-electrons (one
of the three unpaired d electrons pairs with the H s orbital electron) and the
total magnetization decreases by $1\mu_{B}$. In brief, due to the reduction of
a valence electron per H atom, the overall magnetization is expected to
decrease upon H concentration, but probably not completely since it is
difficult to imagine a total passivation of Cr atoms by 3 H atoms each.

O and I electronic configurations are $2s^{2}\ 2p^{4}$ and $5s^{2}\ 5p^{5}$,
respectively. Therefore, based on Hund’s rule, O and I have 2 and 1 unpaired
electrons in their p orbitals, respectively. When the O atom is exposed to the
lower I atom in the CrI3 single-layer, it helps the I atom to break the bonds
with the Cr atoms while making new bonds with the displaced I atom and with
the Cr atoms by itself. As a result, it shares its two unpaired electrons, and
the entire system magnetization remains unchanged. Fig. 8-b illustrates the
schematic Lewis structure that occurs. The PDOSs of the OCrI3 ($2\times 2$)
system is shown in Fig. 8-a. The O p states hybridize with Cr d and I p states
(as they show some identical peaks). To shed more light onto the details of
this hybridization, we have calculated the orbital projected DOS in this case
too (Fig. 9-b). Oxygen has completely destroyed the octahedral symmetry. The d
states are not degenerate in two groups anymore. In the case of O and I, the
main hybridization is between $p_{x}/p_{z}$ orbitals of the two atoms.

Having established the electronic changes induced by the adsorbed species, the
physical and measurable consequences are now addressed. One of the key
measurable quantities that characterizes a ferromagnetic material is its Curie
temperature ($T_{C}$) at which the ferromagnetic-paramagnetic transition
occurs. A common way to study this transition is by mapping our DFT
Hamiltonian to a Heisenberg spin Hamiltonian Smart (1966) from which a single
atom Hamiltonian can be written as:

$\mathcal{H}_{1}=-2\bm{S}_{i}.\sum_{j}{J_{ij}\bm{S}_{j}},$ (1)

where S represents the net magnetic moment at the $N_{i}$ site, i and j stand
for the nearest $N_{i}$ atoms with exchange-coupling parameter J. Replacing
the interaction with other spins by an effective magnetic field
$\textbf{H}_{eff}$ one can obtains:

$\mathcal{H}_{1}=-\vec{\mu}_{i}.\bm{H}_{eff}=-g\mu_{B}\bm{S}_{i}.\bm{H}_{eff},$
(2)

where $g$ is the Landé factor, $\mu_{B}$ is the Boltzmann constant, and

$\bm{H}_{eff}=\frac{2}{g\mu_{B}}\sum_{j}{J_{ij}\bm{S}_{j}}$ (3)

Based on the Weiss approximation, each $\textbf{S}_{j}$ can be replaced by its
average value
$\left\langle\textbf{S}_{j}\right\rangle=\left\langle\textbf{S}\right\rangle$

$\bm{H}_{eff}=\frac{2\left\langle\textbf{S}\right\rangle}{g\mu_{B}}\left(\sum_{j}{J_{ij}}\right)$
(4)

By assumption, all magnetic atoms are identical and equivalent, this implies
that $\left\langle{\textbf{S}_{i}}\right\rangle$ is related to the
magnetization of the crystal by
$\textbf{M}=ng\mu_{B}\left\langle{\textbf{S}}\right\rangle$ where n is the
number of spins per unit volume.

$\bm{H}_{eff}=\frac{2}{ng^{2}\mu^{2}_{B}}\left(\sum_{j}{J_{ij}}\right)\bm{M}=\lambda\bm{M}$
(5)

where $\lambda$ is the Weiss molecular field coefficient. If there are only
nearest-neighbor interaction, $\sum_{j}{J_{ij}}=zJ$ where z is the number of
nearest-neighbors. Then

$\lambda=\frac{2zJ}{ng^{2}\mu^{2}_{B}}$ (6)

If there is an external field $\textbf{H}_{ext}$, then the total field acting
on the $i$th spin is $\textbf{H}=\textbf{H}_{ext}+\textbf{H}_{eff}$. Having
$M=\chi_{C}H$ where $\chi_{C}$ is Curie susceptibility,

$M=\chi_{C}(\textbf{H}_{ext}+\textbf{H}_{eff})=\chi_{C}\textbf{H}_{ext}+\chi_{C}\lambda\textbf{M}$
(7)

since $\chi={\partial M}/{\partial H}$,

$\chi=\frac{\chi_{C}}{1-\lambda\chi_{C}}$ (8)

putting $\chi_{C}=C/T$, where C is Curie constant
($C={n(g\mu_{B})^{2}S(S+1)}/{3K_{B}}$) Reif (2009), gives the Curie-Weiss Law:

$\chi(T)=\frac{C}{T-T_{C}}$ (9)

where $T_{C}=\lambda C$. By substituting $\lambda$ and C in this equation

$T_{C}^{\rm MFT}=\lambda C=\frac{2zJS(S+1)}{3k_{B}}$ (10)

Thus, based on the Weiss molecular-field theory (MFT), the $T_{C}$ can be
simply estimated using the nearest-neighbor exchange-coupling parameter $J$,
being $z=3$ the number of nearest-neighbor Cr atoms in the CrI3 monolayer,
$S=3/2$, and $k_{B}$ is the Boltzmann constant He  _et al._ (2017).

We have computed exchange coupling constants using $jx$, a post-processing
code for OpenMX , Han  _et al._ (2004); Terasawa  _et al._ (2019). We have
first computed $J$ for pristine CrI3, obtaining the nearest-neighbor value and
corresponding Curie temperature that are presented in Table 2. Considering now
the H adatom case, we have computed the Cr-Cr nearest-neighbor coupling
constants between all Cr atoms and taken the average, obtaining a higher value
and a higher Curie temperature, as indicated in Table 2. Repeating this
procedure for the case of O, we obtain, however, a smaller Curie temperature
of 55 K. Therefore, the adsorption of both species, even at low coverage,
translates into a significant change of the Curie temperature with respect to
that of pristine CrI3.

Table 2:  The (average) nearest-neighbor exchange-coupling parameter $J$ and the Curie temperature ($T_{C}$) based on the Weiss molecular-field theory using the $jx$ code and Ising model to compute the coupling constants , along with the rescaled absolute values. |  CrI3 |  HCrI3 |  OCrI3
---|---|---|---
$J({\rm meV}){\rm(jx)}$ | 0.70 | 0.90 | 0.64
$T_{C}(^{\circ}K)\ {\rm(jx)}$ | 61 | 78 | 55
$T_{C}(^{\circ}K)\ {\rm(jx/rescaled)}$ | - | 57 | 41
$J({\rm meV})\ {\rm(Ising)}$ | 0.95 | 1.02 | 1.06
$T_{C}(^{\circ}K)\ {\rm(Ising)}$ | 83 | 89 | 92
$T_{C}(^{\circ}K)\ {\rm(Ising/rescaled)}$ | - | 48 | 50

The MFT $T_{C}$ for the H/O functionalized CrI3 monolayer can be corrected by
an appropriate rescaling Liu  _et al._ (2018); He  _et al._ (2019b)

$\frac{T_{C}\rm{(XCrI_{3})}}{T_{C}^{\rm MFT}\rm{(XCrI_{3})}}=\frac{T_{C}^{\rm
exp}{\rm(CrI_{3})}}{T_{C}^{\rm MFT}{\rm(CrI_{3})}}$ (11)

where $X$ represents H or O, $T_{C}{\rm(XCrI_{3})}$ is the rescaled value,
$T_{C}^{\rm exp}{\rm(CrI_{3})}$ is the experimental one for the clean
monolayer (45K), and $T_{C}^{MFT}{\rm(CrI_{3})}$ represents the Curie
temperature calculated by mean field theory (MFT). Even after this rescaling,
the increase in the Curie temperature compared with that of pristine CrI3 is
noticeable.

Notice, however, that these calculations contain some approximations in
addition to the ones already inherent in MFT. The $J$ coupling constants are
assumed to be equal in the MFT approximation, but the computed ones are not
constant across the cell, with fluctuations of up to 100%. Furthermore, spin-
orbit coupling (SOC) has not been considered and is expected to increase the
Curie temperature due to anisotropy contributions. Therefore, we now turn to
another method. We map the DFT + SOC results for the fully ferromagnetic and
fully antiferromagentic ground state energies to those of an Ising
Hamiltonian:

$\mathcal{H_{I}}=-\sum_{i,j}{J_{ij}\bm{I}_{i}\bm{I}_{j}}.$ (12)

We obtain a value of $J=0.95$ meV for pristine CrI3, which translates into a
Curie temperature of 83 K. These numbers compare rather well with those
obtained with the $jx$ code and, as expected, are slightly higher. The results
for the H and O adatom cases are also presented in Table 2. We see that in
both cases the Curie temperature increases now. We attribute the discrepancy
in the case of O to the fact that the fluctuations of the $J$ coupling
constants are visibly larger than for H, and the MFT is probably less
reliable.

Finally, we fully functionalize the CrI3 surface with H and O atoms. Since the
favorable site for the H atom was found to be on top of the Cr atom, here all
the Cr atoms have H atoms on top. In the case of O, the atoms are placed on VI
sites, as this is the most stable configuration. The TDOS of the fully
functionalized HCrI3 and OCrI3 systems is plotted in Fig. 10-a along with that
of the clean CrI3 for comparative purposes. Here again, as we can see in this
figure, H leaves the TDOS of CrI3 almost without change, while O introduces
visible changes. The full functionalization of the CrI3 system with O atoms
has changed its electronic structure into that of a half-metal. This means
that the spin-up channel crosses the Fermi level, while we have a gap in the
opposite channel. Therefore, it acts as a conductor to spin-up electrons but
as a semiconductor to those of the spin-down electrons. The total
magnetization of the system decreases to 16$\mu_{B}$ in case of HCrI3 (i.e.
its total magnetism decreases by one per H atom) while it remains unchanged
for fully functionalized OCrI3 system.

The increase of the Curie temperature has been confirmed with two different
calculations for one H adatom, so we expect full coverage to have a stronger
impact. The Curie temperature can now be safely computed with the ”jx” code
since all Cr atoms are identical and, expectedly, we see a very significant
increase: $T_{C}^{\rm jx}=192^{\circ}K$ and $T_{C}\ (rescaled)=142^{\circ}K$.
We have ignored SOC again here, although, as we have seen, it may reflect in a
relatively lower increase of the Curie temperature. This calculation confirms
that hydrogenation can serve as a simple mechanism to enhance the Curie
temperature. In addition, the main properties of CrI3 (such as its
ferromagnetic feature and semiconductor band structure) remain unchanged.

Figure 10: TDOS of fully functionalized HCrI3 and OCrI3 systems. The positive
(negative) values correspond to spin-up (spin-down) states. Black lines show
The TDOS of the pure CrI3 (a), and fully functionalized CrI3 ($2\times 2$)
structures with Hydrogen and Oxygen atoms (b).

##  IV Conclusions

We have investigated, through DFT calculations, the electronic and magnetic
structures of the 2D monolayer FM insulator CrI3 resulting from H and O
adsorption. We have found that H adatoms preferentially adsorb on Cr atoms,
but hardly influence the original CrI3 structure although reduce the magnetic
moment of the hosting Cr atoms. Oxygen adatoms present a much higher
adsorption energy (also on Cr atoms) and greatly distort the atomic structure
while introducing an impurity band in the gap. Interestingly, we have
predicted that the intrinsic ferromagnetism and corresponding Curie
temperature can be largely enhanced by the adsorption of H atoms up to approx.
200 % for full H coverage. Since the main properties of the system remain
unchanged, hydrogenation may be a viable route towards room temperature
operation devices. Oxygen adsorption, on the other hand, offers the
possibility of turning CrI3 into a half metal.

##  V Acknowledgements

JJP acknowledges Spanish MINECO through Grant FIS2016-80434-P, the Fundación
Ramón Areces, the María de Maeztu Program for Units of Excellence in R&D
(MDM-2014-0377), the Comunidad Autoónoma de Madrid through S2018/NMT-4321
(NanomagCOST-CM) and the European Union Seventh Framework Programme under
Grant agreement No. 604391 Graphene Flagship. JJP also acknowledges the
computer resources and assistance provided by the Centro de Computación
Científica of the Universidad Autónoma de Madrid and the RES.

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, science 306, 666 (2004).
  * Dean  _et al._ (2010) C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, _et al._ , Nature nanotechnology 5, 722 (2010).
  * Cahangirov  _et al._ (2009) S. Cahangirov, M. Topsakal, E. Aktürk, H. Şahin, and S. Ciraci, Physical review letters 102, 236804 (2009).
  * Mak  _et al._ (2010) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Physical review letters 105, 136805 (2010).
  * Huang  _et al._ (2017) B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, _et al._ , Nature 546, 270 (2017).
  * Wang  _et al._ (2019) G. Wang, X. Wang, and Y. Liao, Applied Surface Science 471, 1011 (2019).
  * Fei  _et al._ (2018) Z. Fei, B. Huang, P. Malinowski, W. Wang, T. Song, J. Sanchez, W. Yao, D. Xiao, X. Zhu, A. F. May, _et al._ , Nature materials 17, 778 (2018).
  * Bonilla  _et al._ (2018) M. Bonilla, S. Kolekar, Y. Ma, H. C. Diaz, V. Kalappattil, R. Das, T. Eggers, H. R. Gutierrez, M.-H. Phan, and M. Batzill, Nature nanotechnology 13, 289 (2018).
  * Gong  _et al._ (2017) C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, _et al._ , Nature 546, 265 (2017).
  * Lee  _et al._ (2016) J.-U. Lee, S. Lee, J. H. Ryoo, S. Kang, T. Y. Kim, P. Kim, C.-H. Park, J.-G. Park, and H. Cheong, Nano letters 16, 7433 (2016).
  * Wang  _et al._ (2016a) X. Wang, K. Du, Y. Y. F. Liu, P. Hu, J. Zhang, Q. Zhang, M. H. S. Owen, X. Lu, C. K. Gan, P. Sengupta, _et al._ , 2D Materials 3, 031009 (2016a).
  * He  _et al._ (2019a) J. He, G. Ding, C. Zhong, S. Li, D. Li, and G. Zhang, Nanoscale 11, 356 (2019a).
  * McGuire  _et al._ (2015) M. A. McGuire, H. Dixit, V. R. Cooper, and B. C. Sales, Chemistry of Materials 27, 612 (2015).
  * Wang  _et al._ (2016b) H. Wang, F. Fan, S. Zhu, and H. Wu, EPL (Europhysics Letters) 114, 47001 (2016b).
  * Gao  _et al._ (2019) Y. Gao, J. Wang, Z. Li, J. Yang, M. Xia, X. Hao, Y. Xu, and F. Gao, physica status solidi (RRL)–Rapid Research Letters 13, 1800410 (2019).
  * Frisk  _et al._ (2018) A. Frisk, L. B. Duffy, S. Zhang, G. van der Laan, and T. Hesjedal, Materials Letters 232, 5 (2018).
  * Klein  _et al._ (2018) D. R. Klein, D. MacNeill, J. L. Lado, D. Soriano, E. Navarro-Moratalla, K. Watanabe, T. Taniguchi, S. Manni, P. Canfield, J. Fernández-Rossier, _et al._ , Science 360, 1218 (2018).
  * Jiang  _et al._ (2018) P. Jiang, L. Li, Z. Liao, Y. Zhao, and Z. Zhong, Nano letters 18, 3844 (2018).
  * Jiang  _et al._ (2019) W. Jiang, S. Li, H. Liu, G. Lu, F. Zheng, and P. Zhang, Physics Letters A 383, 754 (2019).
  * Zhang  _et al._ (2018) J. Zhang, B. Zhao, T. Zhou, Y. Xue, C. Ma, and Z. Yang, Physical Review B 97, 085401 (2018).
  * Song  _et al._ (2018) T. Song, X. Cai, M. W.-Y. Tu, X. Zhang, B. Huang, N. P. Wilson, K. L. Seyler, L. Zhu, T. Taniguchi, K. Watanabe, _et al._ , Science 360, 1214 (2018).
  * Mukherjee  _et al._ (2019) T. Mukherjee, S. Chowdhury, D. Jana, and L. L. Y. Voon, Journal of Physics: Condensed Matter 31, 335802 (2019).
  * Wu  _et al._ (2019) Z. Wu, J. Yu, and S. Yuan, Physical Chemistry Chemical Physics 21, 7750 (2019).
  * Zheng  _et al._ (2018) F. Zheng, J. Zhao, Z. Liu, M. Li, M. Zhou, S. Zhang, and P. Zhang, Nanoscale 10, 14298 (2018).
  * Webster and Yan (2018) L. Webster and J.-A. Yan, Physical Review B 98, 144411 (2018).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Physical review 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Physical review 140, A1133 (1965).
  * (28) “Openmx,” .
  * Bachelet (1982) G. Bachelet, Phys. Rev. B 26, 4199 (1982).
  * Troullier and Martins (1991) N. Troullier and J. L. Martins, Physical review B 43, 1993 (1991).
  * Vanderbilt (1990) D. Vanderbilt, Physical review B 41, 7892 (1990).
  * Blöchl (1990) P. E. Blöchl, Physical Review B 41, 5414 (1990).
  * Morrison  _et al._ (1993) I. Morrison, D. Bylander, and L. Kleinman, Physical Review B 47, 6728 (1993).
  * Ozaki (2003) T. Ozaki, Physical Review B 67, 155108 (2003).
  * Ozaki and Kino (2004) T. Ozaki and H. Kino, Physical Review B 69, 195113 (2004).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Ghosh  _et al._ (2019) S. Ghosh, N. Stojić, and N. Binggeli, Physica B: Condensed Matter (2019).
  * Skomski (2008) R. Skomski,  _Simple models of magnetism_ (Oxford University Press on Demand, 2008).
  * Guo  _et al._ (2018) Y. Guo, S. Yuan, B. Wang, L. Shi, and J. Wang, Journal of Materials Chemistry C 6, 5716 (2018).
  * Lado and Fernández-Rossier (2017) J. L. Lado and J. Fernández-Rossier, 2D Materials 4, 035002 (2017).
  * Zhang  _et al._ (2015) W.-B. Zhang, Q. Qu, P. Zhu, and C.-H. Lam, Journal of Materials Chemistry C 3, 12457 (2015).
  * Hu  _et al._ (2017) T. Hu, Z. Li, M. Hu, J. Wang, Q. Hu, Q. Li, and X. Wang, The Journal of Physical Chemistry C 121, 19254 (2017).
  * (43) “Visualization for electronic and structural analysis,” .
  * (44) “Atomsk,” .
  * He  _et al._ (2017) J. He, X. Li, P. Lyu, and P. Nachtigall, Nanoscale 9, 2246 (2017).
  * Mulliken (1955) R. S. Mulliken, The Journal of Chemical Physics 23, 1833 (1955).
  * Smart (1966) J. S. Smart,  _Effective field theories of magnetism_ (Saunders, 1966).
  * Reif (2009) F. Reif,  _Fundamentals of statistical and thermal physics_ (Waveland Press, 2009).
  * Han  _et al._ (2004) M. J. Han, T. Ozaki, and J. Yu, Physical Review B 70, 184421 (2004).
  * Terasawa  _et al._ (2019) A. Terasawa, M. Matsumoto, T. Ozaki, and Y. Gohda, Journal of the Physical Society of Japan 88, 114706 (2019).
  * Liu  _et al._ (2018) H. Liu, J.-T. Sun, M. Liu, and S. Meng, The journal of physical chemistry letters 9, 6709 (2018).
  * He  _et al._ (2019b) J. He, G. Ding, C. Zhong, S. Li, D. Li, and G. Zhang, Journal of Materials Chemistry C 7, 5084 (2019b).