# Low complexity convergence rate bounds for the synchronous gossip subclass
of push-sum algorithms

Balázs Gerencsér 33footnotemark: 3 B. Gerencsér is with the Alfréd Rényi
Institute of Mathematics, Budapest, Hungary and the Eötvös Loránd University,
Department of Probability and Statistics, Budapest, Hungary,
<EMAIL_ADDRESS>He is supported by the János Bolyai Research
Scholarship of the Hungarian Academy of Sciences. Miklós Kornyik M. Kornyik is
with the Alfréd Rényi Institute of Mathematics, Budapest, Hungary,
<EMAIL_ADDRESS>research was supported by NRDI (National Research,
Development and Innovation Office) grant KKP 137490.

###### Abstract

We develop easily accessible quantities for bounding the almost sure
exponential convergence rate of push-sum algorithms. We analyze the scenario
of i.i.d. synchronous gossip, every agent communicating towards its single
target at every step. Multiple bounding expressions are developed depending on
the generality of the setup, all functions of the spectrum of the network.
While the most general bound awaits further improvement, with more symmetries,
close bounds can be established, as demonstrated by numerical simulations.

##  Introduction

Average consensus algorithms have been around for a while [2], [18], with the
fundamental goal of computing the average of input values on a network in a
distributed manner with only local communication and simple operations. Often
some symmetry is imposed on the communication, in terms of the matrix
describing the linear update of the vector of values to be either doubly
stochastic, or even symmetric. This condition is quite well understood [17],
see the survey [16] also for applications, further discussion and references.

However, the interest for distributed averaging algorithms capable of handling
asynchronous directed communications emerged, naturally driving away the
representing update matrix from being doubly stochastic, still with the intent
to compute the exact average. As a result, the successful scheme of _push-sum_
was proposed [11], later also investigated under the name _ratio consensus_
[7] and joined by variants such as _weighted gossip_ [1]. The goal of these
algorithms is the same, but now using only local, directed communication and
without requiring message passing to happen synchronously or consistently
across the network. Given the simple objective of the algorithm, it also
serves as a building block for more complex tasks, e.g., the spectral analysis
of the network [12] or distributed optimization algorithms [13].

With other real-life communication challenges taken into account, the concept
has been extended in multiple ways to handle such aspects, including packet
loss [9] [15], delay [7] or even the presence of malicious agents [8]. In the
meantime, there is work to better understand the effect of such communication
deficiencies for the reference algorithms. The error of the consensus value
compared to the true average for the push-sum algorithm has been analyzed in
case of packet loss [6], similarly as it has been done for classic (linear)
gossip [3], [14].

An essential question in the analysis for usability and efficiency is
understanding the asymptotics of the processes, their convergence and the rate
at which it happens. In the cases above, the convergence of the push-sum
algorithm (or variants) has been confirmed. Additionally, for the original
push-sum scheme, an exponential convergence has been proven [11]. However, at
that time the focus was not yet on approximating the true rate. An important
step ahead was made in [10] providing a convincing upper bound along an
unspecified, infinite subset of the timeline for the almost sure (a.s.) rate
of convergence. More recently, the _exact rate of a.s. convergence_ has been
identified [5] for stationary ergodic updates as the spectral gap in terms of
the Lyapunov exponents of random matrix updates with generous applicability.
While being a clean representation with the concern that this Lyapunov
spectral gap is known to be _uncomputable_ in general [19]. As a follow-up, it
was possible to combine the inspiration of [10] and the tool-set of [5] to
obtain an actual upper bound on the a.s. rate for the i.i.d. case [4], now
formulated by manipulating the Kronecker square of a single (random) update
matrix, thus leading to a computable quantity. The bounds are solid, however
for a graph on $N$ vertices, matrices of $N^{2}\times N^{2}$ have to be
analyzed, quickly increasing in dimension.

Our goal is to provide even simpler convergence rate estimates. For this
purpose, we focus our attention to the natural setup, where a weighted network
determines the communication scheme driving the consensus process. In
particular, we assume synchronized gossip message passing, i.e. every node
sending a single packet to a single (random) recipient at each time slot.
Convergence of this scheme has been known since the formation of the push-sum
concept [11], ensuring that distributed average computation takes place.

The bounds provided can be computed directly once the standard spectral
description of the network is available. We are to formulate multiple
variants, both to provide general, but more conservative estimates, and also
sharper ones for a more restricted setting with stronger symmetries.

The rest of the paper is structured as follows. In the next section we
formally define the averaging process and state our results. Section 3 builds
a framework for the proof of the theorems, while Section 4 completes the
proofs. Detailed numerical performance analysis and concluding remarks are
provided in Section 5 and Section 6.

##  Main results

Let us introduce the push-sum algorithm along with other concepts that will be
used. Given is a finite graph $G=(V,E)$ with the vertex set
$V=[N]:=\\{1,2,\ldots,N\\}$, having degree sequence $d_{1},\ldots,d_{N}$.
There is an initial vector of values $x(0)\in\mathbb{R}^{N}$ at the vertices
to be averaged. The process is also using an auxiliary vector initialized at
$w(0)=\mathbf{1}\in\mathbb{R}^{N}$.

At each time step, a linear row-stochastic update - representing local
communication - is performed to both vectors as

$\displaystyle x(t)^{\top}$ $\displaystyle=x(t-1)^{\top}K(t),$ $\displaystyle
w(t)^{\top}$ $\displaystyle=w(t-1)^{\top}K(t).$

The average $\bar{x}:=\frac{1}{N}\sum_{i}x_{i}(0)$ is then locally estimated
by $x_{i}(t)/w_{i}(t)$. There is a wide generality of how $(K(t))_{t\geq 0}$
can be chosen. In the current paper, we focus on the scenario of i.i.d.
$K(t)$, when at each step, every vertex sends a single message to a randomly
chosen neighbor with a constant proportion and all these choices independent
from one another. Formally, $K(t)\stackrel{{\scriptstyle
d}}{{=}}K=(1-q)I+q\sum_{i}e_{i}e_{\beta_{i}}^{T}$ with some fixed $q\in[0,1]$
and independent $\beta_{i}$. By setting $p_{ij}:=\mathbb{P}(\beta_{i}=j)$, we
obtain an overall transition probability matrix $P$ which by construction has
to be compatible with the adjacency matrix of $G$.

For convenience, we introduce the notation $P_{q}=(1-q)I+qP$. It is easy to
check that $\mathbb{E}K=P_{q}$. In case $P$ has only real eigenvalues let
$\lambda_{i}$ denote its $i^{th}$ largest eigenvalue and let
$\lambda_{q,i}=(1-q)+q\lambda_{i}$ denote that of $P_{q}$. Following our setup
let us state our main results. Theorem 1 targets scenarios in more general
settings, while Theorem 2 is designed for more symmetric cases.

######  Theorem 1.

Let us consider a push-sum algorithm with message probability matrix $P$. Then

$\limsup\frac{1}{t}\max_{i}\log\left|\frac{x_{i}(t)}{w_{i}(t)}-\bar{x}\right|\leq\frac{1}{2}\log\rho\bigg{(}(I-J)(P_{q}^{\top}P_{q}+q^{2}(\Gamma-P^{\top}P))(I-J)\bigg{)}\quad
a.s.$ (1)

where $\Gamma$ is a diagonal matrix with $\gamma_{ii}=\sum_{j}p_{ji}$ and
$\rho(\cdot)$ denotes the spectral radius. Furthermore, if $P$ is symmetric,
then

$\limsup\frac{1}{t}\max_{i}\log\left|\frac{x_{i}(t)}{w_{i}(t)}-\bar{x}\right|\leq\frac{1}{2}\log((1-q)^{2}+2q(1-q)\lambda_{2}+q^{2})\quad
a.s.$ (2)

######  Remark 1.

In case each vertex chooses a recipient uniformly among its neighbors, the
bounding quantity in Theorem 1 will depend only on the graph structure,
furthermore the diagonal matrix $\Gamma$ takes the form
$\Gamma_{ii}=\sum\limits_{j:(j,i)\in E}d_{j}^{-1}$ with $d_{j}$ denoting the
degree of vertex $j$.

###### Proof.

It is easy to check that in this case $\mathbb{E}K=P_{q}$ with $P=D^{-1}A$,
where $D$ denotes the diagonal matrix consisting of the degrees of the
underlying graph’s vertices, while $A$ denotes the graph adjacency matrix.

∎

A better bound can be obtained for cases with stronger symmetries. A message
probability matrix is said to be transitive if for any pair $(i,j)$ there
exists a permutation matrix $\Pi$ with $\Pi_{ij}=1$ such that $\Pi
P\Pi^{-1}=P$.

######  Theorem 2.

Suppose that the message probability matrix $P$ is symmetric and transitive.
Then

$\limsup\frac{1}{t}\max_{i}\log\left|\frac{x_{i}(t)}{w_{i}(t)}-\bar{x}\right|\leq\frac{1}{2}\log\xi_{1},$
(3)

with $\xi_{1}$ being the largest root of the polynomial

$f(\xi)=\prod_{i>1}(\xi-\lambda_{q,i}^{2})-\frac{q^{2}}{N}\sum_{i>1}(1-\lambda_{i}^{2})\prod_{i\neq
j>1}(\xi-\lambda_{q,i}^{2}).$

######  Remark 2.

If $G$ is a transitive graph and each vertex chooses a recipient uniformly
among its neighbors, then the corresponding message probability matrix
satisfies the assumption of Theorem 2.

A special case has been analyzed in [10], where the underlying topology was
given by the complete graph with $q=1/2$. For this topology and general $q$
Theorem 2 immediately gives

$\displaystyle
f(\xi)=(\xi-(1-q)^{2})^{N-2}\big{(}\xi-(1-q)^{2}-q^{2}(1-N^{-1})\big{)}.$

from which the convergence rate bound $(1-q)^{2}+q^{2}(1-N^{-1})$ can be
easily obtained.

##  Tools

Let us first introduce a framework and corresponding tools in a general
setting. The alignment to the assumptions of the theorems will be carried out
later.

First we remark that the elements $w_{i}(t)$ are all positive because the
diagonal elements of the nonnegative $K(t)$ are strictly positive. Using the
notations $H(t)=K(1)K(2)\cdots K(t)$, $J=\mathbf{1}\mathbf{1}^{\top}/N$ easy
calculation shows

$\displaystyle x(t)^{\top}-\bar{x}w(t)^{\top}$
$\displaystyle=x_{0}^{\top}H(t)-\bar{x}w(t)^{\top}$
$\displaystyle=x_{0}^{\top}(JH(t)+(I-J)H(t))-\bar{x}w(t)^{\top}$
$\displaystyle=x_{0}^{\top}(I-J)H(t)$

meaning that

$\displaystyle\max_{i}\left|\frac{x_{i}(t)}{w_{i}(t)}-\bar{x}\right|$
$\displaystyle\leq
C\\{\min_{i}(w_{i}(t))\\}^{-1}||x_{0}||_{2}||(I-J)H(t)||_{2}$
$\displaystyle\leq
C\\{\min_{i}(w_{i}(t))\\}^{-1}||x_{0}||_{2}||(I-J)H(t)||_{F}.$ (4)

from some constant $C>0$. We are interested in the almost sure convergence
rate of the quantity on the left of (4). As we will see the dominant term will
be $||(I-J)H(t)||_{F}$. To get a handle on this factor let us analyze the
expectation of $||(I-J)H(t)||_{F}^{2}$.

$\displaystyle\mathbb{E}||(I-J)H(t)||_{F}^{2}$ $\displaystyle=\mathbb{E}\
\mathrm{Tr}\\{(I-J)H(t)H(t)^{\top}(I-J)\\}$
$\displaystyle=\mathrm{Tr}\\{(I-J)\mathbb{E}[H(t)H(t)^{\top}](I-J)\\}.$

According to the definition of $H(t)$ we can write

$\displaystyle\mathbb{E}[H(t)H(t)^{\top}]=\mathbb{E}\big{[}\mathbb{E}[K(1)\tilde{H}(t)\tilde{H}(t)^{\top}K(1)^{\top}\big{|}\tilde{H}(t)]\big{]}=\mathbb{E}[\mathbb{E}\big{.}[K(1)XK(1)^{\top}]\big{|}_{X=\tilde{H}(t)\tilde{H}(t)^{\top}}],$
(5)

where $\tilde{H}(t)=K(2)K(3)\cdots K(t)$, thus by the i.i.d. nature of the
updates $\tilde{H}(t)\stackrel{{\scriptstyle d}}{{=}}H(t-1)$. This motivates
the following definition of the linear operator $\Phi:\mathbb{R}^{N\times
N}\to\mathbb{R}^{N\times N}$ acting on matrices:

$\Phi(X):=\mathbb{E}[KXK^{\top}],$

which we will need to understand for further developing (5). For satisfactory
notation, before we progress let us introduce the linear operator
$\Psi:\mathbb{R}^{N\times N}\to\mathbb{R}^{N}$

$(\Psi(X))_{i}=x_{ii}$

and its pseudo-inverse $\Psi^{-}:\mathbb{R}^{N}\to\mathbb{R}^{N\times N}$

$(\Psi^{-}(v))_{ij}=\begin{cases}v_{i}&\mbox{ if }i=j,\\\ 0&\mbox{
otherwise.}\end{cases}$

Following the pattern of (5) we can prove the following

######  Lemma 1.

For any matrix $X$, we have

$\displaystyle\Phi(X)=P_{q}XP_{q}^{\top}+q^{2}\big{\\{}\Psi^{-}[P\
\Psi(X)]-\Psi^{-}\Psi(PXP^{\top})\big{\\}}$

reminding that $P_{q}=(1-q)I+qP$.

###### Proof.

Let $L:=\sum_{i}e_{i}e_{\beta_{i}}^{\top}$ then

