# Magnetic Interactions in BiFeO3: a First-Principles Study

Changsong Xu  Physics Department and Institute for Nanoscience and
Engineering, University of Arkansas, Fayetteville, Arkansas 72701, USA  Bin Xu
<EMAIL_ADDRESS>Physics Department and Institute for Nanoscience and
Engineering, University of Arkansas, Fayetteville, Arkansas 72701, USA  School
of Physical Science and Technology, Soochow University, Suzhou, Jiangsu
215006, China  Bertrand Dupé  Institute of Physics, INSPIRE Group, Johannes
Gutenberg-University Mainz, 55128 Mainz, Germany  L. Bellaiche
<EMAIL_ADDRESS>Physics Department and Institute for Nanoscience and
Engineering, University of Arkansas, Fayetteville, Arkansas 72701, USA

###### Abstract

First-principles calculations, in combination with the four-state energy
mapping method, are performed to extract the magnetic interaction parameters
of multiferroic BiFeO3. Such parameters include the symmetric exchange (SE)
couplings and the Dzyaloshinskii-Moriya (DM) interactions up to second nearest
neighbors, as well as the single ion anisotropy (SIA). All magnetic parameters
are obtained not only for the $R3c$ structural ground state, but also for the
$R3m$ and $R\bar{3}c$ phases in order to determine the effects of
ferroelectricity and antiferrodistortion distortions, respectively, on these
magnetic parameters. In particular, two different second-nearest neighbor
couplings are identified and their origins are discussed in details. Moreover,
Monte-Carlo (MC) simulations using a magnetic Hamiltonian incorporating these
first-principles-derived interaction parameters are further performed. They
result (i) not only in the accurate prediction of the spin-canted G-type
antiferromagnetic structure and of the known magnetic cycloid propagating
along a $<$1$\bar{1}$0$>$ direction, as well as their unusual characteristics
(such as a weak magnetization and spin-density-waves, respectively); (ii) but
also in the finding of another cycloidal state of low-energy and that awaits
to be experimentally confirmed. Turning on and off the different magnetic
interaction parameters in the MC simulations also reveal the precise role of
each of them on magnetism.

##  I Introduction

Bismuth ferrite BiFeO3 (BFO) is one of the most robust room-temperature
multiferroic compounds. Besides its large electric polarization, BFO exhibits
different magnetic phases. For instance, it can possess a long period cycloid
or a canted configuration in which a predominant G-type antiferromagnetism
(AFM) coexists with a weak ferromagnetic vectorSando  _et al._ (2013);
Albrecht  _et al._ (2010). Upon external stimuli, such as temperature, fields,
strain and pressure, such two magnetic states can transform from one to
anotherSando  _et al._ (2013); Popov  _et al._ (1993, 1994); Tokunaga  _et
al._ (2010); Agbelele  _et al._ (2017); Rovillain  _et al._ (2010); Popkov
_et al._ (2015); Sosnowska  _et al._ (2002); Buhot  _et al._ (2015), which
reflects spin-lattice couplings in BFO. More precisely, spins have been
predicted to couple with both ferroelectric (FE) displacements and FeO6
octahedral tiltings (also known as antiferrodistortive (AFD) motions) in BFO,
see, e.g., Ref. Rahmedov  _et al._ (2012) and references therein.

Such spin-lattice couplings form a fundamental and important research
direction, as evidenced by the fact that different models have been proposed
to describe them and the resulting magnetism in BFO. Examples of such models
include the spin current model,Xu  _et al._ (2018a); Katsura  _et al._ (2005);
Raeliarijaona  _et al._ (2013); Rahmedov  _et al._ (2012) theory for
electrical-field control of magnetism from R. de Sousa and collaborators de
Sousa and Moore (2008); de Sousa  _et al._ (2013); de Sousa (2013) and various
models from R. S. Fishman et. al. Fishman  _et al._ (2012, 2013); Fishman
(2018). However, to the best of our knowledge, the magnetic coupling
coefficients, especially the anisotropic ones (that are important to generate
complex magnetic configurations), have never been systematically and
thoroughly studied, especially from direct first principles.

Here, we consider an ab-initio effective Hamiltonian with all its coupling
coefficients being determined from first-principle techniques and adopting the
most general matrix form. Such matrices enable us not only to have a general
idea of the magnetic anisotropy, but also to obtain the individual
isotropic/anisotropic symmetric exchange (SE) couplings, Dzyaloshinskii-Moriya
(DM) interactionsDzyaloshinsky (1958); Moriya (1960) and the single anion
anisotropy (SIA) by decompositions of such matrices. The effect of FE and AFD
distortions on such couplings are also determined and discussed. The paper is
organized as follows. Section II introduces the magnetic matrices and their
decomposition, as well as provides details about our density functional theory
(DFT) calculations and the Monte-Carlo (MC) simulations. Moreover, subsections
III.A, III.B and III.C of Section III focus on first, second nearest neighbor
couplings and SIA, respectively, while Subsection III.D provides results from
MC simulations using the aforementioned ab-initio-based effective Hamiltonian.
A brief conclusion is given in Section IV.

##  II Method

###  II.1 Magnetic effective Hamiltonian

Let us first define our convention for the coordinates as (i) the $x$-, $y$-,
and $z$-axes being along the pseudocubic [100], [010] and [001] directions,
respectively; and (ii) the FE displacements and the AFD axis about which the
FeO6 octahedra rotate being both along the pseudo-cubic [111] direction – as
consistent with the $R3c$ rhombohedral ground state of BiFeO3Wang  _et al._
(2003); Diéguez  _et al._ (2011).

The following magnetic effective Hamiltonian, $\mathcal{H}$, is adopted here:

$\displaystyle\mathcal{H}$
$\displaystyle=\mathcal{H}_{1}^{ex}+\mathcal{H}_{2}^{ex}+\mathcal{H}^{si}$ (1)

with

$\mathcal{H}_{1}^{ex}=\frac{1}{2}\sum_{<i,j>_{1}}\bm{{\rm
S}}_{i}{\cdot}\mathcal{J}_{1,ij}{\cdot}\bm{{\rm S}}_{j}~{}~{}~{},$ (2)

$\displaystyle\mathcal{H}_{2}^{ex}=\frac{1}{2}\sum_{<i,j>_{2}}\bm{{\rm
S}}_{i}{\cdot}\mathcal{J}_{2,ij}{\cdot}$ $\displaystyle\bm{{\rm S}}_{j}$ (3)
$\displaystyle=\frac{1}{2}\sum_{<i,j>_{2}^{1}}\bm{{\rm
S}}_{i}{\cdot}\mathcal{J}_{2,ij}^{1}{\cdot}\bm{{\rm
S}}_{j}+\frac{1}{2}\sum_{<i,j>_{2}^{2}}$ $\displaystyle\bm{{\rm
S}}_{i}{\cdot}\mathcal{J}_{2,ij}^{2}{\cdot}\bm{{\rm S}}_{j}~{}~{}~{},$

and

$\mathcal{H}^{si}=\sum_{i}\bm{{\rm
S}}_{i}{\cdot}\mathcal{A}_{ii}{\cdot}\bm{{\rm S}}_{i}$ (4)

where $\mathcal{H}_{1}^{ex}$ and $\mathcal{H}_{2}^{ex}$ denote the exchange
coupling between first and second nearest neighbors, respectively, and
$\mathcal{H}^{si}$ represents SIA. Note that the sum over first nearest
neighbors $<$$i,j$$>_{1}$ are 6-fold degenerate along $<$100$>$ directions. On
the other hand, the 12 second nearest neighbors $<$$i,j$$>_{2}$ can be
categorized into two types, $<$$i,j$$>_{2}^{1}$ being 6-fold degenerate along
the $<$1$\bar{1}$0$>$ directions that are perpendicular to the [111]
polarization direction versus $<$$i,j$$>_{2}^{2}$ that is also 6-fold
degenerate but along the $<$110$>$ directions that are not perpendicular to
the polarization direction. Moreover, S = 5/2 is used here to be consistent
with the valence state of Fe3+ ions in BFO.

The $\mathcal{J}$ matrices characterizing the magnetic exchange couplings are
calculated in the most general 3$\times$3 matrix form as

$\mathcal{J}=\left({\begin{array}[]{*{20}{c}}{{J_{xx}}}&{{J_{xy}}}&{{J_{xz}}}\\\
{{J_{yx}}}&{{J_{yy}}}&{{J_{yz}}}\\\ {{J_{zx}}}&{{J_{zy}}}&{{J_{zz}}}\\\
\end{array}}\right).$

They can always be decomposed into a symmetric part $\mathcal{J}_{SE}$ and an
antisymmetric part $\mathcal{J}_{DM}$, i.e., $\mathcal{J}$ =
$\mathcal{J}_{SE}$ \+ $\mathcal{J}_{DM}$.

The symmetric $\mathcal{J}_{SE}$ is given by

$\mathcal{J}_{SE}=\left({\begin{array}[]{*{20}{c}}{{J_{xx}}}&{\frac{1}{2}({J_{xy}}+{J_{yx}})}&{\frac{1}{2}({J_{xz}}+{J_{zx}})}\\\
{\frac{1}{2}({J_{xy}}+{J_{yx}})}&{{J_{yy}}}&{\frac{1}{2}({J_{yz}}+{J_{zy}})}\\\
{\frac{1}{2}({J_{xz}}+{J_{zx}})}&{\frac{1}{2}({J_{yz}}+{J_{zy}})}&{{J_{zz}}}\\\
\end{array}}\right).$

The $\mathcal{J}_{SE}$ matrices prefer spins being collinearly aligned. Unless
the fully isotropic case, it prefers an easy axis or an easy plane, whose
direction or normal, respectively, can be determined by the diagonalization of
the $\mathcal{J}_{SE}$ matrices. We numerically found that the off-diagonal
elements of $\mathcal{J}_{SE}$ are negligible and we will thus only focus on
$J_{\alpha\alpha}$ ($\alpha=x,y$ and $z$). Note that $J>0$ favors
antiferromagnetism.

The antisymmetric $\mathcal{J}_{DM}$ matrices (which is related to the DM
interaction) can be obtained as

$\mathcal{J}_{DM}=\left({\begin{array}[]{*{20}{c}}0&{\frac{1}{2}({J_{xy}}-{J_{yx}})}&{\frac{1}{2}({J_{xz}}-{J_{zx}})}\\\
{\frac{1}{2}({J_{yx}}-{J_{xy}})}&0&{\frac{1}{2}({J_{yz}}-{J_{zy}})}\\\
{\frac{1}{2}({J_{zx}}-{J_{xz}})}&{\frac{1}{2}({J_{zy}}-{J_{yz}})}&0\\\
\end{array}}\right).$

Note that, typically, $\mathcal{J}_{DM}$ is written using the vector
$\mathbf{D}$ via
$\mathcal{H}_{DM}=\mathbf{D}{\cdot}(\mathbf{S}_{i}\times\mathbf{S}_{j})$, with

$\mathbf{D}=(D_{x},D_{y},D_{z})$

where $D_{x}=\frac{1}{2}({J_{yz}}-{J_{zy}})$,
$D_{y}=\frac{1}{2}({J_{zx}}-{J_{xz}})$ and
$D_{z}=\frac{1}{2}({J_{xy}}-{J_{yx}})$. $\mathcal{J}_{DM}$, or equivalently
$\mathbf{D}$, favors the spins being perpendicular to each other within the
plane for which the normal vector is parallel to $\mathbf{D}$.

It is necessary to further clarify the term of “exchange coupling”. The
exchange coupling in common sense is of the form $J\bm{{\rm
S}}_{i}{\cdot}\bm{{\rm S}}_{j}$, which leads to isotropic collinear spin
configurations. It is usually considered as an alternative concept to DM
interaction, as in $\mathbf{D}{\cdot}(\mathbf{S}_{i}\times\mathbf{S}_{j})$.
However, in this manuscript, we use a stricter terminology that exchange
coupling refers to the form of $\bm{{\rm
S}}_{i}{\cdot}\mathcal{J}{\cdot}\bm{{\rm S}}_{j}$, with $\mathcal{J}$
including a symmetric part $\mathcal{J}_{SE}$ and an antisymmetric part
$\mathcal{J}_{DM}$ (equivalent to $\mathbf{D}$), both of which can lead to
magnetic anisotropy.

Moreover and according to point group symmetry (3m for $R3c$, $R3m$ and
$\bar{3}$m for $R\bar{3}c$), the $\mathcal{A}$ matrices associated with SIA
for $R3c$, $R3m$ and $R\bar{3}c$ phases all have the form of

$\mathcal{A}=\left({\begin{array}[]{*{20}{c}}0&\Delta&\Delta\\\
\Delta&0&\Delta\\\ \Delta&\Delta&0\\\ \end{array}}\right)$

in the ($x$, $y$, $z$) basis. This $\mathcal{A}$ matrix can be rewritten in
its diagonalizing basis as:

$\mathcal{A}=\left({\begin{array}[]{*{20}{c}}-\Delta&0&0\\\ 0&-\Delta&0\\\
0&0&2\Delta\\\ \end{array}}\right).$

where the third index corresponds to the pseudo-cubic [111] direction, while
indices 1 and 2 are associated with perpendicular directions, such as
[1$\bar{1}$0] and [11$\bar{2}$]. As a result, SIA favors [111] (or
[$\bar{1}$$\bar{1}$$\bar{1}$]) for the spin directions if $\Delta<0$, while it
prefers spins lying inside the (111) plane if $\Delta>0$.

###  II.2 DFT parameters and MC simulations

DFT calculations are performed using the Vienna ab-initio simulation package
(VASP) Kresse and Joubert (1999). The projector augmented wave (PAW) method
Blöchl (1994) is employed with the following electrons being treated as
valence states: Bi 6$s$ and 6$p$, Fe 3$d$ and 4$s$, and O 2$s$ and 2$p$. The
revised Perdew, Burke, and Ernzerhof functional for solids (PBE_sol) Perdew
_et al._ (2008) is used, with a typical effective Hubbard $U$ parameter of 4
eV for the localized 3$d$ electrons of Fe ions Xu  _et al._ (2014); Diéguez
_et al._ (2011). $k$-point meshes are chosen such as they are commensurate
with the choice of 6$\times$6$\times$6 for the 5-atom cubic $Pm\bar{3}m$
phase. For instance, (i) the 10-atom $R3c$ phase is optimized using
4$\times$4$\times$4 $k$-mesh, until the Hellmann-Feynman forces are converged
to be smaller than 0.001 eV/Å on each ion (the $R3m$ and $R\bar{3}c$ phases
are obtained from the decomposition of the optimized $R3c$ phase, that is the
AFD (respectively, FE) displacements of the $R3c$ ground state are left out
when constructing the $R3m$ (respectively, $R\bar{3}c$) state); (ii) the
exchange coupling coefficients are calculated using a 4$\times$4$\times$2
supercell with an 1$\times$1$\times$3 $k$-mesh; and (iii) the SIA parameters
are calculated using a 2$\times$2$\times$2 supercell with a
3$\times$3$\times$3 $k$-mesh. Note that the G-type antiferromagnetism with the
canted ferromagnetism is adopted when optimizing $R3c$ structures. Spin-
orbital coupling and noncolinear magnetic configurations are employed
throughout all calculations (except for the results in Table III, see details
there). The magnetic coefficients are extracted using the four-state energy
mapping method, as detailed in Refs. Xiang  _et al._ (2013, 2011). We
calculate all matrices for different Fe-Fe pairs or Fe sites, and the elements
are displayed to the digit of 0.001 meV through the manuscript.

Monte Carlo simulations are performed using the heat bath algorithmMiyatake
_et al._ (1986). A 12$\times$12$\times$12 supercell are adopted to predict the
Néel temperature ($T_{N}$). The 10-atom primitive cell and 2$\times$2$\times$2
supercells are used to determine the effects of each single magnetic
parameter, while supercells with the form of $\sqrt{2}n\times\sqrt{2}\times 2$
($n$ = 2, 3,…, 240), in which the first axis is along the [1$\bar{1}$0]
direction, and $\sqrt{2}\times\sqrt{2}\times 2n$ ($n$ = 2,3,…,240), in which
the last axis lies along [001], are adopted to determine properties of
cycloidal phases that propagate along $[1\bar{1}0]$ and [001] directions,
respectively (note that we decided to look at cycloids propagating along the
unusual [001] direction because recent effective Hamiltonian computations Xu
_et al._ (2018a) predicted that such cycloids can be very close in energy from
that of the well-known cycloid of BFO propagating along $[1\bar{1}0]$). In
each MC simulation, 2,000 exchange stepsMiyatake  _et al._ (1986) are
performed, with each exchange step containing 200 MC sweeps.

##  III Results

The application of the aforementioned DFT parameters results in the $R3c$
structure with lattice parameters of $a=b=c=$ 5.584 Å and
$\alpha=\beta=\gamma=$ 59.529∘, as well as the internal positions of atoms
being Bi 2a (0.276, 0.276, 0.276), Fe 2a (0, 0, 0) and O 6c (0.672, 0.813,
0.217). Such lattice parameters are within 0.8% difference as compared to
previous calculations and measurementsKubel and Schmid (1990); Diéguez  _et
al._ (2011), which testify the accuracy of our DFT calculations.

###  III.1 First nearest neighbor coupling $\mathcal{J}_{1}$

Table 1: Calculated symmetric exchange parameters and DM interactions for the nearest neighbor Fe-Fe pair along the [100] direction. The isotropic coupling coefficient $J_{1}$ is the average of the diagonal $xx$, $yy$ and $zz$ components. Note that $\mathbf{D_{1}^{a}}$ and $\mathbf{D_{1}^{b}}$ has the form of (0,$\alpha$,-$\alpha)$ and ($\beta$,$\beta$,$\beta)$, respectively. $D_{1}$ is the norm of $\mathbf{D}_{1}$ (unit: meV).  [100] |  $J_{1,xx}$ |  $J_{1,yy}$ |  $J_{1,zz}$ |  $J_{1}$
---|---|---|---|---
$R3c$ |  6.076 |  6.090 |  6.091 |  6.086
$R3m$ |  7.414 |  7.435 |  7.436 |  7.428
$R\bar{3}c$ |  5.847 |  5.858 |  5.860 |  5.855
[100] |  $D_{1,x}$ |  $D_{1,y}$ |  $D_{1,z}$ |  $D_{1}$
$R3c$ |  $\mathbf{D_{1}}$ |  -0.042 |  0.028 |  -0.116 |  0.126
$\mathbf{D_{1}^{a}}$ |  0.000 |  0.072 |  -0.072 |  0.102
$\mathbf{D_{1}^{b}}$ |  -0.043 |  -0.043 |  -0.043 |  0.074
$R3m$ |  0.003 |  0.135 |  -0.136 |  0.192
$R\bar{3}c$ |  -0.077 |  -0.027 |  -0.027 |  0.086

Let us first focus on the nearest neighbor exchange coupling and choose the
Fe-Fe pair along the [100] direction as an example. As shown in Table I, the
isotropic $J_{1}$ (which is the average of $J_{1,xx}$, $J_{1,yy}$ and
$J_{1,zz}$) yields 6.086 meV, whose positive sign indicates that the coupling
is of AFM nature. Such parameter is rather close to the values of 6.48Matsuda
_et al._ (2012), 4.38Jeong  _et al._ (2012) and 4.34Xu  _et al._ (2012) meV
that are estimated from inelastic neutron scattering, which further attests
the accuracy of our calculations. Values of $J_{1}$ are also calculated for
the $R3m$ phase, that only adopts FE displacements, and the $R\bar{3}c$ phase,
that only possesses AFD distortions. The $J_{1}$ value for $R3m$ phase yields
a larger 7.428 meV, while that of $R\bar{3}c$ phase gives a smaller 5.855 meV.
Such comparison indicates that the FE displacements contribute more to the AFM
than the oxygen octahedral tilting does. Taking advantage of the general
$\mathcal{J}$ matrix, SE coupling is found to yield an easy plane that is
perpendicular to the pair direction in the $R3c$ structure, as $J_{1,yy}$
$\approx$ $J_{1,zz}$ = 6.091 meV, while $J_{1,xx}$ = 6.076 meV. Such energy
differences result in an easy plane that is perpendicular to the [111]
direction, when all six nearest neighbors are considered, which is consistent
with proposed directions of the AFM vector in the spin-canted structureSando
_et al._ (2013). Note that such anisotropic SE coupling has been recently
reported to be significant in LaMn3Cr4O12 and is responsible for inducing its
multiferroicity Feng and Xiang (2016). Similar anisotropic SE coupling is also
found in the $R3m$ and $R\bar{3}c$ phases.