$\displaystyle\mathbb{E}[KXK^{\top}]$
$\displaystyle=(1-q)^{2}X+q(1-q)\mathbb{E}[XL^{\top}]+q(1-q)\mathbb{E}[LX]+q^{2}\mathbb{E}[LXL^{\top}]$
$\displaystyle=(1-q)^{2}X+q(1-q)XP^{\top}+q(1-q)PX+q^{2}\mathbb{E}[LXL^{\top}]$
$\displaystyle=P_{q}XP_{q}^{\top}-q^{2}PXP^{\top}+q^{2}\mathbb{E}[LXL^{\top}].$

Next we will compute the term $\mathbb{E}[LXL^{\top}]$ as

$\displaystyle\mathbb{E}[LXL^{\top}]$
$\displaystyle=\sum_{i,i^{\prime}}\mathbb{E}[e_{i}e_{\beta_{i}}^{\top}Xe_{\beta_{i^{\prime}}}e_{i^{\prime}}^{\top}]=\sum_{i\neq
i^{\prime}}\mathbb{E}x_{\beta_{i},\beta_{i^{\prime}}}e_{i}e_{i^{\prime}}^{\top}+\sum_{i}\mathbb{E}x_{\beta_{i},\beta_{i}}e_{i}e_{i}^{\top}$
$\displaystyle=\sum_{\begin{subarray}{c}i\neq i^{\prime}\\\
j,j^{\prime}\end{subarray}}p_{ij}p_{i^{\prime}j^{\prime}}x_{jj^{\prime}}e_{i}e_{i^{\prime}}^{\top}+\sum_{i}p_{ij}x_{jj}e_{i}e_{i}^{\top}$
$\displaystyle=PXP^{\top}-\sum_{i,j,j^{\prime}}p_{ij}p_{ij^{\prime}}x_{jj^{\prime}}e_{i}e_{i}^{\top}+\sum_{i,j}p_{ij}x_{jj}e_{i}e_{i}^{\top}$
$\displaystyle=PXP^{\top}-\Psi^{-}\Psi(PXP^{\top})+\Psi^{-}(P\ \Psi(X)).$

Thus putting together the two parts gives

$\Phi(X)=P_{q}XP_{q}^{\top}+q^{2}\\{\Psi^{-}(P\Psi(X))-\Psi^{-}\Psi(PXP^{\top})\\}$

and this concludes the proof. ∎

In order to obtain a bound on $\mathrm{Tr}\hskip
2.27621pt\\{(I-J)\mathbb{E}[H(t)H(t)^{\top}](I-J)\\}$ it is enough to
understand $\Phi$, since

$\mathrm{Tr}\hskip
2.27621pt\\{(I-J)\mathbb{E}[H(t)H(t)^{\top}](I-J)\\}=\mathrm{Tr}\hskip
2.27621pt\\{(I-J)\Phi^{t}(I)(I-J)\\},$

where $\Phi^{t}(I)$ denotes the application of $\Phi$ on $I$ $t$ times, i.e.
$\underbrace{\Phi\circ\Phi\circ\cdots\circ\Phi}_{t\ times}(I)$.

######  Proposition 1.

The map $\Phi$ has the following fundamental properties:

  1. (P1)

$\Phi$ is linear,

  2. (P2)

$\Phi(X^{\top})=\Phi(X)^{\top}$,

  3. (P3)

for any skew-symmetric matrix $X$, $\Phi(X)=P_{q}XP_{q}^{\top},$

  4. (P4)

if $X\geq 0$ then $\Phi(X)\geq 0$, i.e. $\Phi$ keeps the positive semi-
definite property,

  5. (P5)

if $x_{kl}\geq 0$ $\forall(k,l)$, then $\Phi(X)_{kl}\geq 0$ $\forall(k,l)$,

  6. (P6)

$J$ is an eigenmatrix of $\Phi$ with eigenvalue $1$, i.e. $\Phi(J)=J$,

  7. (P7)

for $X\geq 0$, and $P=P^{\top}$ we have $\mathrm{Tr}\hskip
2.27621pt\Phi(X)\leq\mathrm{Tr}\hskip 2.27621ptX.$

For the adjoint map $\Phi^{*}$ the following observations can be added:

  1. (P*1)

$\Phi^{*}(Y)=P_{q}^{\top}YP_{q}+q^{2}\\{\Psi^{-}[P^{\top}\Psi(Y)]-P^{\top}(\Psi^{-}\Psi
Y)P\\}$,

  2. (P*2)

if $X\geq 0$ then $\Phi^{*}(X)\geq 0$, i.e. $\Phi^{*}$ also keeps the positive
semi-definite property,

  3. (P*3)

if $x_{kl}\geq 0$ $\forall(k,l)$, then $\Phi^{*}(X)\geq 0$ $\forall(k,l)$,

  4. (P*4)

if $X\mathbf{1}=0$ then $\Phi^{*}(X)\mathbf{1}=0$.

###### Proof.

The first three properties follow directly from the definition of $\Phi$,
hence their proofs are left to the respected reader.
Property (P4) can be proved as follows. Let $X\geq 0$ and let $w$ be an
arbitrary vector, then

$\displaystyle w^{\top}\Phi(X)w$
$\displaystyle=w^{\top}\mathbb{E}[KXK^{\top}]w=\mathbb{E}[w^{\top}KXK^{\top}w]\geq
0$

(P5) is analogous to (P4), namely

$\displaystyle\Phi(X)_{kl}=(\mathbb{E}[KXK^{\top}])_{kl}\geq 0.$

Property (P6) is the result of a short series of calculations.

$\displaystyle\Phi(J)$ $\displaystyle=PJP^{\top}+q^{2}(\Psi^{-}P\Psi
J-\Psi^{-}\Psi(PJP^{\top}))$ $\displaystyle=J+q^{2}(1/N\cdot I-1/N\cdot I)=J,$

due to the facts $PJ=JP^{\top}=J$ and $\Psi J=\mathbf{1}/N$.
Before proving (P7) let us note that due to $X\geq 0$ and the linearity of
$\Phi$ it is enough to prove this property for $X=xx^{\top}$. Using the
definition of $K=(1-q)I+q^{2}\sum_{i}e_{i}e_{\beta_{i}}^{\top}$ we have

$\displaystyle\mathrm{Tr}\hskip 2.27621pt\Phi(xx^{\top})$
$\displaystyle=\mathrm{Tr}\hskip
2.27621pt\mathbb{E}[Kxx^{\top}K^{\top}\textasciibreve]=$
$\displaystyle=\mathrm{Tr}\hskip
2.27621pt\mathbb{E}\\{(1-q)^{2}xx^{\top}+q(1-q)(Lxx^{\top}+xx^{\top}L^{\top})+q^{2}Lxx^{\top}L^{\top}\\}$
$\displaystyle=(1-q)^{2}||x||_{2}^{2}+2q(1-q)x^{\top}Px+q^{2}\mathbb{E}||Lx||_{2}^{2}$
$\displaystyle=(1-q)^{2}||x||_{2}^{2}+2q(1-q)x^{\top}Px+q^{2}\sum_{i,j}p_{ij}x_{j}^{2}$
$\displaystyle\leq||x||^{2}_{2}$

where in the last step we used the facts $P=P^{\top}$,
$P\mathbf{1}=\mathbf{1}$ and
$x^{\top}Px\leq\lambda_{1}(P)||x||_{2}^{2}=||x||_{2}^{2}$. Now we proceed with
proving the properties of the adjoint.
The proof of (P*1) is based on the following series of calculations: due to
the equivalences

$\displaystyle\mathrm{Tr}\hskip 2.27621pt(\Psi^{-}(P\ \Psi X))Y^{\top}$
$\displaystyle=\sum_{i,k}p_{ik}x_{kk}y_{ii}=\sum_{k}x_{kk}\sum_{i}p_{ik}y_{ii}=\mathrm{Tr}\hskip
2.27621pt\\{X\Psi^{-}(P^{\top}Y^{\top})\\},$ $\displaystyle\mathrm{Tr}\hskip
2.27621pt\Psi^{-}\Psi(PXP^{\top})Y^{\top}$
$\displaystyle=\sum_{i,k,l}p_{ik}x_{kl}p_{il}y_{ii}=\sum_{k,l}x_{kl}\sum_{i}p_{ik}y_{ii}p_{il}=\mathrm{Tr}\hskip
2.27621pt\\{XP^{\top}\Psi(Y^{\top})P\\}$

we have

$\displaystyle\langle\Phi(X),Y\rangle$ $\displaystyle=\mathrm{Tr}\hskip
2.27621pt\Phi(X)Y^{\top}=\mathrm{Tr}\hskip
2.27621ptP_{q}XP_{q}^{\top}Y^{\top}+q^{2}\mathrm{Tr}\hskip
2.27621pt\bigg{\\{}[\Psi^{-}(P\Psi(X))]Y^{\top}-[\Psi^{-}\Psi(PXP^{\top})]Y^{\top}\bigg{\\}}$
$\displaystyle=\mathrm{Tr}\hskip
2.27621pt(XP_{q}^{\top}Y^{\top}P_{q})+q^{2}\mathrm{Tr}\hskip
2.27621pt\\{X\Psi^{-}(P^{\top}\Psi Y^{\top})-XP^{\top}\Psi(Y^{\top})P\\}$
$\displaystyle=\langle X,\Phi^{*}(Y)\rangle.$

Properties (P*2) and (P*3) can be confirmed analogously to (P4), (P5).
(P*4) is a result of the short derivation

$\displaystyle\Phi^{*}(X)\mathbf{1}$
$\displaystyle=P_{q}^{\top}XP_{q}\mathbf{1}+q^{2}(\Psi^{-}P\Psi
X\mathbf{1}-P^{\top}(\Psi^{-}\Psi X)P\mathbf{1})$
$\displaystyle=0+q^{2}(\Psi^{-}(P^{\top}\Psi(X))\mathbf{1}-P^{\top}(\Psi^{-}\Psi
X)\mathbf{1})$

since $P_{q}\mathbf{1}=P\mathbf{1}=\mathbf{1}$ and we assumed $X\mathbf{1}=0$.
For the second term we have

$\displaystyle(\Psi^{-}(P^{\top}\Psi(X))\mathbf{1})_{i}$
$\displaystyle=\sum_{j}p_{ji}x_{jj}$ $\displaystyle(P^{\top}(\Psi^{-}\Psi
X)\mathbf{1})_{i}$ $\displaystyle=\sum_{j}p_{ji}x_{jj},$

so $\Phi^{*}(X)\mathbf{1}=0$. This concludes the proof. ∎

######  Remark 3.

Properties (P1), (P2), (P4), (P7) mean that $\Phi$ describes a quantum
operation.

According to the previously listed properties we have

######  Corollary 1.

The cone of positive semi-definite matrices is invariant under the action of
$\Phi$.

The following statement is going to help us in proving Theorem 1 as our focus
is on the speed of convergence and not the limit $\bar{x}\mathbf{1}$, which is
of constant order.

######  Corollary 2.

Since $\mathrm{Tr}\hskip 2.27621pt(I-J)H(t)H(t)^{\top}(I-J)=\mathrm{Tr}\hskip
2.27621ptH(t)^{\top}(I-J)H(t)$ and the adjoint of the linear map $f:X\mapsto
AXA^{\top}$ is the map $f^{*}:X\mapsto A^{\top}XA$, it is easy to show that
$\Phi^{*}(X)=\mathbb{E}[K^{\top}XK]$ for any symmetric matrix $X$.

######  Remark 4.

According to the definition of $K$, we have $KJ=J$, and so

$(I-J)K(I-J)=(I-J)(K-J)=K-JK=(I-J)K$

whence

$\mathbb{E}[(I-J)KXK^{T}(I-J)]=\mathbb{E}[(I-J)K(I-J)X(I-J)K^{T}(I-J)].$

######  Corollary 3.

Let us define the operator $\widehat{\Phi}$ as

$\widehat{\Phi}:X\mapsto(I-J)\Phi(X)(I-J)\in\mathrm{End}(\\{Y\in\mathbb{R}^{N\times
N}:Y=Y^{\top},YJ=0\\}).$

According to properties of $\Phi$ combined with Remark 4, we have

$\widehat{\Phi}^{t}(X)=(I-J)\Phi^{t}(X)(I-J),\quad
X\in\\{Y\in\mathbb{R}^{N\times N}:Y=Y^{\top},YJ=0\\}$

furthermore

$\Phi((I-J)X(I-J))=\Phi(X)+JXP_{q}^{T}+P_{q}XJ-JXJ.$ (6)

##  Proofs

The next proposition is the final step before proving Theorem 1.

######  Proposition 2.

Let $P$ be a row stochastic matrix. Then

$\mathrm{Tr}\hskip 2.27621pt(\Phi^{*})^{t}(I-J)\leq
N\rho\big{(}(I-J)B_{q}(I-J)\big{)}^{t}$

where $B_{q}=P_{q}^{\top}P_{q}+q^{2}(\Gamma-P^{\top}P)$ is a positive definite
matrix, and $D$ is a diagonal matrix with $d_{ii}=\sum_{j}p_{ji}$ on its
diagonal.

###### Proof.

Due to its properties
$\Phi^{*}:\mathcal{X}_{0}\to\mathcal{X}_{0}=\\{X\in\mathbb{R}^{N\times
N}:X=X^{\top},XJ=0\\},$ and
$X\geq 0\implies\Phi^{*}(X)\geq 0$, we have

$\displaystyle\mathrm{Tr}\hskip 2.27621pt(\Phi^{*})^{t}(I-J)$
$\displaystyle\leq
N\rho((\Phi^{*})^{t}(I-J))=N\max\\{v^{\top}(\Phi^{*})^{t}(I-J)v:\|v\|=1,v\perp\mathbf{1}\\}$
$\displaystyle=N\|(\Phi^{*})^{t}(I-J)\|_{2}\leq
N\|(\Phi^{*})^{t}\|_{\mathcal{X}_{0}\to\mathcal{X}_{0}}\|(I-J)\|_{2}$
$\displaystyle\leq N\|\Phi^{*}\|^{t}_{\mathcal{X}_{0}\to\mathcal{X}_{0}}$