Moreover, the DM vector for first nearest neighbors and in the ($x$, $y$, $z$)
basis is calculated to be $\mathbf{D_{1}}$ = (-0.042, 0.028, -0.116) for the
$R3c$ state, resulting in a magnitude $D_{1}$ of 0.126 meV – that is about 50
times smaller than $J_{1}$ (note that Ref. Chen  _et al._ (2018) provided a
much larger magnitude of $D_{1}$ that is equal to 0.193, 0.327 and 0.321 meV
for the three different $<001>$ pairs, which is surprising since all these
first nearest-neighbor pairs should have the same magnitude of $D_{1}$ in the
$R3c$ state. The overestimation of the magnitude of $D_{1}$ in Ref. Chen  _et
al._ (2018) with respect to our present results likely lies in the choice of
too small supercells used within the four-state method in Ref. Chen  _et al._
(2018)). As commonly done for magnetic HamiltoniansFishman  _et al._ (2012,
2013); Fishman (2018), $\mathbf{D_{1}}$ can be decomposed into two parts,
$\mathbf{D_{1}^{a}}$ (0, $\alpha$, -$\alpha$) that determines the cycloidal
plane and period $\lambda$Fishman (2018) and $\mathbf{D_{1}^{b}}$ ($\beta$,
$\beta$, $\beta$) that can either create components of spins forming a spin-
density wave and being away from the cycloidal plane Rahmedov  _et al._
(2012); Ramazanoglu  _et al._ (2011a) for the cycloidal configuration or to
the creation of a weak magnetization in the spin-canted structure Albrecht
_et al._ (2010); Bellaiche  _et al._ (2012); Ederer and Spaldin (2005). Here,
we found that $\alpha$ = 0.072 meV and $\beta$ = -0.043 meV. As a result,
$D_{1}^{a}$ possesses a magnitude of 0.102 meV and $D_{1}^{b}$ has a strength
of 0.074 meV. Such parameters are well consistent with the values of 0.18 meV
and 0.06 meV, respectively, which are estimated from previous experiments and
modelsFishman (2018); Sosnowska and Zvezdin (1995); Tokunaga  _et al._ (2010);
Ramazanoglu  _et al._ (2011b); Fishman  _et al._ (2013); Ruette  _et al._
(2004); Jeong  _et al._ (2014). Moreover, the $\mathbf{D_{1}}$ vector of $R3m$
is numerically determined to be (0.003, 0.135, -0.136) meV, that is close to
adopt the form of (0, A, -A). It therefore has mostly a $\mathbf{D_{1}^{a}}$
component, and, consequently, its $\mathbf{D_{1}^{b}}$ component is nearly
vanishing. Such fact implies that the $\mathbf{D_{1}^{b}}$ component in the
$R3c$ phase mostly originates from AFD tiltings. Such finding is consistent
with the expression of the DM effect proposed in Refs.Albrecht  _et al._
(2010); Bellaiche  _et al._ (2012), which involves the tiltings of first-
nearest-neighbors oxygen octahedra and which was suggested to be responsible
for the weak ferromagnetism in the spin canted structure of BFO. Such fact is
further confirmed by the fact that the $\mathbf{D_{1}}$ vector of $R\bar{3}c$
is found to be equal to (-0.077, -0.027, -0.027) meV and has therefore a (B,
C, C) form, which results in a $\mathbf{D_{1}^{b}}$ component that can be be
estimated to be (-0.043, -0.043, -0.043) meV when taking an average $\beta$ to
be equal to (B+2C)/3. Interestingly, this resulting $\mathbf{D_{1}^{b}}$
vector of $R\bar{3}c$ is precisely the one of the $R3c$ structure, which
further confirms that this latter originates from oxygen octahedral tilting
rather than polarization. On the other hand, polarization does contribute to
the $\mathbf{D_{1}^{a}}$ of the $R3c$ phase since the $\mathbf{D_{1}^{a}}$ of
the $R3m$ phase is significant. Such feature is in-line with spin-current
models involving the polarization, ${\bf P}$, and first-nearest neighbors for
the DM effect that has an energy of the form $C_{1}({\bf P}\times{\bf
e_{ij}})\cdot({\bf m_{i}}\times{\bf m_{j}})$, where $C_{1}$ is a material-
dependent coefficient, ${\bf e_{ij}}$ is the unit vector joining site $i$ to
site $j$ and where ${\bf m_{i}}$ and ${\bf m_{j}}$ are the magnetic moments at
these sites $i$ and $j$, respectively Katsura  _et al._ (2005); Rahmedov  _et
al._ (2012). Note that spin-current models have been proposed to be the origin
of magnetic cycloids in BFO Rahmedov  _et al._ (2012); Fishman (2018). Note
also that the $\mathbf{D_{1}}$ vectors of $R3m$ and $R\bar{3}c$ phases do not
add up to that of $R3c$ phase, which implies nonlinear interactions between
polarization and AFD motions in the determination of DM vectors in the $R3c$
state of BFO.

###  III.2 Second nearest neighbor coupling $\mathcal{J}_{2}$

Table 2: Calculated symmetric exchange parameters and DM interactions for the second nearest neighbor Fe-Fe pairs. $J_{2}$ and $D_{2}$ for pairs along $[$1$\bar{1}$0$]$ ($[$110$]$, respectively) directions are marked with superscript 1 (2, respectively). These parameters take into account spin-orbit interactions. (unit: meV) $[$1$\bar{1}$0$]$ |  $J_{2,xx}^{1}$ |  $J_{2,yy}^{1}$ |  $J_{2,zz}^{1}$ |  $J_{2}^{1}$
---|---|---|---|---
$R3c$ |  0.192 |  0.193 |  0.194 |  0.193
$R3m$ |  0.338 |  0.338 |  0.338 |  0.338
$R\bar{3}c$ |  0.049 |  0.048 |  0.049 |  0.049
$[$1$\bar{1}$0$]$ |  $D_{2,x}^{1}$ |  $D_{2,y}^{1}$ |  $D_{2,z}^{1}$ |  $D_{2}^{1}$
$R3c$ |  0.001 |  0.002 |  0.021 |  0.021
$R3m$ |  0.007 |  0.007 |  0.039 |  0.040
$R\bar{3}c$ |  0 |  0 |  0 |  0
$[$110$]$ |  $J_{2,xx}^{2}$ |  $J_{2,yy}^{2}$ |  $J_{2,zz}^{2}$ |  $J_{2}^{2}$
$R3c$ |  0.003 |  0.002 |  0.004 |  0.003
$R3m$ |  -0.105 |  -0.105 |  -0.102 |  -0.104
$R\bar{3}c$ |  0.150 |  0.150 |  0.150 |  0.150
$[$110$]$ |  $D_{2,x}^{2}$ |  $D_{2,y}^{2}$ |  $D_{2,z}^{2}$ |  $D_{2}^{2}$
$R3c$ |  0.000 |  -0.002 |  0.004 |  0.005
$R3m$ |  0.000 |  0.000 |  0.000 |  0.001
$R\bar{3}c$ |  0 |  0 |  0 |  0
Table 3: Calculated isotropic exchange parameters for the second nearest neighbor Fe-Fe pairs with different structures (lattices and atomic patterns). $J_{2}^{1}$ is for Fe-Fe pairs that are along $[$1$\bar{1}$0$]$ directions that are perpendicular to the polarization direction, while $J_{2}^{2}$ is for Fe-Fe pairs that are along $[$110$]$ directions. These parameters are calculated at a collinear level. Struct. |  Distor. |  |  $J_{2}$ |  Distance
---|---|---|---|---
involved |  |  (meV) |  (Å)
Cubic($Pm\bar{3}m$) |  - |  $J_{2}^{1}$,$[$1$\bar{1}$0$]$ |  0.48 |  5.56
$J_{2}^{2}$,$[$110$]$ |  0.48 |  5.56
Rhom.($R3c$) |  FE,AFD |  $J_{2}^{1}$,$[$1$\bar{1}$0$]$ |  0.35 |  5.55
$J_{2}^{2}$,$[$110$]$ |  0.25 |  5.58
Cubic($R3c$) |  FE,AFD |  $J_{2}^{1}$,$[$1$\bar{1}$0$]$ |  0.35 |  5.56
$J_{2}^{2}$,$[$110$]$ |  0.25 |  5.56
Rhom.($Pm\bar{3}m$) |  - |  $J_{2}^{1}$,$[$1$\bar{1}$0$]$ |  0.48 |  5.55
$J_{2}^{2}$,$[$110$]$ |  0.48 |  5.58
Cubic($R3m$) |  FE |  $J_{2}^{1}$,$[$1$\bar{1}$0$]$ |  0.55 |  5.56
$J_{2}^{2}$,$[$110$]$ |  0.28 |  5.56
Cubic($R\bar{3}c$) |  AFD |  $J_{2}^{1}$,$[$1$\bar{1}$0$]$ |  0.31 |  5.56
$J_{2}^{2}$,$[$110$]$ |  0.39 |  5.56

We now look at the second-nearest neighbor couplings. It is found that SE
couplings are nearly isotropic for both pairs along [1$\bar{1}$0] and [110],
since the differences between the $J_{2,\alpha\alpha}$’s (with $\alpha$ = $x$,
$y$ and $z$) are no more than 0.002 meV for both the [1$\bar{1}$0] and [110]
directions, as shown in Table II. The averaged SE coupling for pairs along
[1$\bar{1}$0] yields $J_{2}^{1}$ = 0.193 meV. Such value is very close to the
0.2 meV that is estimated from inelastic neutron scatteringFishman (2018);
Jeong  _et al._ (2012); Matsuda  _et al._ (2012); Xu  _et al._ (2012). On the
other hand, the counterpart interactions for pairs along [110] yield minute
value of $J_{2}^{2}$ $\simeq$ 0.003 meV. Such contrasts between $J_{2}^{1}$
and $J_{2}^{2}$, as well as the nearly vanishing value of $J_{2}^{2}$, are
reported here for the first time, to the best of our knowledge.

Further calculations are performed to determine whether such differences
result from the different Fe-Fe distances, FE displacements and/or AFD
motions. For simplicity, calculations without SOC (that is, we assume spins
being colinearly aligned) are performed, with the outputs being shown in Table
III, for that determination. (Note that the calculations without SOC are
purely for determining the effects of FE displacements and AFD motions and the
resulted $J_{2}$ values may differ from those with SOC.) We first check the
$J_{2}^{1}$ and $J_{2}^{2}$ coefficients for the following two phases: (i) the
cubic $Pm\bar{3}m$ phase, for which Fe-Fe pairs along $[$1$\bar{1}$0$]$ and
$[$110$]$ have the same distance and that yields the same coupling strength as
$J_{2}^{1}$ = $J_{2}^{2}$ = 0.48 meV; and (ii) the rhombohedral $R3c$ phase,
for which Fe-Fe pairs along $[$1$\bar{1}$0$]$ have shorter distance than those
along $[$110$]$, which results in different coupling strength as $J_{2}^{1}$ =
0.35 meV while $J_{2}^{2}$ = 0.25 meV. Moreover, if the internal atomic
positions retain their $R3c$ values while the lattice vectors are changed to
those of the cubic structure, the distances of Fe-Fe pairs along
$[$1$\bar{1}$0$]$ and $[$110$]$ become identical, but the coupling strengths
remain different as $J_{2}^{1}$ = 0.35 meV while $J_{2}^{2}$ = 0.25 meV.
Furthermore, if we force the internal atomic pattern to be that of the
$Pm\bar{3}m$ state while the lattice vectors are changed to those of the
rhombohedral $R3c$ ground state, the distances of Fe-Fe pairs along
$[$1$\bar{1}$0$]$ and $[$110$]$ become different again, but the coupling
strengths $J_{2}^{1}$ and $J_{2}^{2}$ turn out to be the same with the
precision up to 0.01 meV. The comparison among such cases with modified and
unmodified lattice shapes clearly demonstrates that the difference in
$J_{2}^{1}$ and $J_{2}^{2}$ is not related to the different distances (0.02 Å)
of Fe-Fe pairs, but rather if there is a polarization and/or oxygen octahedral
tilting axis in the considered state and if the considered second-nearest
neighbor direction is perpendicular or not to such polarization and/or oxygen
octahedral tilting axis. To investigate the separate effects of FE
displacements and AFD on second-nearest-neighbor couplings, we further checked
two other cases that retain the $R3m$ and $R\bar{3}c$ atomic patterns,
respectively, but with lattice vectors being those of a cubic phase. As also
shown in Table III and with respect to the situation for which both lattice
and atomic displacements are those of a cubic state (and for which $J_{2}^{1}$
= $J_{2}^{2}$ = 0.48 meV), (i) the first other case (i.e., cubic for lattice
and $R3m$ for atomic positions) enhances the couplings among the pairs that
are perpendicular to the [111] direction of polarization with $J_{2}^{1}$ =
0.55 meV, while suppressing the couplings among the pairs that are not
perpendicular to the [111] direction of polarization with $J_{2}^{2}$ = 0.28
meV; and (ii) the second other case (namely, cubic for lattice and $R\bar{3}c$
phase for atomic displacements) suppresses both types of couplings as
$J_{2}^{1}$ = 0.31 meV and $J_{2}^{2}$ = 0.39 meV. These results for these
last two cases also imply that the difference in $J_{2}^{1}$ and $J_{2}^{2}$
in the $R3c$ ground state arises from both FE and AFD displacements (and their
interactions).

Moreover, the SE couplings of second nearest neighbors in $R3m$ and
$R\bar{3}c$ phases are also found to be rather isotropic, as the corresponding
$J_{2,\alpha\alpha}$ ($\alpha=x,y$ and $z$) has the same components along
different directions, as well as that the off-diagonal components of
$\mathcal{J}_{2}$ are all smaller than 0.001 meV (not shown here). As shown in
Table II, it yields an averaged $J_{2}^{1}$ = 0.338 meV in the $R3m$ phase and
an averaged $J_{2}^{1}$ = 0.049 meV in the $R\bar{3}c$ phase for Fe-Fe pairs
along $[$1$\bar{1}$0$]$. Such two quantities work together and lead to the
medium $J_{2}^{1}$ = 0.193 meV in the $R3c$ phase. Furthermore, for Fe-Fe
pairs along $[$110$]$, $R\bar{3}c$ phase has $J_{2}^{2}$ = 0.150 meV, while
$R3m$ surprisingly has $J_{2}^{2}$ = -0.104 meV, which is ferromagnetic in
nature. Such results therefore indicate that the nearly vanishing $J_{2}^{2}$
in $R3c$ phase results from the cancellation between FE displacements and AFD.
Additionally, the facts that the diagonal elements of $J_{1}$, $J_{2}^{1}$ and
$J_{2}^{2}$ are all different when going from $R3c$ to $R3m$ or $R\bar{3}c$ is
consistent with the total energy of the effective Hamiltonian of Refs.
Rahmedov  _et al._ (2012); Kornev  _et al._ (2007) indicating that both FE and
AFD distortions affect the magnetic exchange interactions (note that a recent
study on an hexagonal phase of BFO indicates that complex isotropic
interactions can also lead to long period magnetic structure through
frustrationXu  _et al._ (2017).)

Furthermore, the DM vector between second nearest neighbors is found to nearly
vanish for $<$110$>$ pairs, while being non-negligible and lying nearly along
the $<$001$>$ direction for Fe-Fe pairs being oriented along the
$<$1$\bar{1}$0$>$ directions. In fact and as shown in Table II, such latter DM
is “only” about 6 times smaller than the DM interaction of first nearest
neighbors, and mostly originates solely from FE displacements, since the
inversion centers between second nearest neighbor Fe-Fe pairs in $R\bar{3}c$
prevent the presence of DM interactionMoriya (1960). Such facts are consistent
with a spin-current model involving polarization and magnetic moments of
second-nearest neighbors (in addition to those of first-nearest neighbors), as
done in Refs. Rahmedov  _et al._ (2012); Xu  _et al._ (2018a); Fishman  _et
al._ (2013). However, it is also worthwhile to realize that a spin-current
model for the $[$1$\bar{1}$0$]$ pair provides an energy of the form
$C_{2}({\bf P}\times{\bf e_{ij}})\cdot({\bf m_{i}}\times{\bf m_{j}})$, where
$C_{2}$ is a material-dependent parameter and where ${\bf e_{ij}}$ is the unit
vector along the $[$1$\bar{1}$0$]$ direction, which consequently should give a
${\bf D_{2}^{1}}$ DM vector along the $[\bar{1}\bar{1}2]$ direction and thus
contrasts with the nearly [001] direction found by the DFT calculations and
reported in Table III. As a result, the DFT ${\bf D_{2}^{1}}$ vector contains
effects going beyond the sole spin-current model for second-nearest neighbor
interactions (note, however, that the projection of ${\bf D_{2}^{1}}$ of the
$R3c$ phase into the $[\bar{1}\bar{1}2]$ direction gives a scalar that has a
strength of about 76% of the magnitude of ${\bf D_{2}^{1}}$, implying that
these additional effects are relatively small in comparison with those due the
spin-current model).

###  III.3 Single ion anisotropy $\mathcal{A}$

Table 4: Calculated SIA, as well as the easy axis or easy plane. Note that $3\Delta$ is the total effect of SIA, which indicates the energy difference between spins being along the [111] direction and within the (111) plane. (unit: $\mu$eV) |  $R3c$ |  $R3m$ |  $R\bar{3}c$
---|---|---|---
$\Delta$ |  -2 |  -25 |  19
$3\Delta$ |  -6 |  -75 |  57
Easy axis/plane |  [111] |  [111] |  (111)

As we have analyzed in the method part, the point group symmetry of $R3c$,
$R3m$ and $R\bar{3}c$ requires that the SIA either prefers the [111] direction
or the (111) plane. The sign and magnitude of $3\Delta$ thus defines the total
effect of SIA, which is the energy difference between local moment of one Fe
ion being along the [111] direction and within the (111) plane. As shown in
Table IV, $3\Delta$ = -6 $\mu$eV for $R3c$ phase, which indicates a weak
preference for the [111] direction. Such small value (which is, e.g., about 21
times smaller than the magnitude of the DM vector for first nearest neighbors)
is in good agreement with the experimental value of -6.8 $\mu$eVMatsuda  _et
al._ (2012) and also agrees well with the estimated value of -4 $\mu$eV from
combining different experiments and simulations Fishman (2018); Matsuda  _et
al._ (2012); Zalesskii  _et al._ (2000); Fishman  _et al._ (2013); Jeong  _et
al._ (2014); Ohoyama  _et al._ (2011); Nagel  _et al._ (2013); de Sousa  _et
al._ (2013), as well as being consistent with the neglect of SIA in effective
Hamiltonians of BFO Rahmedov  _et al._ (2012); Kornev  _et al._ (2007). Such
good agreements further attests the accuracy of our presently used four-state
method, as other numerical methods either underestimate SIA to -1.3
$\mu$eVWeingart  _et al._ (2012) or overestimate it to -11 $\mu$eVChen  _et
al._ (2018). Moreover, $3\Delta$ is found to be -75 $\mu$eV for the $R3m$
phase, therefore demonstrating that FE displacements generate an easy axis
along the [111] direction. In contrast, $3\Delta$ = 57 $\mu$eV for the
$R\bar{3}c$ phase, implying that AFD motions favor an easy (111) plane. The FE
displacements and AFD motions both have rather strong effects in determining
the SIA, as evidenced by the fact that $3\Delta$ in $R3m$ and $R\bar{3}c$
phases are an order of magnitude larger than that in the $R3c$ phase.
Interestingly, it is the competition between those two opposite effects that
results in the small SIA of the $R3c$ phase.

###  III.4 Monte-Carlo simulations

Figure 1: Magnetic properties predicted from MC simulations. Panel (a) shows
the specific heat as a function of temperature. The inset of Panel (a) shows
the dependence of the AFM Néel vector $\mathbf{L}$ on temperature, which
further emphasizes a paramagnetic-to-AFM transition taking place at 603 K;
Panel (b) displays the energy per Fe ion with respect to the period of
$[1\bar{1}0]$ and [001] cycloids; Panel (c) is the energy per Fe ion with
respect to the period of the $[1\bar{1}0]$ cycloid, using selected magnetic
parameters; and Panel (d) demonstrates the tilting angles at different
phases/positions along the propagation direction of the $[1\bar{1}0]$ cycloid.
The direction notations above the horizontal axis in Panel (d) mark the
approximate directions that the magnetic moments are parallel to. Note that
the energy of the collinear G-type AFM state is set to be energy reference
(zero) in both Panels (b) and (c).

MC simulations, using the aforementioned DFT-determined parameters and
Hamiltonian of Eq. (1), are first performed on a 12$\times$12$\times$12
supercell, therefore containing 1728 Fe atoms. As shown in Fig. 1(a), the
specific heat-versus-temperature curve shows a clear peak at 603 K, which is
indicative of a magnetic transition. We further define the AFM Néel vector
$\mathbf{L}=\frac{1}{2}|\mathbf{S_{1}}$-$\mathbf{S_{2}}|$ as the difference
between spins of the two sublattices that are represented by the two Fe sites
in the primitive cell. As shown in the inset of Fig. 1(a), the AFM Néel vector
$\mathbf{L}$ reaches the saturated value of about 2.5, showing that such
transition is from paramagnetic to the dominant G-type AFM phase. Further
analysis indicates that such G-type AFM phase in the 12$\times$12$\times$12
supercell is associated with a canted weak ferromagnetism of 0.025 $\mu$B/Fe.
The presently predicted Néel temperature $T_{N}$ = 603 K agrees rather well
with the measured value of about 643 KMoreau  _et al._ (1971); Blaauw and Van
der Woude (1973), which attests the accuracy of our magnetic parameters, as
well as the MC simulations.