where

$\|\Phi^{*}\|_{\mathcal{X}_{0}\to\mathcal{X}_{0}}=\max\\{\|\Phi^{*}(X)\|_{2}:\|X\|_{2}\leq
1,X\in\mathcal{X}_{0}\\}.$

It is not hard to show that $\arg\max\|\Phi^{*}(X)\|_{2}\geq 0$, since let
$X=\arg\max\|\Phi^{*}(X)\|_{2}$ and let us consider its decomposition
$X=X^{*}-X^{-}$ with $X^{+},X^{-}\geq 0$. Then

$v^{\top}\Phi^{*}(X)v=v^{\top}\Phi^{*}(X^{+})v-v^{\top}\Phi^{*}(X^{-})v\leq
v^{\top}\Phi^{*}(X^{+})v+v^{\top}\Phi^{*}(X^{-})v,$

and this would lead to a contradiction if $X^{-}$ was not $0$. This implies

$\|\Phi^{*}\|_{\mathcal{X}_{0}\to\mathcal{X}_{0}}=\max\\{v^{\top}\Phi^{*}(X)v:X\in\mathcal{X}_{0},X\geq
0,\|X\|_{2}\leq 1,\|v\|\leq 1,v\perp\mathbf{1}\\}.$

meaning that it is enough to bound $v^{\top}\Phi^{*}(X)v$ from above. Let
$v\perp\mathbf{1}$ and $X\in\mathcal{X}_{0}$ with $\|X\|_{2}\leq 1$, then

$\displaystyle v^{\top}\Phi^{*}(X)v$
$\displaystyle=(P_{q}v)^{\top}XP_{q}v+q^{2}(\sum_{i,j}p_{ij}x_{ii}v_{j}^{2}-\sum_{i}(Pv)_{i}^{2}x_{ii}).$

Due to the conditions imposed on $X$, we have $x_{ii}\in[0,1]$, furthermore
$\sum_{j}p_{ij}-(Pv)_{i}^{2}\geq 0$ for any $i$, thus