The simulations on small cells (primitive cell or 2$\times$2$\times$2
supercell) are also performed, which predict not only the dominant collinear
G-type AFM configuration, but also a canting moment that further lowers the
energy by 0.09 meV/Fe, as shown in Fig. 1(b). Such canting moment results from
the $\mathbf{D_{1}^{b}}$ parameter, which originates from the oxygen
octahedral tiltings among first-nearest neighbors. The resulting magnetization
in the 2$\times$2$\times$2 supercell is determined to be 0.031 $\mu$B/Fe
(corresponding to an canting angle of 0.36∘), which agrees very well with the
value of 0.027 $\mu$B/Fe reported in previous MC effective Hamiltonian-based
simulationsAlbrecht  _et al._ (2010) and the value $\approx$0.02 $\mu$B/Fe of
the measured weak ferromagnetism Wardecki  _et al._ (2008).

We have also explored the possibility of stabilizing a spin spiral in the
$[-110]$ direction. For that we have used $\sqrt{2}n\times\sqrt{2}\times 2$
($n$ = 2, 3,…, 240) supercells, containing 4$n$ Fe ions and with its first
axis being along the [1$\bar{1}$0] direction, to determine the period of the
cycloid state along that direction. It is found that the [1$\bar{1}$0] cycloid
phase becomes lower in energy than the canted G-type AFM state, when the
cycloid period is longer than 47 nm. The minimum in the energy-versus-period
curve further indicates that the cycloid period is predicted to be $\lambda$ =
83 nm, which is slightly larger but of the same order of magnitude than the
measured 62 nm cycloidal period Sosnowska  _et al._ (1982). Note that, in
order to obtain the measured period (62 $\pm$ 3 nm), one can, for instance,
increase the magnitude of $\mathbf{D_{1}^{a}}$ from 0.102 to 0.184 meV, or
slightly increase the strength of $D_{2}^{1}$ from 0.021 to 0.032 meV and that
of $D_{2}^{2}$ from 0.005 to 0.008 meV (note also that using all parameters
directly obtained from DFT gives a critical magnetic field (aligned along the
$[11\bar{2}]$ direction) of 5.4 T associated with the magnetic-field induced
transition from the [1$\bar{1}$0] cycloid phase to canted G-type AFM state,
while increasing $\mathbf{D_{1}^{a}}$ to 0.184 meV provides a critical field
of 18.4 T – which is very close to the measured value 18 T Tokunaga  _et al._
(2010). Alternatively, if $D_{2}^{1}$ is increased to 0.032 meV and
$D_{2}^{2}$ to 0.008 meV, the critical field yields 7.1 T. It therefore
appears that having the best comparisons with different experimental data
require the choice of $\mathbf{D_{1}^{a}}$ to be 0.184 meV.) Furthermore, the
[001] cycloid is also investigated to compare with the [1$\bar{1}$0] cycloid.
It is found that (i) the [001] cycloid always has slightly higher energy than
the [1$\bar{1}$0] cycloid in all investigated range and (ii) its energy has a
minimum at $\lambda$ = 102 nm which is even lower than the energies of the
pure G-AFM state and of the spin-canted G-AFM structure, as shown in Fig.
1(b). Our predictions that the [1$\bar{1}$0] cycloid is the ground state and
that the [001] cycloid can be very close in energy is fully consistent with a
recent study using spin current model involving first and second nearest
neighbors Xu  _et al._ (2018a).

We now further look at, and report, the effects of individual magnetic
parameters in determining the stability of the magnetic configurations. (1)
The dominant isotropic first nearest neighbor magnetic exchange interaction
$J_{1}$ favors the collinear G-type AFM. The isotropic second nearest neighbor
magnetic exchange interaction parameter $J_{2}$, favors also an AFM coupling.
Therefore, $J_{1}$ and $J_{2}$ compete with each other and disfavor the
stabilization of a collinear G-type magnetic state. (2) Considering
$J_{1,\alpha\alpha}$, $J_{2,\alpha\alpha}^{1}$ and $J_{2,\alpha\alpha}^{2}$
($\alpha$ = $x$, $y$ and $z$) favors a collinear AFM within the (111) plane.
Such (111) easy plane is determined through a weak competition among pairs
along different directions. Specifically, Fe-Fe pairs along [100] ([010] and
[001], respectively) direction prefer (100) ((010) and (001), respectively)
plane, which lead to an overall effect in favor of the (111) plane. Such
competition/frustration effect is similar to the determination of the easy
axis in CrI3 and CrGeTe3 systemsXu  _et al._ (2018b). (3) The SIA favors an
easy axis along the [111] direction but the small value of 3$\Delta$ = -6
$\mu$eV is scarcely influencing magnetic properties. Specifically, when the
SIA is turned off in the MC simulations, the weakly canted G-type AFM remains
the ground state in small cells and the [1$\bar{1}$0] cycloid state remains
unchanged (aside a small increase of 1 nm of its period). Such results further
validate the neglect of SIA in effective Hamiltonians of BFO in previous
worksRahmedov  _et al._ (2012); Kornev  _et al._ (2007). (4) The DM
interactions, including $\mathbf{D_{1}^{a}}$, $\mathbf{D_{2}^{1}}$ and
$\mathbf{D_{2}^{2}}$, all contribute to generate a cycloid. Such effect is
evidenced by the facts that (i) if only isotropic $J_{1}$ and
$\mathbf{D_{1}^{a}}$ are used (all other parameters are set to be zero), it
results in a [1$\bar{1}$0] cycloid with a period of $\lambda\approx$ 122 nm;
while (ii) if $\mathbf{D_{2}}$ is also incorporated, it further stabilizes the
[1$\bar{1}$0] cycloid (by decreasing its energy) and consequently shortens the
period to $\lambda\approx$ 89 nm, as shown in Fig. 1(c). (5) The DM
interaction $\mathbf{D_{1}^{b}}$ creates spin canting in the (111) plane for
the nearest neighbor moments that have components in the (111) plane. As a
result, for a small 2$\times$2$\times$2 supercell, it leads to a homogenous
canting angle $\tau$ with the aforementioned value of 0.36∘ for the spin-
canted G-type AFM configuration. For the [1$\bar{1}$0] cycloid, there is no
canting when magnetic moments are along the [111] or [$\bar{1}\bar{1}\bar{1}$]
directions and the canting angle reaches a maximum magnitude of 0.36∘ when
moments are near the [1$\bar{1}$0] or [$\bar{1}$10] directions, as shown in
Fig. 1(d). Such modulated canting corresponds to a spin-density wave that is
formed by components of magnetic moments that are away from the plane spanned
by the [111] polarization direction and the [1$\bar{1}$0] propagation
direction, and that has been experimentally seen in Ref. Ramazanoglu  _et al._
(2011a). The maximal $|\tau|$ = 0.36∘ agrees well with the estimated 0.3∘ and
1∘values provided in Ref. Fishman (2018).

##  IV Conclusion

To conclude, the magnetic interaction parameters of multiferroic BiFeO3 are
obtained using first-principles calculations, in combination with the four-
state energy mapping method. We explicitly considered symmetric exchange
couplings (i.e., $J_{xx}$, $J_{yy}$, $J_{zz}$), DM interactions up to the
second nearest neighbor (for the first time, to the best of our knowledge), as
well as the SIA. MC simulations with those parameters successfully reproduce,
and explain, the energy hierarchy between the ground state and excited states.
The resulting [1$\bar{1}$0] cycloid has a period of 83 nm, which is in
reasonable agreement with the value of 62 nm measured in experiments. We also
predict a magnetic cycloid propagating along a $<$100$>$ direction which has a
low energy, and may thus appear in some future experiments when varying
external parameters. We are thus confident that the present work is of
interest to the scientific community, in general, and can be used as basis for
future phenomenological or ab-initio-based simulations, in particular.

###### Acknowledgements.

We thank Hongjun Xiang for useful discussion. C.X. thanks the financial
support of the Department of Energy, Office of Basic Energy Sciences, under
contract ER-46612. B.X. acknowledges funding from Air Force Office of
Scientific Research under Grant No. FA9550-16-1-0065, and L.B. thanks DARPA
Grant No. HR0011727183-D18AP00010 (TEE Program). B.D. thanks the financial
support of the Alexander von Humboldt Foundation and the Transregional
Collaborative Research Center (SFB/TRR) SPIN+X. We acknowledge the Arkansas
High Performance Computing Center (AHPCC) University of Arkansas, for using
their computing facilities.