$\displaystyle v^{\top}\Phi^{*}(X)v$
$\displaystyle\leq\|P_{q}v\|^{2}+q^{2}(\sum_{i,j}p_{ij}v_{j}^{2}-\sum_{i}(Pv)_{i}^{2})$
$\displaystyle=v^{\top}\big{(}P_{q}^{\top}P_{q}+q^{2}(\Gamma-P^{\top}P)\big{)}v\leq\rho\bigg{(}(I-J)(P^{\top}_{q}P_{q}+q^{2}(\Gamma-P^{\top}P)(I-J)\bigg{)}$

where $\Gamma$ is a diagonal matrix with diagonal elements
$\gamma_{ii}=\sum_{j}p_{ji}$. Note that for symmetric $P$ we have $\Gamma=I$
and so the upper bound above becomes $(1-q)^{2}+2q(1-q)\lambda_{2}+q^{2}.$ ∎

With all the tools at our hands we can prove Theorem 1.

###### Proof of Theorem 1.

According to Lemma 10 in [4] whose assumptions are clearly satisfied we have

$\limsup_{t}\frac{1}{t}\log\frac{1}{\min_{i}w_{i}(t)}\leq 0$

thus considering the quantity in (4) we can infer that

$\limsup_{t}\frac{1}{t}\log\big{(}\\{\min_{i}w_{i}(t)\\}^{-1}||x_{0}||\cdot||(I-J)H(t)||_{F}\big{)}\leq\limsup_{t}\frac{1}{t}\log||(I-J)H(t)||_{F}.$

The matrix $(I-J)H(t)$ can we written as a product of the i.i.d. random
matrices $(I-J)K(t)$, therefore due to the Fürstenberg-Kesten theorem it
follows that

$\displaystyle\limsup_{t}\frac{1}{t}\log||(I-J)H(t)||_{F}$
$\displaystyle=\lim_{t}\frac{1}{t}\log||(I-J)H(t)||_{F}=\lim_{t}\frac{1}{t}\mathbb{E}\log||(I-J)H(t)||_{F}$
$\displaystyle=\lim_{t}\frac{1}{2t}\mathbb{E}\log\|(I-J)H(t)\|_{F}^{2}.$

Using Jensen’s inequality yields

$\mathbb{E}\log||(I-J)H(t)||^{2}_{F}\leq\log\mathbb{E}||(I-J)H(t)||_{F}^{2}$

thus

$\limsup_{t}\frac{1}{t}\log\|(I-J)H(t)\|_{F}\leq\lim_{t}\frac{1}{2t}\log\mathbb{E}\|(I-J)H(t)\|_{F}^{2}.$

By taking the expectation we obtain

$\mathbb{E}||(I-J)H(t)||_{F}^{2}=\mathrm{Tr}\hskip
2.27621pt(\Phi^{*})^{t}(I-J).$

Combining the series of calculations above with Proposition 2 we arrive at

$\limsup_{t}\frac{1}{2t}\log\mathrm{Tr}\hskip
2.27621pt(\Phi^{*})^{t}(I-J)\leq\frac{1}{2}\begin{cases}\rho\bigg{(}(I-J)(P_{q}^{\top}P_{q}+q^{2}(\Gamma-P^{\top}P))(I-J)\bigg{)}&\mbox{in
general,}\\\ (1-q)^{2}+2q(1-q)\lambda_{2}+q^{2}&\mbox{ if
}P=P^{\top}.\end{cases}$

∎

Now we will turn to the case when the underlying graph is transitive. In this
scenario it is possible to give stronger bounds, but in order to do this we
need to reformulate the problem. Rearranging our main quantity of interest as

$\mathbb{E}||(I-J)H(t)||_{F}^{2}=\mathrm{Tr}\hskip
2.27621pt\mathbb{E}[H(t)^{\top}(I-J)H(t)],$

therefore

$\mathbb{E}||(I-J)H(t)||_{F}^{2}=\mathrm{Tr}\hskip
2.27621pt\\{(\Phi^{*})^{t}(I-J)\\}.$

Let us define $X_{t}=(\Phi^{*})^{t}(I-J).$ The following two lemmas will help
us in our progress.

######  Lemma 2.

Assume that $P$ is a kernel of the transitive Markov chain, implying that it
is symmetric and any diagonal element $p^{(k)}_{ii}$ of $P^{k}$ depends solely
on $k$ and not on $i$. Then we have

  1. 1.

$X_{t}\in\mathcal{P}_{t}:=\mathrm{Span}\\{P^{k},J\ ;\ 0\leq k\leq 2t\\}$,
hence the diagonal of $X_{t}$ is also constant,

  2. 2.

$X_{t+1}=P_{q}X_{t}P_{q}+q^{2}r_{t}(I-P^{2})$, where $r_{t}$ denotes the
common diagonal element of $X_{t}$.

###### Proof.

We prove by induction. For $t=0$ $X_{0}=I-J$, which trivially is a polynomial
of $P$ and $J$. For the induction step $t\to t+1$ assume that
$X_{t}\in\mathcal{P}_{t}$, then

$\displaystyle X_{t+1}$
$\displaystyle=\Phi^{*}(X_{t})=P_{q}X_{t}P_{q}+q^{2}\\{(\Psi^{-}P\Psi
X_{t})-P(\Psi^{-}\Psi X_{t})P\\}$

since $X_{t}$ is a polynomial of $P$ and $J$ it is also transitive. Noting
$\Psi X_{t}=r_{t}\mathbf{1}$ and $\Psi^{-}\Psi X_{t}=r_{t}I$, where $r_{t}$
denotes the common diagonal element of $X_{t}$, we can derive the recursion

$\displaystyle X_{t+1}=P_{q}X_{t}P_{q}+q^{2}r_{t}(I-P^{2})$ (7)

showing that $X_{t+1}\in\mathcal{P}_{t+1}$. ∎

######  Proposition 3.

If $P$ is symmetric and transitive then $X_{t}$ and $P$ possess the same
eigenvectors. If $v$ is an eigenvector to $P$ and $X_{t}$ corresponding to the
eigenvalue $\lambda$ and $\mu_{t}$ respectively then the following recursion
holds for $\mu_{t}$:

$\mu_{t+1}=\lambda_{q}^{2}\mu_{t}+q^{2}r_{t}(1-\lambda^{2}).$ (8)

Recall that $r_{t}$ is the common diagonal element of $X_{t}$. Furthermore the
largest eigenvalue of $X_{t}$ is asymptotically bounded by the second largest
root $\xi_{2}$ of the polynomial

$p(x)=\prod_{i}(x-\lambda_{q,i}^{2})\bigg{(}1+\frac{q^{2}}{N}\sum_{j>1}\frac{1-\lambda_{j}^{2}}{x-\lambda_{q,j}^{2}}\bigg{)}$

in the following sense:

$\limsup_{t}\frac{1}{t}\log\max_{i}\mu_{t,i}\leq\log\xi_{2}.$

###### Proof.

Using $r_{t}=\frac{1}{N}\mathrm{Tr}\hskip
2.27621ptX_{t}=\frac{1}{N}\sum_{i}\mu_{t,i}$ we can write the recursion
described in (8) as

$\boldsymbol{y}_{t+1}=\bigg{(}D\mathbf{+}\frac{q^{2}}{N}\mathbf{b}\mathbf{1}^{\top}\bigg{)}\boldsymbol{y}_{t}$
(9)

with the vectors $(\boldsymbol{y}_{t})_{i}=\mu_{t,i}$,
$\mathbf{b}_{i}=1-\lambda_{i}^{2}$, and $D=\mathrm{diag}\hskip
2.84526pt(\lambda_{q,1}^{2},\ldots,\lambda_{q,N}^{2})$, where $\lambda_{q,i}$
denotes the $i^{th}$ largest eigenvalue of $P_{q}$ and $\mu_{t,i}$ denotes the
eigenvalue of $X_{t}$ corresponding to $\lambda_{q,i}$. By correspondence we
mean in the sense of defined by the recursion (8), in which case
$\mu_{t,1}=0$, since $\lambda_{q,1}=1$ and $X_{t}\mathbf{1}=0$ for any $t$.
Let $(\lambda,v)$ denote an eigen-pair of $P$. According to the previous lemma
$X_{t}$ is a polynomial of $P$ and $J$, hence $v$ is an eigenvector of
$X_{t}$, furthermore, due to the symmetry of $P$, we have $JP=PJ=J$. Recursion
(7) then yields

$\displaystyle X_{t+1}v$ $\displaystyle=P_{q}X_{t}P_{q}v+q^{2}r_{t}(I-P^{2})v$
$\displaystyle=\lambda_{q}^{2}\mu_{t}v+q^{2}r_{t}(1-\lambda^{2})v$

and we have

$\mu_{t+1}=\lambda_{q}^{2}\mu_{t}+q^{2}r_{t}(1-\lambda^{2}),$

proving the first part.

Before continuing with the proof of the second part, let us note that $e_{1}$
is a left-eigenvector of the matrix $D+q^{2}\mathbf{b}\mathbf{1}^{\top}/N$
corresponding to the eigenvalue 1, meaning that the right-eigenvectors
corresponding to a different eigenvalue are orthogonal to $e_{1}$. Now let us
choose the following vectors as basis: $f_{1}=e_{1}$, $f_{i}=e_{i}-e_{1},i>1$,
then for $i>1$ we have $f_{i}\perp\mathbf{1}$. Writing equation (9) in basis
$\\{f_{i}\\}$ yields

$\displaystyle Df_{1}$ $\displaystyle=De_{1}=\lambda_{q,1}^{2}e_{1},$
$\displaystyle Df_{i}$
$\displaystyle=De_{i}-De_{1}=\lambda_{q,i}^{2}e_{i}-\lambda_{q,1}^{2}e_{1}=$
$\displaystyle=\lambda_{q,1}^{2}f_{i}+(\lambda_{q,i}^{2}-\lambda_{q,1}^{2})f_{1},\quad
i>1,$ $\displaystyle\mathbf{b}\mathbf{1}^{\top}f_{1}$
$\displaystyle=\mathbf{b}=b_{1}f_{1}+\sum_{i>1}b_{i}(f_{i}+f_{1}),$
$\displaystyle\mathbf{b}\mathbf{1}^{\top}f_{i}$ $\displaystyle=0\quad i>1,$

thus rewriting the matrix $D+q^{2}\mathbf{b}\mathbf{1}^{\top}/N$ in this new
basis gives us

$\begin{bmatrix}(q^{2}/N)\sum_{i}b_{i}+\lambda_{q,1}^{2}&-(\lambda_{q,1}^{2}-\lambda_{q,2}^{2})&-(\lambda_{q,1}^{2}-\lambda_{q,3}^{2})&\ldots&-(\lambda_{q,1}^{2}-\lambda_{q,N}^{2})\\\
(q^{2}/N)b_{2}&\lambda_{q,2}^{2}&0&\ldots&0\\\
(q^{2}/N)b_{3}&0&\lambda_{q,3}^{2}&\ldots&0\\\ \vdots&&&\ddots&\\\
(q^{2}/N)b_{N}&0&\ldots&0&\lambda_{q,N}^{2}\end{bmatrix}$ (10)

The characteristic polynomial of the matrix in (10) can be computed via
expanding along the first column, leading to

$\displaystyle p(x)$
$\displaystyle=\bigg{(}x-\frac{q^{2}}{N}\sum_{i}b_{i}-\lambda_{q,1}^{2}\bigg{)}\prod_{i>1}(x-\lambda_{q,i}^{2})+\sum_{i>1}\frac{q^{2}}{N}b_{i}(\lambda_{q,1}^{2}-\lambda_{q,i}^{2})\prod_{1<j\neq
i}(x-\lambda_{q,j}^{2})$
$\displaystyle=\prod_{i>1}(x-\lambda_{q,i}^{2})\bigg{\\{}x-\lambda_{q,1}^{2}-\frac{q^{2}}{N}\sum_{i}(1-\lambda_{i}^{2})+\frac{q^{2}}{N}\sum_{j>1}(1-\lambda_{j}^{2})\frac{(1-\lambda_{q,j}^{2})}{x-\lambda_{q,j}^{2}}\bigg{\\}}.$

Exploiting the fact $\lambda_{1}=\lambda_{q,1}=1$ we obtain

$\displaystyle p(x)$
$\displaystyle=\prod_{i>1}(x-\lambda_{q,i}^{2})\bigg{\\{}x-1+\frac{q^{2}}{N}\sum_{j>1}\bigg{[}(1-\lambda_{j}^{2})\bigg{(}\frac{1-\lambda_{q,j}^{2}}{x-\lambda_{q,j}^{2}}-1\bigg{)}\bigg{]}\bigg{\\}}$
$\displaystyle=\prod_{i>1}(x-\lambda_{q,i}^{2})\bigg{\\{}x-1+\frac{q^{2}}{N}\sum_{j>1}(1-\lambda_{j}^{2})\frac{1-x}{x-\lambda_{q,j}^{2}}\bigg{\\}}$
$\displaystyle=\prod_{i}(x-\lambda_{q,i}^{2})\bigg{(}1-\frac{q^{2}}{N}\sum_{j>1}\frac{1-\lambda_{j}^{2}}{x-\lambda_{q,j}^{2}}\bigg{)}.$

We note here that due the initialization $\mu_{0}=(0,1,\ldots,1)$ and the fact
that $e_{1}$ is a left-eigenvector of the recursion, we have for each
subsequent vector $\mu_{t}\perp e_{1}$. This means that we are only interested
in the second largest root $\eta_{2}$ of the characteristic polynomial to
obtain the asymptotic growth rate of $\mu_{t}$, i.e.

$||\mu_{t}||_{\infty}\leq C\eta_{2}^{t}\ ||\mu_{0}||_{\infty}$

for some $C>0$. ∎

###### Proof of Theorem 2.

The first part of the proof, up until the point where we have to bound the
quantity $\mathbb{E}||(I-J)H(t)||_{F}^{2}$ from above, is analogous to the
proof of Theorem 1, hence it will be omitted here. Before proceeding with the
actual proof we note that

$||(I-J)H(t)||_{F}^{2}=\mathrm{Tr}\hskip
2.27621pt(I-J)H(t)H(t)^{\top}(I-J)=\mathrm{Tr}\hskip
2.27621ptH(t)^{\top}(I-J)H(t),$

whence

$\mathbb{E}||(I-J)H(t)||_{F}^{2}=\mathrm{Tr}\hskip
2.27621pt\\{(\Phi^{*})^{t}(I-J)\\}=\sum_{j}\mu_{t,j}\leq N\max_{j}\mu_{t,j},$

where $\mu_{t,j}$ denotes the eigenvalue of $X_{t}=(\Phi^{*})^{t}(I-J)$
corresponding to the $j^{th}$ largest eigenvalue $\lambda_{j}$ of $P$. Due to
$X_{0}=I-J$ and $\lambda_{1}=1$ with $v_{1}=c\mathbf{1}$, we have
$\mu_{0,1}=0$ and due to the recursion (8), $\mu_{t,1}=0$ for any $t>0$. In
Proposition 3 it was shown that

$||\mu_{t}||_{\infty}\leq C\eta_{2}^{t}||\mu_{0}||=C\eta_{2}^{t},$

where $\eta_{2}$ denoted the second largest root of the polynomial

$\displaystyle
p(x)=\prod_{i}(x-\lambda_{q,i}^{2})-\frac{q^{2}}{N}\sum_{i>1}(1-\lambda_{i}^{2})\prod_{i\neq
j>1}(x-\lambda_{q,i}^{2}).$

Exploiting the fact that $1$ is a root of $p$ we have

$\displaystyle
p_{1}(x):=p(x)/(x-1)=\prod_{i>1}(x-\lambda_{q,i}^{2})-\frac{q^{2}}{N}\sum_{i>1}(1-\lambda_{i}^{2})\prod_{i\neq
j>1}(x-\lambda_{q,i}^{2})$

thus by denoting $\xi_{1}$ the largest root of $p_{1}(x)$ we can write
$||\mu_{t}||_{\infty}\leq C\xi_{1}^{t}$. Altogether

$\limsup_{t}\frac{1}{2t}\log\mathbb{E}||(I-J)H(t)||_{F}^{2}\leq\frac{1}{2}\log\xi_{1},$

and this concludes the proof. ∎

##  Numerical experiments

In order to evaluate the performance of the bounds obtained on the convergence
rate, we perform a detailed numerical comparison of the available quantities.
For setups with $N\leq 120$ the synchronous gossip process is realized for
$t=500$ steps and the approximate rate is expressed as

$\frac{1}{t}\log\bigg{\|}\frac{1}{\sqrt{N}}(I-J)H(t)\Psi^{-}(w(t))^{-1}\bigg{\|}_{F}.$

For $N>120$ we set $t=1000$ and, for complexity and memory usage reduction, we
use the modified approximation

$\frac{1}{t}\log\bigg{\|}\frac{1}{\sqrt{M}}X^{M}(0)H(t)\Psi^{-}(w(t))^{-1}\bigg{\|}_{F}$

where $X^{M}(0)$ is an $M\times N$ matrix of uniform random independent rows
in $1^{\perp}$ with unit norm, representing various initializations, reaching
the principal rate with high probability. We choose
$M=\lfloor\sqrt{N}\rfloor$.

This is compared both with $\eta:=\frac{1}{2}\log\rho((I-J)^{\otimes
2}\mathbb{E}(K(1)^{\otimes 2}))$ which is the bound of [4] and the bounds
developed in the current paper.

For the context of Theorem 1, more precisely its specialized version Corollary
1 we build a Barabási-Albert random graph with 2 edges added with each new
vertex, and assign uniform probabilities for choosing gossip recipients. See
Figure 1 for the resulting approximate rates and bounds as $q$ varies in
$(0,1)$. We can see that $\eta$ provides a very close fit, the current bound
is only moderately usable for small $q$. However, note that for $N=100$
handling the $N^{2}\times N^{2}$ matrices needed for the tensor product based
bound is getting computationally heavy, thus we only plot the simulations
versus the current bound. We also conducted simulation experiments for
$N=5000$ in which case our bound was mostly nonnegative and hence did not
carry any useful information.

(a) $N=24$ and $k=2$.

(b) $N=100$ and $k=2$

Figure 1: Simulations and bounds (1) for general graph without symmetries,
generated by the Barabási-Albert algorithm.

In the case of symmetric $P$, the second part of Theorem 1 provides a more
usable bound as shown in Figure 2. We use random regular graphs and uniform
recipient probabilities again. The resulting bound captures the linear trend
when $q$ is near 0, it is a factor $\approx 2$ off from the numerical value
when $q\leq 0.5$ and then deteriorates.

(a) $N=24$ and $d=4$.

(b) $N=100$ and $d=5$.

(c) $N=5000$ and $d=6$.

Figure 2: Simulations and bounds (2) for random regular graphs.

For the more refined bound of Theorem 2, we generate a random Cayley-graph of
$S_{k}$. This time we see a substantially better fit as shown in Figure 3. In
fact, it recovers $\eta$ which is expected from the exact analysis carried out
in its proof, but circumvents the necessity of working with the large matrices
of the Kronecker-product.

(a) $S_{4}$ with 2 generators
$N=24,d=4$

(b) $S_{5}$ with 3 generators
$N=120,d=5$

(c) $S_{7}$ with 3 generators,
$N=5040$, $d=6$

Figure 3: Simulations and bounds (3) for random Cayley-graphs.

##  Summary and future plans

In this work we have presented bounds of various accuracy depending on the
level of symmetry of the underlying topology. Along the way proving our main
results we have developed a framework as described in Section 3 relying on
matrix operators $\Phi,\Phi^{*}$ that we hope can be useful when analyzing
similar dynamics.

The computational cost of these bounds is orders of magnitude less than that
of the simulations or computation of $\eta$ from [4] confirming their
usefulness in assessing the efficiency of various push-sum algorithms.

Concerning our future plans, we observed in some of our numerical experiments
that the expression in (3) was also a valid bound for regular graphs. This
lead us to state the following

######  Conjecture 1.

The bound given in 2 remains valid in the simple symmetric case, e.g. for
regular graphs without transitivity.

## References

  * [1] F. Bénézit, V. Blondel, P. Thiran, J. N. Tsitsiklis, and M. Vetterli, Weighted gossip: Distributed averaging using non-doubly stochastic matrices, in Proceedings of 2010 IEEE International Symposium on Information Theory (ISIT), 2010, pp. 1753–1757.
  * [2] V. Blondel, J. M. Hendrickx, A. Olshevsky, and J. N. Tsitsiklis, Convergence in multiagent coordination, consensus, and flocking, in 44th IEEE Conference on Decision and Control, 2005, pp. 2996–3000.
  * [3] P. Frasca and J. M. Hendrickx, Large network consensus is robust to packet losses and interferences, in 2013 European Control Conference (ECC), IEEE, 2013, pp. 1782–1787.
  * [4] B. Gerencsér, Computable convergence rate bound for ratio consensus algorithms, IEEE Control Systems Letters, 6 (2022), pp. 3307–3312.
  * [5] B. Gerencsér and L. Gerencsér, Tight bounds on the convergence rate of generalized ratio consensus algorithms, IEEE Transactions on Automatic Control, 67 (2022), pp. 1669–1684.
  * [6] B. Gerencsér and J. M. Hendrickx, Push sum with transmission failures, IEEE Transactions on Automatic Control, 64 (2018), pp. 1019–1033.
  * [7] C. N. Hadjicostis and T. Charalambous, Average consensus in the presence of delays in directed graph topologies, IEEE Transactions on Automatic Control, 59 (2014), pp. 763–768.
  * [8] C. N. Hadjicostis and A. D. Domínguez-García, Trustworthy distributed average consensus, in 2022 IEEE 61st Conference on Decision and Control (CDC), IEEE, 2022, pp. 7403–7408.
  * [9] C. N. Hadjicostis, N. H. Vaidya, and A. D. Domínguez-García, Robust distributed average consensus via exchange of running sums, IEEE Transactions on Automatic Control, 61 (2015), pp. 1492–1507.
  * [10] F. Iutzeler, P. Ciblat, and W. Hachem, Analysis of sum-weight-like algorithms for averaging in wireless sensor networks, IEEE Transactions on Signal Processing, 61 (2013), pp. 2802–2814.
  * [11] D. Kempe, A. Dobra, and J. Gehrke, Gossip-based computation of aggregate information, in Proceedings of 44th Annual IEEE Symposium on Foundations of Computer Science, 2003, pp. 482–491.
  * [12] D. Kempe and F. McSherry, A decentralized algorithm for spectral analysis, Journal of Computer and System Sciences, 74 (2008), pp. 70–83.
  * [13] A. Nedić and A. Olshevsky, Distributed optimization over time-varying directed graphs, IEEE Transactions on Automatic Control, 60 (2014), pp. 601–615.
  * [14] R. Olfati-Saber and R. M. Murray, Consensus problems in networks of agents with switching topology and time-delays, IEEE Transactions on Automatic Control, 49 (2004), pp. 1520–1533.
  * [15] A. Olshevsky, I. C. Paschalidis, and A. Spiridonoff, Fully asynchronous push-sum with growing intercommunication intervals, in 2018 Annual American Control Conference (ACC), IEEE, 2018, pp. 591–596.
  * [16] A. H. Sayed, Adaptation, learning, and optimization over networks, Foundations and Trends® in Machine Learning, 7 (2014), pp. 311–801.
  * [17] A. Tahbaz-Salehi and A. Jadbabaie, Consensus over ergodic stationary graph processes, IEEE Transactions on automatic Control, 55 (2009), pp. 225–230.
  * [18] J. N. Tsitsiklis, Problems in decentralized decision making and computation, PhD thesis, Massachusetts Institute of Technology, 1984.
  * [19] J. N. Tsitsiklis and V. Blondel, The Lyapunov exponent and joint spectral radius of pairs of matrices are hard—when not impossible—to compute and to approximate, Mathematics of Control, Signals and Systems, 10 (1997), pp. 31–40.


# Computational Discovery of Two-Dimensional Rare-Earth Iodides: Promising
Ferrovalley Materials for Valleytronics

Abhishek Sharan  Department of Physics, Khalifa University of Science and
Technology, Abu Dhabi-127788, United Arab Emirates (UAE)  Stephan Lany
National Renewable Energy Laboratory, Golden, CO 80401, USA  Nirpendra Singh
<EMAIL_ADDRESS>Department of Physics, and Center for Catalysis and
Separation (CeCaS), Khalifa University of Science and Technology, Abu
Dhabi-127788, United Arab Emirates (UAE)

###### Abstract

Two-dimensional Ferrovalley materials with intrinsic valley polarization are
rare but highly promising for valley-based nonvolatile random access memory
and valley filter. Using Kinetically Limited Minimization (KLM), an
unconstrained crystal structure prediction algorithm, and prototype sampling
based on first-principles calculations, we have discovered 17 new Ferrovalley
materials (rare-earth iodides RI2, where R is a rare-earth element belonging
to Sc, Y, or La-Lu, and I is Iodine). The rare-earth iodides are layered and
demonstrate 2H, 1T, or 1Td phase as the ground-state in bulk, analogous to
transition metal dichalcogenides (TMDCs). The calculated exfoliation energy of
monolayers is comparable to that of graphene and TMDCs, suggesting possible
experimental synthesis. The monolayers in the 2H phase exhibit two-dimensional
ferromagnetism due to unpaired electrons in $d$ and $f$ orbitals. Throughout
the rare-earth series, $d$ bands show valley polarization at $K$ and
$\overline{K}$ points in the Brillouin zone near the Fermi level. Due to
strong magnetic exchange interaction and spin-orbit coupling, large intrinsic
valley polarization in the range of 15-143 meV without external stimuli is
observed, which can be tuned and enhanced by applying a biaxial strain. These
valleys can selectively be probed and manipulated for information storage and
processing, potentially offering superior performance beyond conventional
electronics and spintronics. We further show that the 2H ferromagnetic phase
of RI2 monolayers possesses non-zero Berry curvature and exhibits the valley
Hall effect with considerable anomalous Hall conductivity. Our work will
incite exploratory synthesis of the predicted Ferrovalley materials and their
application in valleytronics and beyond.

##  I Introduction

Two-dimensional (2D) magnetism was a massive challenge until it was
successfully demonstrated in layered Cr2Ge2Te6 [1] and CrI3 [2]. Magnetism
combined with electronic and optical properties can lead to novel gateways for
exploring magneto-electrical, magneto-optical, and valley properties of
materials [3, 4, 5, 6, 7]. In the last few years, 2D magnets have been used in
novel magnetic devices to enhance the performance of existing devices, such as
development of novel magnetic tunnel junctions (Fe3GeTe2/hBN/Fe3GeTe2) [5],
gate tunable Curie temperature in Fe3GeTe2 over the extensive range 100-300 K
[6], and giant tunnel magnetoresistance in CrI3 based heterostructure [7] to
name a few. Certain 2D magnetic materials exhibit an extra valley degree of
freedom for electrons in addition to charge and spin, which can be leveraged
for memory devices leading to faster and more efficient logic systems and
next-generation storage devices beyond conventional electronics and
spintronics. Valley polarization is central to successfully controlling and
manipulating valley degree of freedom, which requires breaking time-reversal
and inversion symmetries. Since it was first observed and exploited in
graphene [8], several approaches, such as optical pumping [9], magnetic
proximity effect [10, 11], magnetic doping [12], and external magnetic [13]
and electric field [14] have been used to achieve valley polarization in
transition metal dichalcogenides (TMDCs). However, these methods involve
additional external complexity, such as strong interaction with magnetic
substrates, dynamical processes for optical pumping, and external magnetic
field, which is cumbersome for robust manipulation of valleys in device
architecture.

2D Ferrovalley materials possess exciting properties like 2D magnetism and
intrinsic valley polarization and can be used in valley-based devices. These
materials allow control and manipulation of valleys without any external
perturbation. Till date only a handful of Ferrovalley materials have been
reported in literature, LaI2, PrI2 [15], VSe2 [16], GeSe [17], VAgP2Se6 [18],
CeI2 [19], LaBr2 [20], LaBrI [21], CeI2 [19], NbX2 (X=S, Se) [22], Nb3I8 [23],
TiVI6 [24] and GdI2 [25]. Many more promising candidate materials are yet to
be explored, and the expansion of the family of Ferrovalley materials will
accelerate the research interest in fabricating valley-based electronic
devices. Furthermore, newly discovered materials can be combined with other 2D
materials to form novel heterostructures exhibiting long-range ferromagnetism,
complex spin textures, optoelectronic properties, topological bands, and
ferroelectricity. This work explores the chemically rich space of rare-earth
iodide compounds for their Ferrovalley properties and valleytronics
applications. Using unconstrained ground-state crystal structure prediction
algorithm Kinetically Limited Minimization (KLM) and prototype sampling
approach combined with first-principles calculations, we computationally
predict several new two-dimensional Ferrovalley materials. Rare-earth iodides
(in RI2 stoichiometry, where R ranges from Sc, Y, or La-Lu) are stable layered
materials with ground-state structure either in 2H, 1T, or 1Td phase,
analogous to TMDCs. These materials are magnetic in bulk and monolayer form,
except for YI2, which is non-magnetic. The 2H phase in monolayer form exhibits
valley polarization in the 15-143 meV range due to strong magnetic exchange
interaction and spin-orbit coupling. In the next few sections, we discuss the
crystal structure, relative phase stability, exfoliation energy, electronic,
magnetic, and anomalous Hall properties of predicted 2D RI2 monolayers.

##  II Methods

The Kinetically Limited Minimization (KLM), an unconstrained crystal structure
prediction algorithm, which combines unbiased random sampling and sequential
atomic perturbations based on ground-state energy derived from the first
principle calculations, is employed to predict the crystal structure of RI2
materials. The KLM method has successfully been employed for the discovery of
novel ternary nitrides [26, 27], ternary oxynitrides [28], and recently 2-D
materials [15, 29]. The crystal structure search is performed over 100-150
independent seed structures with varying formula units from 1 to 4. The
initial random seed structures and subsequent perturbations enforce the
minimum interatomic distance between different atomic species, significantly
reducing the configurational space and prohibiting the inclusion of unphysical
structures. The list of minimum inter-atomic distances used in the crystal
structure search is included in the Supplementary Information (SI). The KLM
approach demonstrates that the ground-state is either in 2H, 1T, and 1Td phase
for all the RI2 materials. Additionally, prototype structure sampling of RI2
materials in 2H, 1T, and 1Td phases are performed to determine the low energy
competing for meta-stable phases, which in some cases is not captured by the
KLM method. The possible magnetic configurations, ferromagnetic (FM), anti-
ferromagnetic (AFM), and non-magnetic (NM) are also taken into account for the
ground-state crystal structure search. To address the intrinsic magnetism in
these materials, a $2\times 2\times 1$ supercell for the 2H phase, $2\times
2\times 2$ supercell for the 1T phase, and $1\times 2\times 1$ supercell for
the 1Td phase is constructed such that each unit cell contains two layers. For
AFM configuration, the magnetic moment in two layers is initialized such that
the nearest rare-earth atoms in two layers possess opposite magnetic moments.
For FM and NM configurations, the respective unit cell is used to estimate the
total energy.

The first-principles calculations are performed based on density functional
theory with meta-GGA SCAN functional [30] as implemented in the Vienna ab
initio simulation package (VASP) [31, 32]. For faster convergence, generalized
gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) functional
[33] is used for crystal structure sampling. Then further analysis is
performed using SCAN functional on the determined structures. The interactions
between ionic cores and valence electrons are described with projector-
augmented-wave (PAW) potentials [34, 35]. The $3d$, $4s$ orbitals of Sc, $4s$,
$4p$, $4d$, $5s$ orbitals of Y, and $4f$, $5s$, $5p$, $5d$, $6s$ orbitals of
La-Lu rare earth elements are treated as valence electrons. The van der Waal
corrections are incorporated using the DFT-D3 method [36]. Energy cutoff of
520 eV for plane-wave basis set expansion and $\Gamma$ centered $k$-mesh of
$11\times 11\times 3$ for bulk and $11\times 11\times 1$ for monolayers for 2H
and 1T phases, and $7\times 12\times 2$ for bulk and $7\times 12\times 1$ for
monolayer for 1Td phase is used for integration over Brillouin zone. The
atomic positions and lattice constant are optimized until forces are less than
0.01 eV/ Å on each atom and total energy converged to less than 10-6 eV. A
vacuum of 20 Å is used along the out-of-plane direction for monolayers to
avoid interactions with its periodic images. The phonon band structure is
computed using the density functional perturbation theory [37] as implemented
in Phonopy code [38]. The Berry curvature and anomalous Hall conductivity are
calculated with WannierTools package [39] with tight binding Hamiltonian
extracted from WANNIER90 [40].