## References

  * Sando  _et al._ (2013) D. Sando, A. Agbelele, D. Rahmedov, J. Liu, P. Rovillain, C. Toulouse, I. Infante, A. Pyatakov, S. Fusil, E. Jacquet, C. Carrétéro, C. Deranlot, S. Lisenkov, D. Wang, J. Le Breton, M. Cazayous, A. Sacuto, J. Juraszek, A. Zvezdin, L. Bellaiche, B. Dkhil, A. Barthelemy, and M. Bibes, Nature materials 12, 641 (2013).
  * Albrecht  _et al._ (2010) D. Albrecht, S. Lisenkov, W. Ren, D. Rahmedov, I. A. Kornev, and L. Bellaiche, Physical Review B 81, 140401 (2010).
  * Popov  _et al._ (1993) Y. F. Popov, A. Zvezdin, G. Vorob’Ev, A. Kadomtseva, V. Murashev, and D. Rakov, ZhETF Pisma Redaktsiiu 57, 65 (1993).
  * Popov  _et al._ (1994) Y. F. Popov, A. Kadomtseva, G. Vorob’Ev, and A. Zvezdin, Ferroelectrics 162, 135 (1994).
  * Tokunaga  _et al._ (2010) M. Tokunaga, M. Azuma, and Y. Shimakawa, Journal of the Physical Society of Japan 79, 064713 (2010).
  * Agbelele  _et al._ (2017) A. Agbelele, D. Sando, C. Toulouse, C. Paillard, R. Johnson, R. Rüffer, A. Popkov, C. Carrétéro, P. Rovillain, J.-M. Le Breton, _et al._ , Advanced Materials 29, 1602327 (2017).
  * Rovillain  _et al._ (2010) P. Rovillain, R. De Sousa, Y. Gallais, A. Sacuto, M. Méasson, D. Colson, A. Forget, M. Bibes, A. Barthélémy, and M. Cazayous, Nature materials 9, 975 (2010).
  * Popkov  _et al._ (2015) A. Popkov, N. Kulagin, S. Soloviov, K. Sukmanova, Z. Gareeva, and A. Zvezdin, Physical Review B 92, 140414 (2015).
  * Sosnowska  _et al._ (2002) I. Sosnowska, W. Schäfer, W. Kockelmann, K. Andersen, and I. Troyanchuk, Applied Physics A 74, s1040 (2002).
  * Buhot  _et al._ (2015) J. Buhot, C. Toulouse, Y. Gallais, A. Sacuto, R. De Sousa, D. Wang, L. Bellaiche, M. Bibes, A. Barthélémy, A. Forget, _et al._ , Physical review letters 115, 267204 (2015).
  * Rahmedov  _et al._ (2012) D. Rahmedov, D. Wang, J. Íñiguez, and L. Bellaiche, Physical review letters 109, 037207 (2012).
  * Xu  _et al._ (2018a) B. Xu, B. Dupé, C. Xu, H. Xiang, and L. Bellaiche, Phys. Rev. B 98, 184420 (2018a).
  * Katsura  _et al._ (2005) H. Katsura, N. Nagaosa, and A. V. Balatsky, Physical review letters 95, 057205 (2005).
  * Raeliarijaona  _et al._ (2013) A. Raeliarijaona, S. Singh, H. Fu, and L. Bellaiche, Physical review letters 110, 137205 (2013).
  * de Sousa and Moore (2008) R. de Sousa and J. E. Moore, Applied Physics Letters 92, 022514 (2008).
  * de Sousa  _et al._ (2013) R. de Sousa, M. Allen, and M. Cazayous, Physical review letters 110, 267202 (2013).
  * de Sousa (2013) R. de Sousa, in  _Spintronics VI_ , Vol. 8813 (International Society for Optics and Photonics, 2013) p. 88131L.
  * Fishman  _et al._ (2012) R. S. Fishman, N. Furukawa, J. T. Haraldsen, M. Matsuda, and S. Miyahara, Physical Review B 86, 220402 (2012).
  * Fishman  _et al._ (2013) R. S. Fishman, J. T. Haraldsen, N. Furukawa, and S. Miyahara, Physical Review B 87, 134416 (2013).
  * Fishman (2018) R. S. Fishman, Physica B: Condensed Matter 536, 115 (2018).
  * Dzyaloshinsky (1958) I. Dzyaloshinsky, Journal of Physics and Chemistry of Solids 4, 241 (1958).
  * Moriya (1960) T. Moriya, Physical Review 120, 91 (1960).
  * Wang  _et al._ (2003) J. Wang, J. Neaton, H. Zheng, V. Nagarajan, S. Ogale, B. Liu, D. Viehland, V. Vaithyanathan, D. Schlom, U. Waghmare, _et al._ , science 299, 1719 (2003).
  * Diéguez  _et al._ (2011) O. Diéguez, O. González-Vázquez, J. C. Wojdeł, and J. Íñiguez, Physical Review B 83, 094105 (2011).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical Review Letters 100, 136406 (2008).
  * Xu  _et al._ (2014) C. Xu, Y. Yang, S. Wang, W. Duan, B. Gu, and L. Bellaiche, Physical Review B 89, 205122 (2014).
  * Xiang  _et al._ (2013) H. Xiang, C. Lee, H.-J. Koo, X. Gong, and M.-H. Whangbo, Dalton Transactions 42, 823 (2013).
  * Xiang  _et al._ (2011) H. Xiang, E. Kan, S.-H. Wei, M.-H. Whangbo, and X. Gong, Physical Review B 84, 224429 (2011).
  * Miyatake  _et al._ (1986) Y. Miyatake, M. Yamamoto, J. Kim, M. Toyonaga, and O. Nagai, Journal of Physics C: Solid State Physics 19, 2539 (1986).
  * Kubel and Schmid (1990) F. Kubel and H. Schmid, Acta crystallographica section B 46, 698 (1990).
  * Matsuda  _et al._ (2012) M. Matsuda, R. S. Fishman, T. Hong, C. Lee, T. Ushiyama, Y. Yanagisawa, Y. Tomioka, and T. Ito, Physical review letters 109, 067205 (2012).
  * Jeong  _et al._ (2012) J. Jeong, E. Goremychkin, T. Guidi, K. Nakajima, G. S. Jeon, S.-A. Kim, S. Furukawa, Y. B. Kim, S. Lee, V. Kiryukhin, _et al._ , Physical review letters 108, 077202 (2012).
  * Xu  _et al._ (2012) Z. Xu, J. Wen, T. Berlijn, P. M. Gehring, C. Stock, M. B. Stone, W. Ku, G. Gu, S. M. Shapiro, R. Birgeneau, _et al._ , Physical Review B 86, 174419 (2012).
  * Feng and Xiang (2016) J. Feng and H. Xiang, Physical Review B 93, 174416 (2016).
  * Chen  _et al._ (2018) Z. Chen, Z. Chen, C.-Y. Kuo, Y. Tang, L. R. Dedon, Q. Li, L. Zhang, C. Klewe, Y.-L. Huang, B. Prasad, _et al._ , Nature communications 9, 3764 (2018).
  * Ramazanoglu  _et al._ (2011a) M. Ramazanoglu, M. Laver, I. W Ratcliff, S. Watson, W. Chen, A. Jackson, K. Kothapalli, S. Lee, S.-W. Cheong, and V. Kiryukhin, Physical review letters 107, 207206 (2011a).
  * Bellaiche  _et al._ (2012) L. Bellaiche, Z. Gui, and I. A. Kornev, Journal of Physics: Condensed Matter 24, 312201 (2012).
  * Ederer and Spaldin (2005) C. Ederer and N. A. Spaldin, Physical Review B 71, 060401 (2005).
  * Sosnowska and Zvezdin (1995) I. Sosnowska and A. Zvezdin, Journal of magnetism and magnetic materials 140, 167 (1995).
  * Ramazanoglu  _et al._ (2011b) M. Ramazanoglu, M. Laver, I. W Ratcliff, S. Watson, W. Chen, A. Jackson, K. Kothapalli, S. Lee, S.-W. Cheong, and V. Kiryukhin, Physical review letters 107, 207206 (2011b).
  * Ruette  _et al._ (2004) B. Ruette, S. Zvyagin, A. P. Pyatakov, A. Bush, J. Li, V. Belotelov, A. Zvezdin, and D. Viehland, Physical Review B 69, 064114 (2004).
  * Jeong  _et al._ (2014) J. Jeong, M. D. Le, P. Bourges, S. Petit, S. Furukawa, S.-A. Kim, S. Lee, S. Cheong, and J.-G. Park, Physical review letters 113, 107202 (2014).
  * Kornev  _et al._ (2007) I. A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, and L. Bellaiche, Physical Review Letters 99, 227602 (2007).
  * Xu  _et al._ (2017) C. Xu, H. Xiang, and L. Bellaiche, Advanced Electronic Materials 3, 1700332 (2017).
  * Zalesskii  _et al._ (2000) A. Zalesskii, A. Zvezdin, A. Frolov, and A. Bush, Journal of Experimental and Theoretical Physics Letters 71, 465 (2000).
  * Ohoyama  _et al._ (2011) K. Ohoyama, S. Lee, S. Yoshii, Y. Narumi, T. Morioka, H. Nojiri, G. Sang Jeon, S.-W. Cheong, and J.-G. Park, Journal of the Physical Society of Japan 80, 125001 (2011).
  * Nagel  _et al._ (2013) U. Nagel, R. S. Fishman, T. Katuwal, H. Engelkamp, D. Talbayev, H. T. Yi, S.-W. Cheong, and T. Room, Physical review letters 110, 257201 (2013).
  * Weingart  _et al._ (2012) C. Weingart, N. Spaldin, and E. Bousquet, Physical Review B 86, 094413 (2012).
  * Moreau  _et al._ (1971) J.-M. Moreau, C. Michel, R. Gerson, and W. J. James, Journal of Physics and Chemistry of Solids 32, 1315 (1971).
  * Blaauw and Van der Woude (1973) C. Blaauw and F. Van der Woude, Journal of Physics C: Solid State Physics 6, 1422 (1973).
  * Wardecki  _et al._ (2008) D. Wardecki, R. Przenioslo, I. Sosnowska, Y. Skourski, and M. Loewenhaupt, Journal of the Physical Society of Japan 77, 103709 (2008).
  * Sosnowska  _et al._ (1982) I. Sosnowska, T. P. Neumaier, and E. Steichele, Journal of Physics C: Solid State Physics 15, 4835 (1982).
  * Xu  _et al._ (2018b) C. Xu, J. Feng, H. Xiang, and L. Bellaiche, npj Computational Materials 4, 57 (2018b).


# ICS-CTM2: Industrial Control System Cybersecurity Testbed Maturity Model

Souradeep Bhattacharya, Burhan Hyder, Manimaran Govindarasu  Department of
Electrical and Computer Engineering
Iowa State University, Ames, IA 50011
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Industrial Control System (ICS) testbeds serve as a platform for evaluating
and validating control system performances, cybersecurity tools and
technologies. In order to build or enhance an ICS testbed, it is vital to have
a deeper understanding of its design specifications and characteristic
attributes. Satisfying this prerequisite involves examination and assessment
of these attributes for existing testbeds. To further increase confidence in a
testbed’s functionality, it is important to perform a comparative analysis of
its specifications with other ICS testbeds. However, at present, there is no
standardized methodology available to provide a comparative assessment of
different testbeds. In this paper, we propose a methodology for analyzing ICS
testbeds, inspired by the Cybersecurity Capability Maturity Model (C2M2). In
particular, we then define a ICS Cybersecurity Testbed Maturity Model, its
domains, and the associated maturity indicator levels. To demonstrate the
benefit of the model, we have conducted a case study analysis for several ICS
testbeds, representing different industrial sectors. Our analysis provides
deeper insights into the relative strengths and limitations of these testbeds,
together with scope for future enhancements, with respect to the domains
defined by the model.

###### Index Terms:

Industrial Control System (ICS), Cybersecurity, Testbed maturity model

##  I Introduction

In recent years, with the emergence of Industry 4.0 criterion, there has been
a considerable increase in digitization of legacy Industrial Control Systems
(ICS). This accelerated process has resulted in growing convergence between
information technologies (IT) and operational technologies (OT), thereby
increasing the inter-connectivity of different systems. Owing to advancement
in communication, information, and operational technologies, new approaches in
the industrial chain have risen enabling tremendous growth for ICS, such as
industrial internet of things, big data analytics, cloud computing, robotics
and human-machine interactions [1]. However, the paradigm shift towards
digitization and inter-connectivity exposes ICSs to new vulnerabilities and
challenges. This in turn increases the attack surfaces for adversaries to
exploit. Therefore, academia, industry and government entities are
increasingly proposing development of ICS testbeds for evaluating tools and
technologies designed to assess system performance metrics and enhance
cybersecurity capabilities.

ICSs interconnect, monitor, and control processes in a variety of industries,
such as electric power generation and distribution, transportation, water,
gas, oil, chemical production, and manufacturing. ICS consist of two major
components - the OT network, and the IT network. The OT network comprises of
hardware and software that are used to monitor and manage industrial
equipment, processes, assets, and events, e.g., Programmable Logic Controllers
(PLC), Supervisory Control and Data Acquisition (SCADA), Distributed Control
System (DCS), sensors and actuators. The IT network comprises of workstations,
servers, databases, routers and network switches, which are used to control
and manipulate the flow of information.

This paper is structured as follows. In Section-I, we provide a brief
introduction of ICS, highlighting the differences between IT and OT networks.
In Section-II, we discuss why it is crucial to develop testbeds for ICS and
what are the requirements to consider. We also refer some of the related works
in this field. Lastly, in this section we have also discussed the need for
developing a systematic testbed assessment methodology for analyzing different
ICS testbeds. In Section-III we describes the design and implementation
methodology for the proposed ICS Cybersecurity Testbed Maturity Model (ICS-
CTM2) framework. Finally, Section-IV provides a case study for ICS testbed
analysis using the proposed ICS-CTM2 framework, and also presents a
visualization of the results.

##  II Overview of ICS Testbeds

A testbed can be defined as a testing environment or a platform that simulates
real-world activities under supervision, in order to conduct experimental
evaluations. When researchers need to work with real-world ICS environments
for examining new technologies, testbeds provide a validated ecosystem
representing the real-world process.

###  II-A Why Do We Need A Testbed For ICS

Testing IT security systems is a mature, well-defined field with a plethora of
best practices, procedures and information. However, testing industrial
control systems and OT networks is a comparatively newer frontier [2].
Evolution of today’s ICS has been a causal effect of amalgamation of IT
capabilities into legacy physical systems, often replacing or augmenting
existing physical control operations. However, it also comes with various
challenges that need to be mitigated - cybersecurity being critical. As per
annual reports published by Kaspersky, it can be observed that there is a
considerable increase in vulnerabilities detected in ICS every year (19 in
2010, 189 in 2015, 415 in 2018, and 509 in 2019) [3]. Therefore, developing
custom testing environments imitating real-world solutions are a crucial
component in ICS research. In addition to vulnerability research, cyber-
physical ICS testbeds provide several applications such as impact analysis
(capturing risk posed by security events); mitigation evaluation (to increase
robustness of cyber infrastructure); cyber-physical metrics (functional tests
for system reliability and performance); data and model development (models
and datasets to facilitate more accurate analysis); security validation
(assessment of cybersecurity compliance requirements); cyber forensics
(analysis of cyber attacks); operator training and education [4].

###  II-B Requirements For Developing An ICS Testbed

Building an ICS testbed is a cost-dependent, resource-intensive, and complex
task, involving collaboration of multiple disciplines. As it is essential to
perform experimental activities and generate data from different industrial
scenarios, it is imperative to develop custom testbeds accurately replicating
the respective industrial practices.

####  II-B1 Testbed Classification

Depending on the process, protocol, and infrastructure elements involved in a
testbed, it can be classified into one of the four categories - Physical,
Simulation, Virtual, or Hybrid [5]. In a physical testbed, real hardware and
software are deployed in both the cyber and physical layers. Although, it aims
to be as close to the real system as possible providing highest fidelity,
physical testbeds are usually expensive and have a long development cycle.
Simulation testbeds are based on software simulations, whereas virtual
testbeds involve both simulations and emulation of hardware components.
Simulation and virtual testbeds provide a low-cost alternative to physical
testbeds, and are more flexible to upgrades. However, these testbeds may lack
precision and often does not provide real-time operation. The most commonly
used approach is the hybrid testbed. These testbeds are designed through a
combination of physical devices, simulation, and virtualized components.
Hybrid testbeds provides lower fidelity than physical testbeds, but allows for
lower development cost and higher scalability.

####  II-B2 Testbed Properties

Any practical testbed is required to satisfy certain clearly defined
objectives. These objectives constitute the properties of an ICS testbed which
indicate its reliability [6]. Although, it may be challenging to meet every
objective, it is vital for researchers to ascertain an effective and optimized
trade-off criteria during the design phase. All properties pertaining to an
ICS testbed can be grouped under three major inter-related properties as
follows:

Fidelity: Fidelity of a testbed can be defined as how closely and accurately
the testbed replicates a real-world ICS [6], [7]. The properties which
directly impact fidelity are: Measurability (allowing assessment of
experiments), Measurement Accuracy (degree of correctness of experimental
results), Repeatability (ability to replicate experiments), Reproducibility
(ability to produce consistent results), Usability (ability to be utilized for
defined purposes), Complexity (transparency of architecture), and Safety (safe
operation of testbed).

Scale: Scalability of a testbed can be defined as a testbed’s ability to be
expanded in functionality, architecture, and, infrastructure [4], [6]. The
properties which directly impact scale are: Cyber-Physical Integration,
Reconfigurability (ability to accommodate new design or components),
Extensibility (ability to broaden functional scope), Adaptability (ability to
support new or alternative test cases), Interoperability (ability to support
any combination of simulation or hybrid approaches), Federation (inter-
connectivity of multiple testbeds), and Critical Infrastructure.

Cost: Cost of a testbed depends on the approximate implementation cost of the
process, hardware, software, and licenses [4], [7]. The properties which
directly impact cost are: Diversity (ability to utilize products and services
from multiple vendors), Openness (ability to support open-source products and
remote access), Training & Development (for education of operators), and
Classification (Physical, Simulation, Virtual, or Hybrid).

####  II-B3 Compliance of Industrial Standards

NIST SP 800-82 Rev.2 [8] provides a guideline for ICS Security. It specifies
that the basic architecture for an ICS environment must include four core
components: Physical Process, Field Devices, Communication Architecture, and
Control Center. Therefore, adherence to standardized design attributes
contribute significantly to the reliability and trustworthiness of ICS
testbeds.

###  II-C Related Work

In this paper, we directed our efforts towards performing an extensive review
of the current state in the field of Industrial Control System (ICS) testbeds.
Among the existing literature, survey article [5] focuses on exploring the
primary objectives and application possibilities of 30 different ICS testbeds.
Another survey [7] provides a comprehensive classification of ICS testbeds
into three groups (physical, virtual, or hybrid), along with an analysis of
ICS datasets available for developing new IDS security mechanisms. [9]
provides an outline for developing a SCADA Security Testbed. Some existing
literature (e.g. [10], [11]) provides detailed description of testbed
architecture and design, implementation and experiments conducted. For
instance, [10] discusses the PowerCyber testbed at Iowa State University,
which models a smart electric grid, and [11] discusses the SMART testbed at
University of Michigan, which models a warehouse production line.

Figure 1: The Proposed ICS Cybersecurity Testbed Maturity Model (ICS-CTM2)
Framework

###  II-D Assessment of ICS Testbeds

Although survey articles provide a compilation of different testbeds, they
lack any quantitative assessment methodology which makes it difficult to
easily recognise relevant variations in the listed testbeds. Literature
pertaining to design and implementation of ICS testbeds focus only on the
individual testbeds, without any comparative analysis with other testbeds. To
the best of our knowledge, there is no standardized methodology available at
present, to provide a comparative assessment of different testbeds. Therefore,
a systematic methodology is needed which can provide a focal guide for
researchers developing future ICS testbeds or enhancing current testbeds.

##  III Proposed ICS Cybersecurity Testbed Maturity Model

The central objective of the ICS Cybersecurity Testbed Maturity Model (ICS-
CTM2) framework is to provide a standardized and systematic method for
comprehensive and comparative analysis of ICS testbed design specifications.
The ICS-CTM2 framework is inspired from the Energy Subsector-Cybersecurity
Capability Maturity Model (ES-C2M2) developed by the US energy sector and
Department of Energy (DOE) [12].

Figure 1 shows the conceptual diagram for the proposed ICS-CTM2 framework.
Table I provides a list of ICS testbeds that we have considered for the case
study analysis of ICS-CTM2 and is further described in detail under Section
IV.

ES-C2M2: The ES-C2M2 [12] is a tool which allows electric utilities and grid
operators to assess, evaluate, and improve their cybersecurity capabilities,
and also to prioritize their actions and investments for enhancing
cybersecurity. The ES-C2M2 was developed in 2012 as part of a White House
initiative led by the Department of Energy (DOE) in partnership with the
Department of Homeland Security (DHS), involving collaboration with industry
members, other Federal agencies and stakeholders. ES-C2M2 provides a toolkit
(available on request from DOE) that can be utilized by any organization to
measure and improve its cybersecurity program. Owing to its self-evaluation
methodology and industry focused application which makes ES-C2M2 an easily
scalable tool, in this paper we adapted ES-C2M2 to develop our ICS-CTM2
framework.

###  III-A Methodology

Maturity model can be defined as an organized and structured method to convey
a path of knowledge, experience, and learning. Such assessment models enable
researchers and organizations to evaluate and make improvements to their ICS
testbeds, by providing standardized benchmarks. Inspired by ES-C2M2, the ICS-
CTM2 model arises from a combination of existing industrial standards,
frameworks, and procedures adapted for developing an ICS testbed. The ICS-CTM2
architecture comprises of the following three sections, also illustrated in
Figure-1:

####  III-A1 Domains

ICS-CTM2 domains represent a structured set of parameters that include
essential information related to an ICS testbed. ICS-CTM2 framework consists
of five domains (Architecture, Fidelity, Scale, Cost, and Application) based
on which each testbed can be evaluated and its capabilities can be determined.
For each domain, the model provides a description, which summarizes the
objective of the domain along with the evaluation criterion considered within
the context of the domain. The defined objective is then assessed by
determining the extent of implementation of each evaluation criteria.

Figure 2 describes the evaluation criteria for each domain in ICS-CTM2
framework. Based on this assessment, a maturity level is assigned to each
domain for a testbed. The five domains for ICS-CTM2 model are:

TABLE I: list of selected testbeds for ICS-CTM2 case study analysis Sl.No. | Institute | Industry
---|---|---
1 | Iowa State University (ISU) [10] | Smart Grid
2 | SUTD, Singapore [13] | Water Treatment
3 | University of Michigan (U-M) [11] | Manufacturing
4 | ORNL [14] | HVAC/Cooling
5 | Univ. of Alabama, Huntsville (UAH) [15] | Gas Pipeline
6 | Univ. of Tennessee (U.T) [16] | Nuclear Plant
7 | Ohio State University (OSU) [17] | Automobile
8 | NIST [18] | Multi-scenario
Figure 2: Case Study: Maturity Indicator Levels of each ICS-CTM2 domain for
selected ICS Testbeds. Figure 3: Case Study: Radar Analysis for ICS Testbeds
(1)-(4) Figure 4: Case Study: Radar Analysis for ICS Testbeds (5)-(8)