##  III Results

###  III.1 Crystal Structure

Unconstrained crystal structure search taking into account various possible
magnetic configurations, is executed for all possible rare-earth iodides with
the general formula RI2, where R is the rare-earth element from Sc, Y, or La-
Lu and I is Iodine. We find that rare-earth iodides in RI2 stoichiometry are
layered materials either in trigonal prismatic (2H) phase belonging to
$P6_{3}/mmc$ space group, octahedral (1T) phase belonging to $P\overline{3}m1$
space group, or distorted octahedral (1Td) phase belonging to $Pnm2_{1}$ space
group, as shown in Fig 1 (a), (b) and (c), respectively, analogous to TMDCs.
Additionally, prototype structure sampling of RI2 in 2H, 1T, and 1Td phase
shows the relative stability of the polytypes (see Fig 1(d)). The 2H and 1T
phases are composed of a hexagonal lattice with 2 and 1 formula unit per unit
cell, respectively, as shown by dashed lines in Fig 1(a) and (b). 1Td phase
can be generated from 1T phase by reconstruction in a 2$\times$1 orthorhombic
cell and is composed of 4 formula units per unit cell. The reconstruction
originated from a dimerized line of rare-earth atoms that leads to a lower
energy structure, driven by distortion due to the Peierls transition mechanism
[41, 42]. Similar distortion leading to the stability of the 1Td phase in
TMDCs has also been observed [43, 44]. All three phases are layered structures
with van der Waals interactions between two subsequent layers. The equilibrium
lattice parameters for each of the RI2 in the three phases are determined and
listed in Table 1. The lattice parameters lie in a very narrow range
throughout the series, (in-plane $a$ = 3.93 Å - 4.60 Å and out-of-plane per
layer $c$ = 6.86 Å - 8.18 Å ), owing to similar atomic radius of the rare-
earth elements. The distortion of the 1Td phase is not locally stable for all
materials. In the case of PrI2, NdI2, SmI2, EuI2, DyI2, ErI2 and TmI2, the
initial 1Td structure relaxes into the 1T structure, suggesting that the 1Td
phase is not even meta-stable for these materials.

###  III.2 Phase Stability

Relative phase stability of RI2 bulk compounds is shown in Fig 1(d), where the
phase stability is plotted in the space of $\mathrm{\Delta E_{2H-1T_{d}}}$ on
$x$-axis and $\mathrm{\Delta E_{1T-1T_{d}}}$ on $y$-axis. The possible
magnetic configurations, FM, AFM, and NM, are also taken into account for each
RI2 compound and the magnetic ground-state of each case is given in Table 1.
Out of the 17 RI2 compounds, ScI2, CeI2 and YbI2 have 1Td as the ground-state,
while 2H and 1T phases are meta-stable. The lowest energy structure of YI2,
LaI2, PmI2, GdI2, TbI2, ErI2 and LuI2 is 2H phase, while 1T and 1Td are meta-
stable phases except for ErI2 where 1Td phase is unstable. For the case of
PrI2, NdI2, SmI2, EuI2, DyI2, HoI2, and TmI2 the ground-state is 1T phase
while 2H phase makes meta-stable phase and 1Td is meta-stable only for HoI2.
For PrI2, NdI2, SmI2, EuI2, DyI2, ErI2, and TmI2 that lie on $x$-axis in phase
stability plot, Fig 1(d), the 1Td is unstable. The 1Td phase transforms into
lower energy 1T phase upon structural relaxation, hence $\mathrm{\Delta
E_{1T-1T_{d}}}=0$.

All the bulk RI2 compounds have FM as the ground-state for magnetic
configuration, except PmI2, DyI2, HoI2, and ErI2, which exhibit AFM as the
ground-state, and YbI2 is NM due to filled $4f$ orbitals in Yb+2. The $\Delta
E_{\mathrm{{2H-1T_{d}}}}$ and $\Delta E_{\mathrm{{1T-1T_{d}}}}$ in Table 1
show slight energy differences between different structural phases, which
suggests that synthesis of metastable phases may be possible at finite
temperatures under appropriate experimental conditions. The thermodynamic
stability of RI2 against their competing phases, including the RI3 phase, is
investigated using convex hull analysis and computing the decomposition energy
($\Delta E\mathrm{{}_{d}}$), shown in Fig 1(e) and (f). The $\Delta
E\mathrm{{}_{d}}$ is defined as the energy difference of the lowest energy RI2
phase from the convex hull formed by all the known phases (including elemental
phase) between respective R and I atoms, shown in Fig 1(e) for Nd-I system. A
negative $\Delta E\mathrm{{}_{d}}$ suggests that the energy of the lowest
energy RI2 phase falls below the convex hull, as for NdI2 in Fig 1(e), and
that it is stable against decomposition into its competing phases and can be
experimentally synthesized, which is the case for all the RI2 compounds
predicted in this work (see Fig 1(f)). It is to be noted that a three-
dimensional phase of EuI2 (space group 62 $Pnma$), which is already known
[45], is lower in energy than the 1T-FM phase by 32 meV per formula unit.
Towards I rich conditions, the RI3 phase is an important composition for
competing phases, as illustrated in Fig 1(e). The ground-state energy of all
the known phases in different stoichiometries is computed for the convex hull
analysis. RI3 compounds are known for all the rare-earth elements except for
Eu, Gd, and Yb. Therefore, we perform crystal structure sampling for EuI3,
GdI3, and YbI3, using the KLM approach. We find that all three materials are
stable with layered crystal structures belonging to space groups 2
($P\overline{1}$), 162 ($P\overline{3}1m$), and 12 ($C2/m$) for EuI3, GdI3 and
YbI3 respectively. The crystallographic information file for EuI3, GdI3, and
YbI3 structures are available in SI.

Rare-earth iodides in RI3 stoichiometry are also stable as layered materials,
where rare-earth elements are in a +3 oxidation state. However, However, this
work shows that rare-earth elements in a low-valent +2 oxidation state are
also stable. R+2 has the electron configuration [Ar]$3d^{1}4s^{0}$ for Sc+2,
[Kr]$4d^{1}5s^{0}$ for Y+2, and [Xe]$4f^{0-14}5d^{0-1}6s^{0}$ for R+2 where R
is La-Lu. From the electron configuration, one can say that the rare-earth
elements are also early transition elements with zero or one $d$ electron in
the valence shell. As with other transition elements which exhibit a range of
oxidation states, rare-earth elements in R0, R+1, R+2 and R+3 oxidation states
are also expected to exist [46, 47]. Considering other transition metal ions
with similar electron configurations, for instance, Hf+3 and La+2 have the
same electron configuration of [Xe]$5d^{1}6s^{0}$, Hf+3 is known in several
compounds [48]. Therefore, the stability of rare-earth elements in a less
common +2 oxidation state is also expected, as shown in this work for RI2
compounds.