Architecture: The Purdue Enterprise Reference Architecture (PERA) [19] is the
most widely accepted and utilized ICS network architecture model in the
industry. Since ICS consists of different equipments and protocols across
different industries and regions, PERA serves as a common reference model for
evaluating this domain of ICS-CTM2 framework. PERA was developed by Theodore
J. Williams in 1993, as a collaborative effort between Purdue University and
members of the industry. This model has been adopted by several industry
security standards such as ISA-99 (ISA/IEC 62443) [20] and is used as the
foundation for ICS IT/OT network design. PERA is segmented into five levels as
described below.
Level 0 - Field Level: Includes the physical process to be implemented.
Level 1 - Control Level: Includes field controllers and devices (PLC, RTU)
implemented for controlling the physical process.
Level 2 - Supervisory Level: Includes the supervisory level controllers and
devices (SCADA, HMI) implemented for controlling and monitoring of the field
level controllers.
Level 3 - Operations/Monitoring Level: Includes the components implemented for
information sharing between the IT and OT layer (such as Workstation, HMI,
DMZ).
Level 4 - Management/Enterprise Level: Includes implementation of the
management and planning resources at the enterprise level, such as resource
management or risk assessment.

Fidelity, Scale and, Cost: These three inter-related domains represent the
several properties of an ICS testbed, as described earlier in Section II-B.
Evaluation of the fidelity domain for ICS-CTM2 framework is based on the
testbed construction methodology adopted for implementation of three key
components - Process, Protocol, and Infrastructure (PPI). Implementation of
multiple process control scenarios and cyber-physical integration influence
the scalability of a testbed. Cost of design is determined by factors such as
physical components used, open-sourced or licensed products, training and
development expenses.

Application: We have evaluated ICS testbeds based on the practical
applications and use cases implemented [4]. The application scenarios have
been divided into four categories in the order of increasing criticality: (1)
Functionality Testing, (2) Education, (3) Cybersecurity Analysis and Research,
and (4) Development of Standards .

####  III-A2 Maturity Indicator Level (MIL)

ICS-CTM2 framework defines four Maturity Indicator Levels - MIL0 thru MIL3. A
fifth maturity level (MIL4) is reserved for future use. The MILs helps to
identify the degree of implementation in a domain. The MILs are independent to
each domain, and to earn a MIL in a given domain, a testbed must satisfy all
the evaluation criterion specified in that level. For instance, if a testbed
satisfies all criteria in MIL1 and MIL2, but not in MIL3, then the testbed
would be assigned a level of MIL2 in that domain.

Figure 5: Evaluation Criterion considered within the ICS-CTM2 framework
domains (case study performed for PowerCyber Testbed [10]).

####  III-A3 Diagnosis and Assessment:

For a testbed, each domain is assessed depending on the evaluation criterion
and a MIL is assigned. First, we analyze a testbed with respect to each of the
five domains. Analyzing a domain involves assessment of each evaluation
criterion within that domain. Once all criteria are evaluated, a MIL rating is
assigned to that domain. Similarly, all remaining domains are assessed.
Finally, the ICS-CTM2 diagnosis results are visualized though Radar Analysis
and Ring Analysis.

Radar Analysis: The radar analysis visualizes the MIL ratings obtained by a
testbed in each of the five domains. Radar analysis can be performed for a
single as well as multiple ICS testbeds, and the results can be aggregated and
displayed together for better comparison. All domains are arranged in a radial
chart and the respective maturity ratings are recorded after implementation of
ICS-CTM2 framework. Figure 3 and Figure 4 demonstrates the Radar Analysis.

Ring Analysis: The ring analysis provides a more in-depth assessment by
illustrating the implementation degree of each evaluation criterion for all
the five domains of the ICS-CTM2 framework. Ring analysis is a higher level
evaluation, hence it must be performed for each testbed individually. As such,
results of different testbeds cannot be combined for ring analysis. Figure 6
demonstrates the Ring Analysis.

Figure 6: Case Study: Ring Analysis of Evaluation Criterion for all five
domains (case study performed for PowerCyber Testbed [10]).

##  IV Case Study Analysis

This section describes a case study implemented to analyze different testbeds,
mentioned in Table I, using the ICS-CTM2 framework. In our experiment, we have
considered one ICS testbed from 8 different industrial sectors, to demonstrate
their differences. However, the ICS-CTM2 framework is not limited to this
comparison only, and can be applied to any combination of ICS testbeds as
deemed necessary.

Figure 2 provides the MIL ratings for the ICS-CTM2 domains after analyzing the
selected testbeds. Figure 3 and Figure 4 represents the Radar Analysis
comparing the performance of testbeds in each domain utilizing results
demonstrated in Figure 2. The results of the selected testbeds have been
divided into two plots for ease of illustration. As shown in the figures, from
the Radar Analysis the researchers can easily identify the maturity of the
testbed in each domain.

Figure 6 demonstrates the Ring Analysis for the ICS-CTM2 framework. We have
performed this evaluation only for the PowerCyber Testbed [10] at Iowa State
University. The evaluation criterion for each domain have been specified in
Figure 5. The number specified in the inner ring indicates the total
evaluation criterion considered at that maturity level. The number specified
in the outer ring summarizes the total evaluation criterion satisfied or
implemented (either fully, partially, or not implemented). For example, in
Figure 6 Ring Analysis, the domain Architecture has 5 evaluation criterion for
MIL1; 4 criterion are added in the next level making the total 9 for MIL2, and
12 criterion are added in the next level making the total 21 for MIL3. The
testbed satisfies all criterion for MIL1 and MIL2, but not for MIL3.
Therefore, the testbed achieves an overall rating of MIL2 in this domain.

##  V Conclusion and Future Work

The need for a standardized method to perform comparative assessment of ICS
Testbeds is quite compelling. The design goal of the ICS-CTM2 framework,
described in this paper, is to enable researchers to perform a self-evaluation
of ICS Testbeds with respect to their design specifications and characteristic
attributes, to determine their strengths and limitations. The ICS-CTM2
framework enables ICS testbeds to be analyzed on the basis of standardized
domains which serve as benchmarks. Each evaluation criterion within the
domains can be assessed to gauge the maturity in the design of a testbed. This
paper also provides a case study analysis performed on 8 ICS testbeds to
demonstrate the functionality of ICS-CTM2 framework and visualization of
results.

Future Work: The implementation methodology and the framework domains
described in this paper is being actively refined, to enhance the assessment
parameters. The future work includes: (1) a comparative assessment of a larger
number of ICS testbeds, applying ICS-CTM2 framework, based on the data about
testbed specifications collected through scientific survey, (2) the proposed
ICS-CTM2 model can be further refined by updating or expanding the domains,
evaluation criterion, or the analysis method.

##  VI Acknowledgement

This work is funded in part by the NSF CPS grant ECCS 1739969.

## References

  * [1] C. Alcaraz, “Secure Interconnection of IT-OT networks in Industry 4.0 in Critical Infrastructure Security and Resilience. Springer, 2019, pp.201–217.,” 2019.
  * [2] F. Khorrami, P. Krishnamurthy, and R. Karri, “Cybersecurity for Control Systems: A Process-Aware Perspective,” 10 2016.
  * [3] T. Menze, “Kaspersky, The State of Industrial Cybersecurity,” September 2020.
  * [4] A. Hahn, A. Ashok, S. Sridhar, and M. Govindarasu, “Cyber-Physical Security Testbeds: Architecture, Application, and Evaluation for Smart Grid,” IEEE Transactions on Smart Grid, vol. 4, pp. 847–855, 2013.
  * [5] H. Holm, M. Karresand, A. Vidström, and E. Westring, “A survey of industrial control system testbeds. In Proceedings of the Nordic Conference on Secure IT Systems, Stockholm, Sweden, 19–21 October 2015; Springer: Berlin/Heidelberg, Germany, 2015; pp. 11–26.,” vol. 9417, Springer International Publishing, 2015.
  * [6] U. P. D. Ani and J. Watson, “What Makes an Industrial Control System Security Testbed Credible and Acceptable? Towards a Design Consideration Framework,” pp. 181–190, July 2021.
  * [7] M. Conti, D. Donadel, and F. Turrin, “A Survey on Industrial Control System Testbeds and Datasets for Security Research,” CoRR, vol. abs/2102.05631, 2021.
  * [8] S. K., V. Pillitteri, S. Lightman, M. Abrams, and A. Hahn, “NIST SP 800-82 Rev. 2, Guide to Industrial Control Systems (ICS) Security,” 2015\.
  * [9] E. Christiansson and H. Luiijf, “Creating a European SCADA Security Testbed,” pp. 237–247, Springer US, 2008.
  * [10] A. Ashok, S. Krishnaswamy, and M. Govindarasu, “Powercyber: A Remotely Accessible Testbed for Cyber Physical Security of the Smart Grid,” Institute of Electrical and Electronics Engineers Inc., 12 2016.
  * [11] I. Kovalenko, M. Saez, K. Barton, and D. Tilbury, “SMART: A System-level Manufacturing and Automation Research Testbed,” 2017.
  * [12] “Department of Energy (DOE), Electricity Subsector Cybersecurity Capability Maturity Model (ES-C2M2),” 2014.
  * [13] A. P. Mathur and N. O. Tippenhauer, “SWaT: A Water Treatment Testbed for Research and Training on Ics Security,” pp. 31–36, 2016\.
  * [14] R. E. Gillen, L. A. Anderson, C. Craig, J. Johnson, A. Columbia, R. Anderson, A. Craig, and S. L. Scott, “Design and Implementation of Full-Scale Industrial Control System Testbed for Assessing Cyber-Security Defenses,” pp. 341–346, 2020.
  * [15] T. Alves, R. Das, and T. Morris, “Virtualization of Industrial Control System Testbeds for Cybersecurity,” pp. 10–14, Association for Computing Machinery, 2016.
  * [16] F. Zhang, H. A. D. E. Kodituwakku, J. W. Hines, and J. Coble, “Multilayer Data-Driven Cyber-Attack Detection System for Industrial Control Systems Based on Network, System, and Process Data,” IEEE Transactions on Industrial Informatics, vol. 15, pp. 4362–4369, 2019\.
  * [17] P. S. Oruganti, M. Appel, and Q. Ahmed, “Hardware-in-loop based Automotive Embedded Systems Cybersecurity Evaluation Testbed,” pp. 41–44, Association for Computing Machinery, Inc, March 2019.
  * [18] R. Candell, D. M. Anand, and K. Stouffer, “A Cybersecurity Testbed for Industrial Control Systems,” 2014.
  * [19] T. J. Williams, “The Purdue Enterprise Reference Architecture and Methodology (PERA),” 1994.
  * [20] ISA, ”ANSI/ISA-62443-1-1 (99.01.01)-2007, Security for Industrial Automation and Control Systems: Part 1 Terminology, Concepts, and Models. American National Standards Institute, 2007.