###  III.3 Exfoliation Energy

The exfoliation energy ($E_{\mathrm{exf}}$) of monolayers in different phases
is calculated and listed in Table 1 and also shown in Fig 2 (a). The
exfoliation energy is computed by the relation [49],

$E_{\mathrm{exf}}=(E_{\mathrm{mono}}-E_{\mathrm{bulk}})/A$

where $E_{\mathrm{mono}}$ and $E_{\mathrm{bulk}}$ are the energy of monolayer
and bulk per formula unit, respectively, and $A$ is the in-plane area of bulk.
The calculated $E_{\mathrm{exf}}$ of monolayers are comparable to the
exfoliation energy of graphene (12 meV/$\AA^{2}$) [50] and MoS2 (26
meV/$\AA^{2}$) [51], advocating their mechanical exfoliation. The
$E_{\mathrm{exf}}$ of 2H and 1T phases ranges between 16-23 meV/$\AA^{2}$,
which is close to the exfoliation energies of TMDCs. The $E_{\mathrm{exf}}$ of
distorted octahedral (1Td) phase is significantly lower than that in 2H or 1T
phase and ranges between 7-10 meV/$\AA^{2}$, which suggests that the
monolayers in 1Td phase are weakly bound to each other and are easily
exfoliable. The 1Td phase is derived from the 1T phase, where two neighbor
metal atoms distort the octahedra around them, thereby increasing the out-of-
plane lattice parameter due to Peierls distortion [41, 42], which effectively
reduces the binding of adjacent layers, thereby decreasing the exfoliation
energy.

Figure 1: Lateral and top view of (a) trigonal prismatic (2H), (b) octahedral (1T), and (c) distorted octahedral (1Td) phases of RI2 compound. The unit cell is shown with a dashed line. (d) shows the relative stability of three phases for different RI2 bulk compounds. $\Delta E_{\mathrm{{2H-1T_{d}}}}$ and $\Delta E_{\mathrm{{1T-1T_{d}}}}$ are plotted on the $x$ and $y$ axis respectively. The stability region for each phase is highlighted. (e) convex hull diagram of Nd-I system schematically showing the definition of decomposition energy, $\Delta E\mathrm{{}_{d}}$. (f) shows $\Delta E\mathrm{{}_{d}}$ in meV per atom, for different RI2 bulk compounds. Table 1: Calculated lattice constants, ($a_{0}$, $b_{0}$ and $c_{0}$) in Å and exfoliation energy ($E_{\mathrm{exf}}$) in $\mathrm{meV/\AA^{2}}$ for 2H, 1T and 1Td phases of RI2 in bulk and monolayer form (for 2H and 1T phase $a_{0}$ = $b_{0}$). The energy difference between 2H and 1Td ($\Delta E_{\mathrm{2H-1T_{d}}}$), and between 1T and 1Td ($\Delta E_{\mathrm{1T-1T_{d}}}$) phase of bulk RI2 is given in meV per formula unit (fu). The ground-state (GS) and corresponding magnetic configuration for bulk RI2 (FM for ferromagnetic, AFM for anti-ferromagnetic, and NM for non-magnetic) is included. The decomposition energy, $\Delta E_{\mathrm{d}}$ in meV per atom, for each compound is also listed. |  | 2H Phase | 1T Phase | 1Td Phase | $\Delta E_{\mathrm{{2H-1T_{d}}}}$ | $\Delta E_{\mathrm{{1T-1T_{d}}}}$ | GS | $\Delta E_{\mathrm{d}}$
---|---|---|---|---|---|---|---|---
|  | $a_{0}$ | $c_{0}$ | Eexf | $a_{0}$ | $c_{0}$ | Eexf | $a_{0}$ | $b_{0}$ | $c_{0}$ | Eexf | (meV/fu) | (meV/fu) |  | (meV/at)
ScI2 | bulk | 3.93 | 14.60 |  | 4.01 | 7.05 |  | 4.02 | 6.90 | 14.35 |  | 64 | 93 | 1Td \- FM | -185
ML | 3.93 |  | 21 | 4.01 |  | 22 | 4.02 | 6.90 |  | 10 |  |  |
YI2 | bulk | 4.06 | 15.10 |  | 4.12 | 7.31 |  | 4.09 | 7.18 | 15.03 |  | -79 | 123 | 2H - FM | -56
ML | 4.06 |  | 20 | 4.12 |  | 21 | 4.09 | 7.18 |  | 9 |  |  |
LaI2 | bulk | 4.26 | 15.35 |  | 4.26 | 7.50 |  | 4.25 | 7.55 | 15.27 |  | -96 | 126 | 2H - FM | -24
ML | 4.26 |  | 19 | 4.26 |  | 19 | 4.25 | 7.55 |  | 9 |  |  |
CeI2 | bulk | 4.19 | 15.30 |  | 4.12 | 7.64 |  | 4.25 | 7.42 | 15.11 |  | 8 | 143 | 1Td \- FM | -29
ML | 4.19 |  | 16 | 4.12 |  | 10 | 4.25 | 7.42 |  | 9 |  |  |
PrI2 | bulk | 4.28 | 15.14 |  | 4.45 | 7.24 |  | unstable | 24 | 0 | 1T - FM | -42
ML | 4.28 |  | 18 | 4.45 |  | 18 |  |  |  |  |  |  |
NdI2 | bulk | 4.35 | 14.93 |  | 4.56 | 7.00 |  | unstable | 119 | 0 | 1T - FM | -345
ML | 4.35 |  | 16 | 4.56 |  | 17 |  |  |  |  |  |  |
PmI2 | bulk | 4.22 | 14.77 |  | 4.36 | 7.22 |  | 4.28 | 7.40 | 14.91 |  | -102 | -53 | 2H - AFM | -224
ML | 4.22 |  | 19 | 4.36 |  | 21 | 4.28 | 7.40 |  | 8 |  |  |
SmI2 | bulk | 4.35 | 14.81 |  | 4.57 | 6.94 |  | unstable | 177 | 0 | 1T - FM | -342
ML | 4.35 |  | 17 | 4.57 |  | 17 |  |  |  |  |  |  |
EuI2 | bulk | 4.43 | 14.78 |  | 4.60 | 6.89 |  | unstable | 275 | 0 | 1T - FM | -629
ML | 4.43 |  | 16 | 4.60 |  | 16 |  |  |  |  |  |  |
GdI2 | bulk | 4.09 | 14.99 |  | 4.06 | 7.45 |  | 4.11 | 7.18 | 14.99 |  | -60 | 204 | 2H - FM | -116
ML | 4.09 |  | 21 | 4.06 |  | 22 | 4.11 | 7.18 |  | 10 |  |  |
TbI2 | bulk | 4.05 | 15.01 |  | 4.09 | 7.33 |  | 4.29 | 7.56 | 16.32 |  | -97 | 34 | 2H - FM | -198
ML | 4.05 |  | 22 | 4.09 |  | 22 | 4.29 | 7.56 |  | 8 |  |  |
DyI2 | bulk | 4.22 | 14.81 |  | 4.37 | 7.05 |  | unstable | 86 | 0 | 1T - AFM | -276
ML | 4.22 |  | 18 | 4.37 |  | 19 |  |  |  |  |  |  |
HoI2 | bulk | 4.25 | 14.73 |  | 4.32 | 6.86 |  | 4.38 | 7.47 | 14.46 |  | 142 | -77 | 1T - AFM | -348
ML | 4.25 |  | 17 | 4.32 |  | 18 | 4.38 | 7.47 |  | 7 |  |  |
ErI2 | bulk | 4.11 | 14.56 |  | 4.22 | 6.99 |  | unstable | -42 | 0 | 2H - AFM | -91
ML | 4.11 |  | 23 | 4.22 |  | 21 |  |  |  |  |  |  |
TmI2 | bulk | 4.13 | 14.66 |  | 4.26 | 7.09 |  | unstable | 135 | 0 | 1T - FM | -290
ML | 4.13 |  | 20 | 4.26 |  | 19 |  |  |  |  |  |  |
YbI2 | bulk | 3.96 | 14.79 |  | 4.02 | 7.20 |  | 3.99 | 7.00 | 14.85 |  | 488 | 219 | 1Td \- NM | -769
ML | 3.96 |  | 21 | 4.02 |  | 21 | 3.99 | 7.00 |  | 8 |  |  |
LuI2 | bulk | 3.97 | 14.91 |  | 4.02 | 7.16 |  | 3.99 | 6.98 | 14.95 |  | -43 | 141 | 2H - FM | -54
ML | 3.97 |  | 22 | 4.02 |  | 23 | 3.99 | 6.98 |  | 10 |  |  |
Figure 2: (a) Exfoliation energy ($E_{\mathrm{exf}}$) of monolayer RI2 in 2H,
1T and 1Td phase, (b) splitting of $d$ orbitals in trigonal prismatic crystal
field of 2H phase, and (c) Valley polarization $\Delta E_{\mathrm{valley}}$
for different RI2 monolayers in 2H-FM phase.

###  III.4 Electronic Properties of 2H-monolayer

The electronic band structure of RI2 monolayer (ML) in the 2H phase shows
valley polarization and is interesting for valleytronics. Therefore, the
discussion will focus mainly on the 2H-ML phase of RI2 compounds unless
otherwise stated. The 2H phase in ML form is the ground-state for YI2, LaI2,
CeI2, PmI2, GdI2, TbI2, ErI2, TmI2 and LuI2 MLs, while 1T phase of PrI2, NdI2,
SmI2, EuI2, DyI2, HoI2 and YbI2 MLs is the ground-state and 1Td phase of ScI2
ML has the lowest energy. All the RI2 MLs favor FM configuration, except for
DyI2 (AFM), ErI2 (AFM), and YbI2 (NM), as mentioned in Table 2. For a few
cases, the magnetic configuration differs from the respective bulk ground-
state of RI2 compounds, as discussed in section III.2 and listed in Table 1.
The energy difference between the lowest energy phase and the 2H$-$FM
configuration for MLs, $\Delta E_{\mathrm{2H,FM}}$, is shown for comparison in
Table 2. The energy difference between stable and metastable phases
considering different magnetic configurations is comparable to room
temperature for several cases, suggesting that, under suitable experimental
conditions, the 2H phase in monolayer form can potentially be synthesized for
all the RI2 compounds. The dynamical stability of RI2 MLs in the 2H phase is
demonstrated by the absence of negative frequencies in phonon dispersion
plots, shown in Fig S1.

In the trigonal prismatic crystal field of 2H phase, each R atom is
coordinated with 6 I atoms. The $d$ orbitals of R element split into three
groups $a$($d_{z2}$), $e_{1}$($d_{x2-y2,xy}$) and $e_{2}$($d_{xz,yz}$) with
increasing energy as shown in Fig 2(b). The $3d_{z2}$ of Sc+2 and $4d_{z2}$ of
Y+2 are singly occupied in spin-up channel giving rise to a magnetic moment of
1 $\mu_{B}$ per formula unit with FM ordering. For other RI2 monolayers, the
partial occupancy of $4f$ and $5d$ orbitals give rise to magnetism in these
monolayers with magnetic moments given in Table 2. Depending on the number of
electrons in $4f$ orbitals, the $5d_{z2}$ orbitals is either occupied or
unoccupied but lies near the Fermi energy ($E_{F}$). Irrespective of the
occupation of $d_{z2}$ orbitals, the dispersion of $d_{z2}$ state near the
$E_{F}$ gives rise to the valleys at $K$ and $\overline{K}$ points in the
Brillouin zone. The electronic band structure of ScI2, YI2, NdI2, and LuI2
monolayers in 2H$-$FM configuration, with and without spin-orbit coupling
(SOC) is shown in Fig 3. With SOC, ScI2, YI2, NdI2, and LuI2 monolayers in 2H
phase are FM semiconductors with band gap of 0.23 eV, 0.32 eV, 0.10 eV, and
0.19 eV, and magnetic moments of 0.7 $\mu_{B}$, 0.5 $\mu_{B}$, 3.5 $\mu_{B}$
and 0.5 $\mu_{B}$ per formula unit, respectively. The band structure of other
RI2 MLs in 2H phase are shown in Fig S3-S6. LaI2, CeI2, PmI2, SmI2, EuI2, GdI2
and LuI2 monolayers are FM semiconductors with small band gap varying from
0.08 eV (LaI2) to 0.77 eV (EuI2). PrI2, TbI2, DyI2, HoI2, ErI2 and TmI2 are FM
metallic due to partial occupation of $4f$ and $5d$ orbitals. YbI2 is a NM
semiconductor due to fully filled $4f$ orbitals of Yb+2. The indirect band gap
$E_{g}$ of all the RI2 MLs in 2H phase is given in Table 2. The direct band
gap at $K$, $E_{g,K}$ is also given in Table 2 and schematically shown in Fig
3(e). Without SOC, the valleys at $K$ and $\overline{K}$ are degenerate, as
highlighted by dashed lines in Fig 3 (a)-(d). The SOC lifts the degeneracy,
leading to a large spontaneous valley polarization ($\Delta
E_{\mathrm{valley}}$) of 89 meV, 103 meV, 81 meV and 118 meV in ScI2, YI2,
NdI2, and LuI2 MLs respectively, also shown in Fig 3(e)-(f). A large valley
polarization is also seen in other RI2 MLs varying in the range of 15-143 meV,
as shown in their electronic band structures Fig S3-S6, listed in Table 2 and
shown Fig 2(c). The largest value of 143 meV for $\Delta E_{\mathrm{valley}}$
is for GdI2 monolayer which also possess largest magnetic moment of 7.4
$\mu_{B}$ in FM configuration. Although the indirect band gap in these MLs is
small, there is a large separation between the valence band maximum at $K$ and
conduction band minimum at $M$, ensuring the robustness of valleys against
scattering by long-wavelength phonons. Therefore, these materials exhibit
intrinsic magnetism in two dimensions and show large spontaneous valley
polarization without any external stimuli. Additionally, these valleys are
spin-polarized and coupled to the magnetic ordering, i.e., the polarity of
valleys can be reversed by changing the direction of intrinsic magnetism.

In certain cases, due to the relative ordering of $5d_{z2}$ and $4f$ orbitals,
the $5d_{z2}$ orbitals lie slightly away from $E_{F}$. Such as, in case of
PrI2 monolayer the valley band is $\sim$ 1 eV below the $E_{F}$ (Fig S3(c) and
(g)), while for NdI2 monolayer the valley band is $\sim$ 0.5 eV above $E_{F}$
(Fig 3(c) and (g)). The $E_{F}$ needs to be tuned to access the valley states
in such cases. Such $E_{F}$ tuning has been successfully achieved in Bi2Se3,
and Ca / Sn doped Sb2Te2Se [52], without modifying the states near the
$E_{F}$. The valleys in SmI2, TbI2, DyI2 and TmI2 MLs are difficult to access
due to the presence of $4f$ states near the $E_{F}$ (Fig S3-S6). Nevertheless,
we have included the valley polarization in Table 2 and Fig 2(c) for all the
RI2 monolayers in 2H-FM configuration. The valley polarization can further be
enhanced and tuned by the application of biaxial strain [15]. The spin-
polarized valleys give an extra degree of freedom for the electrons to
possess, which can potentially be utilized for enhanced functionalities of
next-generation devices. The traditional valley materials require an external
source such as magnetic proximity effect [53], optical pumping [9], or
magnetic field [13] to achieve valley polarization. In contrast, RI2
Ferrovalley materials eliminate the need for any external effect and are
easier to implement in device architecture, making them emerging materials for
next-generation nonvolatile memory devices [54, 55].

It is worth mentioning that the value of valley polarization and position of
$4f$ and $5d$ states of a rare-earth element is dependent on the exchange-
correlation functional. For instance, when using Generalized Gradient
Approximation (GGA) functional of Perdew-Burke-Ernzerhof (PBE) [33], valley
polarization of 10 meV and 55 meV is computed for LaI2 and PrI2 MLs
respectively, compared to 28 meV and 35 meV using meta-GGA SCAN functional
[30]. The Hubbard potential parameter ($U$) correctly describes crystal
lattice, and electronic properties of rare-earth or transition metal [56].
However, the $U$ value is often used as a fitting parameter, and its value is
ambiguous in describing the properties of a novel system. To avoid such
ambiguity, we have reported the results based on SCAN functional, which better
describes the properties of correlated systems than the PBE functional [57,
58]. However, it is critical to mention here that although the quantitative
description of valley polarization reported in this work may lie within a
certain energy window, due to the limitation of exchange-correlation
functionals in the accurate description of $4f$ and $5d$ state of rare-earth
elements, a qualitative description of the discovery of magnetism in two
dimensions and Ferrovalley property in these materials and trends in valley
polarization across the rare-earth series certainly holds true.

Figure 3: Calculated band structure of (a) ScI2, (b) YI2, (c) NdI2, and (d) LuI2 monolayers without SOC in 2H-FM configuration. Spin up and spin down bands are plotted in red and blue colors, respectively. Valley degeneracy at $K$ and $\overline{K}$ points are indicated by a dashed line in each case. The inset in (a) shows the high symmetry lines in the Brillouin zone along which the bandstructures are plotted. The $3d_{z2}$, $4d_{z2}$, $5d_{z2}$ and $4f$ states for ScI2, YI2, NdI2 and Lu2 is also shown in (a)-(d). (e), (f), (g) and (h) shows the SOC band structure of ScI2, YI2, NdI2, and LuI2 in 2H-FM configuration. The broken valley degeneracy is highlighted in the zoomed in inset in each case. Valley polarization values are also mentioned. The $E_{\mathrm{g}}$ and $E_{\mathrm{g,K}}$ are schematically shown in (e). Table 2: Valley polarization ($\Delta E_{\mathrm{valley}}$), global band gap ($E_{\mathrm{g}}$), direct band gap at $K$ point ($E_{\mathrm{g,K}}$), and magnetic moment ($\mu_{B}$) per formula unit for RI2 monolayers in 2H-FM configuration. The ground-state ($GS_{\mathrm{ML}}$) of each monolayer is given and ground-state energy with respect to 2H-FM configuration ($\Delta E_{\mathrm{GS-2H,FM}}$) in meV per formula unit is also shown. |  $\Delta E_{\mathrm{valley(2H,FM)}}$ (meV) |  $E_{\mathrm{g(2H,FM)}}$ (eV) |  $E_{\mathrm{g,K(2H,FM)}}$ (eV) |  $(\mu_{B}$ per fu)2H,FM | $GS_{\mathrm{ML}}$ |  $\Delta E_{\mathrm{GS-2H,FM}}$ (meV/fu)
---|---|---|---|---|---|---
ScI2 | 89 | 0.23 | 0.86 | 0.7 | 1Td \- FM | -76
YI2 | 103 | 0.32 | 1.24 | 0.5 | 2H - FM | 0
LaI2 | 28 | 0.08 | 0.74 | 0.4 | 2H - FM | 0
CeI2 | 15 | 0.39 | 0.50 | 1.4 | 2H - FM | 0
PrI2 | 35 | metallic | 0.76 | 2.5 | 1T - FM | -8
NdI2 | 81 | 0.10 | 0.79 | 3.5 | 1T - FM | -86
PmI2 | 68 | 0.22 | 1.00 | 4.5 | 2H - FM | 0
SmI2 | 40 | 0.18 | 0.14 | 5.5 | 1T - FM | -160
EuI2 | 115 | 0.77 | 0.60 | 6.8 | 1T - FM | -258
GdI2 | 143 | 0.60 | 0.58 | 7.4 | 2H - FM | 0
TbI2 | 114 | metallic | 0.05 | 6.1 | 2H - FM | 0
DyI2 | 132 | metallic | 0.15 | 4.5 | 1T - AFM | -70
HoI2 | 116 | metallic | 0.69 | 3.4 | 1T - FM | -184
ErI2 | 109 | metallic | 0.86 | 2.7 | 2H - AFM | -668
TmI2 | 51 | metallic | 0.39 | 1.2 | 2H - FM | 0
YbI2 | 0 | 0.31 | 0.00 | 0.0 | 1T - NM | -265
LuI2 | 118 | 0.19 | 1.25 | 0.5 | 2H - FM | 0

###  III.5 Valley Hall Effect

In 2D hexagonal magnetic systems with broken inversion and time-reversal
symmetries, the charge carriers possess non-zero Berry curvature
$\Omega_{z}(k)$ at $K$ and $\overline{K}$ points in the Brillouin zone,
leading to valley contrasting features in electron transport, known as the
valley Hall effect [59, 8, 60, 61]. The impact of an external electric field
on electron transport is investigated by calculating the out of plane ($z$)
component of the Berry curvature using the Kubo formula [62]:

$\Omega_{z}(k)=-\sum_{m}\sum_{m\neq
n}f(m)\frac{\mathrm{2Im}\bra{\phi_{mk}}\nu_{x}\ket{\phi_{nk}}\bra{\phi_{nk}}\nu_{y}\ket{\phi_{mk}}}{(E_{m}(k)-E_{n}(k))^{2}},$
(1)

where $m$ and $n$ are the band index for the occupied orbitals, $f(m)$ is the
Fermi-Dirac distribution function for the $m^{th}$ band, $\nu_{x}(y)$ is the
velocity operator of electrons along $x(y)$ directions, and $\ket{\phi_{mk}}$
and $\ket{\phi_{nk}}$ are the periodic part of the Bloch wavefunctions with
eigenvalues $E_{m}(k)$ and $E_{n}(k)$ respectively. The valley contrasting
features in $\Omega_{z}(k)$ for ScI2, YI2, NdI2, and LuI2 2H MLs is shown in
Fig 4. The peaks of $\Omega_{z}(k)$ that are opposite in sign and unequal in
magnitude can be seen at $K$ and $\overline{K}$ points. Under the effect of an
external in-plane electric field, $E$, the electrons in these monolayers
acquire velocity along transverse direction, $v_{\bot}\sim E_{\parallel}$
$\times$ $\Omega_{z}(k)$, giving rise to anomalous Hall conductivity
($\sigma_{xy}$). The $\sigma_{xy}$ is a function of the $E_{F}$ and is
computed as follows:

$\sigma_{xy}=\frac{e^{2}}{\hbar}\int_{BZ}\frac{d^{2}k}{(2\pi)^{2}}f(k)\Omega_{z}(k)$
(2)

where $f(k)$ is the Fermi-Dirac function as a function of reciprocal lattice
wave vector $k$.

The $\sigma_{xy}$ of ScI2, YI2, NdI2, and LuI2 monolayers as a function of the
$E_{F}$ is shown in Fig 4(e). The zero level means that the $d^{\downarrow}$
bands up to the valley (also highlighted in Fig 3) are occupied, and
$d^{\uparrow}$ bands are unoccupied. $\sigma_{xy}$ attains a large value when
the $E_{F}$ is negative and crosses the valleys. The similarity in the
variation of $\sigma_{xy}$ in ScI2, YI2, and LuI2 monolayers is a direct
consequence of the similar nature of $d_{z2}$ orbitals near the $E_{F}$.
$\sigma_{xy}$ attains large value near the near the $E_{F}$ in NdI2 than ScI2,
YI2 or LuI2 due to the contribution of $4f$ states near the $E_{F}$ to
$\sigma_{xy}$. The variation of $\sigma_{xy}$ of all other RI2 monolayers is
shown in Fig. S7 of SI. It is to be noted that valleys at the $\Gamma$ point
do not contribute to $\Omega_{z}(k)$ or $\sigma_{xy}$, and therefore, these
charge carriers will be undeflected in the transverse direction in a Hall
setup. Due to significant anomalous Hall conductivity in these monolayers, the
charge carriers along valleys can selectively be controlled and manipulated
for valleytronic applications.

Figure 4: Distribution of out of plane ($z$) component of Berry curvature
($\Omega_{z}$) for 2H-FM phase of (a) ScI2, (b) YI2, (c) NdI2, and (d) LuI2
monolayers within the reciprocal space. First Brillouin zone and high symmetry
points $K$ and $\overline{K}$ are highlighted. (e) shows the variation of
Anomalous Hall conductivity, $\sigma_{xy}$ as a function of $E_{F}$. The
region with $E_{F}$ from 0 to -0.25 eV is highlighted.

##  IV Conclusions

The unconstrained ground-state crystal structure prediction algorithm and
prototype sampling have been employed to discover new 2D rare earth iodides
(RI2 materials where R is a rare-earth element from Sc, Y, and La-Lu, and I is
Iodine). All the bulk materials are layered and either have trigonal prismatic
(2H-$P6_{3}/mmc$ space group), octahedral (1T- $P\overline{3}m1$ space group),
or distorted octahedral (1Td-$Pnm2_{1}$ space group) phase as the lowest
energy state. The thermodynamic stability analysis shows that all the
compounds are thermodynamically stable. The exfoliation energy of monolayers
is comparable to graphene and TMDCs, suggesting possible mechanical
exfoliation in the experiment. ScI2, YI2, LaI2, CeI2, NdI2, PmI2, SmI2, EuI2,
GdI2 and LuI2 monolayers in 2H phase are FM semiconductors with small band gap
varying from 0.08 eV (LaI2) to 0.77 eV (EuI2). PrI2, TbI2, DyI2, HoI2, ErI2
and TmI2 monolayers in 2H phase are FM metals, while YbI2 is a NM
semiconductor. RI2 monolayers exhibit large intrinsic valley polarization in
the range of 15-143 meV in the 2H-FM phase, which can further be tuned and
enhanced on the application of biaxial strain. These materials are novel
Ferrovalley materials, where electrons possess the additional valley degree of
freedom along with charge and spin. Beyond conventional electronics and
spintronics, these materials are emerging for next-generation electronic
devices for information processing and data storage. We further show that
these monolayers exhibit the valley Hall effect with considerable anomalous
Hall conductivity up to 30 S/cm with proper tuning of Fermi level. Because of
the rapid development of 2D materials and their applications, discovering new
2D materials with magnetism and large intrinsic valley polarization without
external stimuli is significant for valleytronics application. We believe that
incite exploratory synthesis, and the related applications of 2D rare earth
iodides in valleytronics will be highly sought after.

## Supplementary Material

See supplementary material for the structure sampling parameters, phonon
dispersion, details on convex hull analysis, electronic bandstructres of RI2
materials and crystallographic information file for EuI3, GdI3, and YbI3
structures. Table S1: Minimum inter-atomic distances $d_{min}$ used in KLM
structure sampling. Figure S1: Phonon dispersion of RI2 monolayers in 2H
phase. Figure S2: Convex Hull diagrams. Figure S3-S6: Electronic
bandstructures of RI2 monolayers in 2H-FM phase. Figure S7: Anomalous Hall
conductivity as a function of the Fermi level for RI2 monolayers in 2H-FM
phase.

## Acknowledgements

A. S. and N. S. acknowledge the financial support from Khalifa University of
Science and Technology under the startup grant FSU-2020-11/2020.

## References

  * Gong  _et al._ [2017] C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia, and X. Zhang, Discovery of intrinsic ferromagnetism in two-dimensional van der waals crystals, Nature 546, 265 (2017).
  * 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, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Layer-dependent ferromagnetism in a van der waals crystal down to the monolayer limit, Nature 546, 270 (2017).
  * Lan  _et al._ [2020] T. Lan, B. Ding, and B. Liu, Magneto-optic effect of two-dimensional materials and related applications, Nano Select 1, 298 (2020), https://onlinelibrary.wiley.com/doi/pdf/10.1002/nano.202000032 .
  * Burch [2018] K. S. Burch, Electric switching of magnetism in 2d, Nature Nanotechnology 13, 532 (2018).
  * Wang  _et al._ [2018] Z. Wang, D. Sapkota, T. Taniguchi, K. Watanabe, D. Mandrus, and A. F. Morpurgo, Tunneling spin valves based on Fe3GeTe2/hBN/Fe3GeTe2 van der waals heterostructures, Nano Letters 18, 4303 (2018).
  * Deng  _et al._ [2018] Y. Deng, Y. Yu, Y. Song, J. Zhang, N. Z. Wang, Z. Sun, Y. Yi, Y. Z. Wu, S. Wu, J. Zhu, J. Wang, X. H. Chen, and Y. Zhang, Gate-tunable room-temperature ferromagnetism in two-dimensional Fe3GeTe2, Nature 563, 94 (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, and P. Jarillo-Herrero, Probing magnetism in 2d van der waals crystalline insulators via electron tunneling, Science 360, 1218 (2018).
  * Rycerz  _et al._ [2007] A. Rycerz, J. Tworzydło, and C. W. J. Beenakker, Valley filter and valley valve in graphene, Nature Physics 3, 172 (2007).
  * Zeng  _et al._ [2012] H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Valley polarization in MoS2 monolayers by optical pumping, Nature Nanotechnology 7, 490 (2012).
  * Subhan and Hong [2020] F. Subhan and J. Hong, Large valley splitting and enhancement of curie temperature in a two-dimensional VI3/CrI3 heterostructure, The Journal of Physical Chemistry C 124, 7156 (2020).
  * Zhao  _et al._ [2017] C. Zhao, T. Norden, P. Zhang, P. Zhao, Y. Cheng, F. Sun, J. P. Parry, P. Taheri, J. Wang, Y. Yang, T. Scrace, K. Kang, S. Yang, G.-x. Miao, R. Sabirianov, G. Kioseoglou, W. Huang, A. Petrou, and H. Zeng, Enhanced valley splitting in monolayer WSe2 due to magnetic exchange field, Nature Nanotechnology 12, 757 (2017).
  * Singh and Schwingenschlögl [2017] N. Singh and U. Schwingenschlögl, A route to permanent valley polarization in monolayer MoS2, Advanced Materials 29, 1600970 (2017), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201600970 .
  * MacNeill  _et al._ [2015] D. MacNeill, C. Heikes, K. F. Mak, Z. Anderson, A. Kormányos, V. Zólyomi, J. Park, and D. C. Ralph, Breaking of valley degeneracy by magnetic field in monolayer MoSe2, Phys. Rev. Lett. 114, 037401 (2015).
  * Khan  _et al._ [2021] I. Khan, B. Marfoua, and J. Hong, Electric field induced giant valley polarization in two dimensional ferromagnetic WSe2/CrSnSe3 heterostructure, npj 2D Materials and Applications 5, 10 (2021).
  * [15] A. Sharan and N. Singh, Intrinsic valley polarization in computationally discovered two-dimensional ferrovalley materials: LaI2 and PrI2 monolayers, Advanced Theory and Simulations n/a, 2100476, https://onlinelibrary.wiley.com/doi/pdf/10.1002/adts.202100476 .
  * Tong  _et al._ [2016] W.-Y. Tong, S.-J. Gong, X. Wan, and C.-G. Duan, Concepts of ferrovalley material and anomalous valley hall effect, Nature Communications 7, 10.1038/ncomms13612 (2016).
  * Shen  _et al._ [2017] X.-W. Shen, W.-Y. Tong, S.-J. Gong, and C.-G. Duan, Electrically tunable polarizer based on 2d orthorhombic ferrovalley materials, 2D Materials 5, 011001 (2017).
  * Song  _et al._ [2018] Z. Song, X. Sun, J. Zheng, F. Pan, Y. Hou, M.-H. Yung, J. Yang, and J. Lu, Spontaneous valley splitting and valley pseudospin field effect transistors of monolayer VAgP2Se6, Nanoscale 10, 13986 (2018).
  * Sheng  _et al._ [2022] K. Sheng, Q. Chen, H.-K. Yuan, and Z.-Y. Wang, Monolayer CeI2: An intrinsic room-temperature ferrovalley semiconductor, Phys. Rev. B 105, 075304 (2022).
  * Zhao  _et al._ [2019] P. Zhao, Y. Ma, C. Lei, H. Wang, B. Huang, and Y. Dai, Single-layer LaBr2: Two-dimensional valleytronic semiconductor with spontaneous spin and valley polarizations, Applied Physics Letters 115, 261605 (2019), https://doi.org/10.1063/1.5129311 .
  * Jiang  _et al._ [2021] P. Jiang, L. Kang, Y.-L. Li, X. Zheng, Z. Zeng, and S. Sanvito, Prediction of the two-dimensional janus ferrovalley material LaBrI, Phys. Rev. B 104, 035430 (2021).
  * Zang  _et al._ [2020] Y. Zang, Y. Ma, R. Peng, H. Wang, B. Huang, and Y. Dai, Large valley-polarized state in single-layer NbX2 (X = S, Se): Theoretical prediction, Nano Research 14, 834 (2020).
  * Peng  _et al._ [2020] R. Peng, Y. Ma, X. Xu, Z. He, B. Huang, and Y. Dai, Intrinsic anomalous valley hall effect in single-layer $\mathrm{N}{\mathrm{b}}_{3}\mathrm{{I}}_{8}$, Phys. Rev. B 102, 035412 (2020).
  * Du  _et al._ [2020] W. Du, Y. Ma, R. Peng, H. Wang, B. Huang, and Y. Dai, Prediction of single-layer TiVI6 as a promising two-dimensional valleytronic semiconductor with spontaneous valley polarization, J. Mater. Chem. C 8, 13220 (2020).
  * Cheng  _et al._ [2021] H.-X. Cheng, J. Zhou, W. Ji, Y.-N. Zhang, and Y.-P. Feng, Two-dimensional intrinsic ferrovalley GdI2 with large valley polarization, Phys. Rev. B 103, 125121 (2021).
  * Sun  _et al._ [2019] W. Sun, C. J. Bartel, E. Arca, S. R. Bauers, B. Matthews, B. Orvañanos, B.-R. Chen, M. F. Toney, L. T. Schelhas, W. Tumas, J. Tate, A. Zakutayev, S. Lany, A. M. Holder, and G. Ceder, A map of the inorganic ternary metal nitrides, Nature Materials 18, 732 (2019).
  * Arca  _et al._ [2018] E. Arca, S. Lany, J. D. Perkins, C. Bartel, J. Mangum, W. Sun, A. Holder, G. Ceder, B. Gorman, G. Teeter, W. Tumas, and A. Zakutayev, Redox-mediated stabilization in zinc molybdenum nitrides, Journal of the American Chemical Society 140, 4293 (2018).
  * Sharan and Lany [2021] A. Sharan and S. Lany, Computational discovery of stable and metastable ternary oxynitrides, The Journal of Chemical Physics 154, 234706 (2021), https://doi.org/10.1063/5.0050356 .
  * Lany and Sharan [2021] S. Lany and A. Sharan, First principles predictions of SnO2/CdTe and SnO2/CdCl2/CdTe interface structures, in  _2021 IEEE 48th Photovoltaic Specialists Conference (PVSC)_ (2021) pp. 1464–1465.
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115, 036402 (2015).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Grimme  _et al._ [2011] S. Grimme, S. Ehrlich, and L. Goerigk, Effect of the damping function in dispersion corrected density functional theory, Journal of Computational Chemistry 32, 1456 (2011), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.21759 .
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scripta Materialia 108, 1 (2015).
  * Wu  _et al._ [2018] Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, WannierTools: An open-source software package for novel topological materials, Computer Physics Communications 224, 405 (2018).
  * Mostofi  _et al._ [2008] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, wannier90: A tool for obtaining maximally-localised wannier functions, Computer Physics Communications 178, 685 (2008).
  * Guyot and Dorn [1967] P. Guyot and J. E. Dorn, A critical review of the peierls mechanism, Canadian Journal of Physics 45, 983 (1967), https://doi.org/10.1139/p67-073 .
  * Lee  _et al._ [1973] P. A. Lee, T. M. Rice, and P. W. Anderson, Fluctuation effects at a peierls transition, Phys. Rev. Lett. 31, 462 (1973).
  * Besse  _et al._ [2018] R. Besse, N. A. M. S. Caturello, C. M. O. Bastos, D. Guedes-Sobrinho, M. P. Lima, G. M. Sipahi, and J. L. F. D. Silva, Size-induced phase evolution of MoSe2 nanoflakes revealed by density functional theory, The Journal of Physical Chemistry C 122, 20483 (2018).
  * Keum  _et al._ [2015] D. H. Keum, S. Cho, J. H. Kim, D.-H. Choe, H.-J. Sung, M. Kan, H. Kang, J.-Y. Hwang, S. W. Kim, H. Yang, K. J. Chang, and Y. H. Lee, Bandgap opening in few-layered monoclinic MoTe2, Nature Physics 11, 482 (2015).
  * Krings  _et al._ [2009] M. Krings, M. Wessel, and R. Dronskowski, EuI2, a low-temperature europium(II) iodide phase, Acta Crystallographica Section C 65, i66 (2009).
  * Woen and Evans [2016] D. H. Woen and W. J. Evans, Chapter 293 - expanding the +2 oxidation state of the rare-earth metals, uranium, and thorium in molecular complexes, in  _Handbook on the Physics and Chemistry of Rare Earths_, Including Actinides, Vol. 50 (Elsevier, 2016) pp. 337–394.
  * Nief [2010] F. Nief, Chapter 246 molecular chemistry of the rare-earth elements in uncommon low-valent states (Elsevier, 2010) pp. 241–300.
  * Fryzuk  _et al._ [1996] M. D. Fryzuk, M. Mylvaganam, M. J. Zawarotko, and L. R. MacGillivray, Synthesis, structure and reactivity of mononuclear paramagnetic complexes of zirconium and hafnium, Polyhedron 15, 689 (1996).
  * Jung  _et al._ [2018] J. H. Jung, C.-H. Park, and J. Ihm, A rigorous method of calculating exfoliation energies from first principles, Nano Letters 18, 2759 (2018).
  * Liu  _et al._ [2012] Z. Liu, J. Z. Liu, Y. Cheng, Z. Li, L. Wang, and Q. Zheng, Interlayer binding energy of graphite: A mesoscopic determination from deformation, Phys. Rev. B 85, 205418 (2012).
  * Zhou  _et al._ [2014] J.-J. Zhou, W. Feng, C.-C. Liu, S. Guan, and Y. Yao, Large-gap quantum spin hall insulator in single layer bismuth monobromide Bi4Br4, Nano Letters 14, 4767 (2014).
  * yong Wang  _et al._ [2011] Z. yong Wang, P. Wei, and J. Shi, Tuning the fermi level in Bi2Se3 bulk materials and transport devices, Frontiers of Physics 7, 160 (2011).
  * Qi  _et al._ [2015] J. Qi, X. Li, Q. Niu, and J. Feng, Giant and tunable valley degeneracy splitting in MoTe2, Phys. Rev. B 92, 121403 (2015).
  * Sakamoto  _et al._ [2013] K. Sakamoto, T.-H. Kim, T. Kuzumaki, B. Müller, Y. Yamamoto, M. Ohtaka, J. R. Osiecki, K. Miyamoto, Y. Takeichi, A. Harasawa, S. D. Stolwijk, A. B. Schmidt, J. Fujii, R. I. G. Uhrberg, M. Donath, H. W. Yeom, and T. Oda, Valley spin polarization by using the extraordinary rashba effect on silicon, Nature Communications 4, 2073 (2013).
  * Ahammed and De Sarkar [2022] R. Ahammed and A. De Sarkar, Valley spin polarization in two-dimensional $h\text{$-$}m\mathrm{N}(m=\mathrm{Nb},{Ta})$ monolayers: Merger of valleytronics with spintronics, Phys. Rev. B 105, 045426 (2022).
  * Li  _et al._ [2020] S. Li, Y. Li, M. Bäumer, and L. V. Moskaleva, Assessment of PBE+U and HSE06 methods and determination of optimal parameter U for the structural and energetic properties of rare earth oxides, The Journal of Chemical Physics 153, 164710 (2020), https://doi.org/10.1063/5.0024499 .
  * Kitchaev  _et al._ [2016] D. A. Kitchaev, H. Peng, Y. Liu, J. Sun, J. P. Perdew, and G. Ceder, Energetics of ${\mathrm{mno}}_{2}$ polymorphs in density functional theory, Phys. Rev. B 93, 045132 (2016).
  * Chakraborty  _et al._ [2018] A. Chakraborty, M. Dixit, D. Aurbach, and D. T. Major, Predicting accurate cathode properties of layered oxide materials using the SCAN meta-GGA density functional, npj Computational Materials 4, 60 (2018).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of MoS2 and other group-vi dichalcogenides, Phys. Rev. Lett. 108, 196802 (2012).
  * Mak  _et al._ [2014] K. F. Mak, K. L. McGill, J. Park, and P. L. McEuen, The valley hall effect in MoS2 transistors, Science 344, 1489 (2014).
  * Pan  _et al._ [2014] H. Pan, Z. Li, C.-C. Liu, G. Zhu, Z. Qiao, and Y. Yao, Valley-Polarized Quantum Anomalous Hall Effect in Silicene, Phys. Rev. Lett. 112, 106802 (2014).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized hall conductance in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405 (1982).