# Towards automated extraction and characterization of scaling regions in
dynamical systems

Varad Deshmukh  Department of Computer Science, University of Colorado,
Boulder, CO, USA  Elizabeth Bradley  Department of Computer Science,
University of Colorado, Boulder, CO, USA  Santa Fe Institute, Santa Fe, NM,
USA  Joshua Garland  Santa Fe Institute, Santa Fe, NM, USA  James D. Meiss
Department of Applied Mathematics, University of Colorado, Boulder, CO, USA

###### Abstract

Scaling regions—intervals on a graph where the dependent variable depends
linearly on the independent variable—abound in dynamical systems, notably in
calculations of invariants like the correlation dimension or a Lyapunov
exponent. In these applications, scaling regions are generally selected by
hand, a process that is subjective and often challenging due to noise,
algorithmic effects, and confirmation bias. In this paper, we propose an
automated technique for extracting and characterizing such regions. Starting
with a two-dimensional plot—e.g., the values of the correlation integral,
calculated using the Grassberger-Procaccia algorithm over a range of scales—we
create an ensemble of intervals by considering all possible combinations of
endpoints, generating a distribution of slopes from least-squares fits
weighted by the length of the fitting line and the inverse square of the fit
error. The mode of this distribution gives an estimate of the slope of the
scaling region (if it exists). The endpoints of the intervals that correspond
to the mode provide an estimate for the extent of that region. When there is
no scaling region, the distributions will be wide and the resulting error
estimates for the slope will be large. We demonstrate this method for
computations of dimension and Lyapunov exponent for several dynamical systems,
and show that it can be useful in selecting values for the parameters in time-
delay reconstructions.

††preprint: Currently in review at Chaos: An Interdisciplinary Journal of
Nonlinear Science

Lead Paragraph: Many problems in nonlinear dynamics involve identification of
regions of constant slope in plots of some quantity as a function of a length
or time scale. The slopes of these scaling regions, if they exist, allow one
to estimate quantities such as a fractal dimension or a Lyapunov exponent. In
practice, identifying scaling regions is not always straightforward. Issues
with the data and/or the algorithms often cause the linear relationship to be
valid only for a limited range in such a plot, if it exists at all. Noise,
geometry and data quality can disturb its shape in various ways, and even
create multiple scaling regions. Often the presence of a scaling region, and
the endpoints of its range, are determined by eye, Kantz and Schreiber (1997);
Bradley and Kantz (2015) a process that is subjective and not immune to
confirmation bias.Ji, Zhu, and Jiang (2011); Chen, Liu, and Zhou (2019);
Clauset, Shalizi, and Newman (2009); Zhou and Wang (2018) Worse yet, we know
of no formal results about the relationship between the width of the scaling
region and the validity of the resulting estimate; often practitioners simply
use the heuristic notion that “wider is better.” In this work, we propose an
ensemble-based method to objectively identify and quantify scaling regions,
thereby addressing some of the issues raised above.

##  I Introduction

To identify a scaling region in a plot, our method begins by generating
multiple fits using intervals of different lengths and positions along the
entire range of the data set. Penalizing the lower-quality fits—those that are
shorter or have higher error— we obtain a weighted distribution of the slopes
across all possible intervals on the curve. As we will demonstrate by example,
the slope of a scaling region that is broad and straight manifests as a
dominant mode in the weighted distribution, allowing easy identification of an
optimal slope. Moreover, the extent of the scaling region is represented by
the modes of a similar distribution that assesses the fit as a function of the
interval endpoints. As we will show, for a long, straight scaling region,
markers closer to the endpoints of the scaling region are more frequently
sampled (due to the combinatorial nature of sampling) and have a higher
weighting (due to longer lengths of the fits). Thus, the method gives the
largest reasonable scaling interval for the optimal fit along with error
bounds on the estimate as computed from the distributions.

Section II covers the details of this ensemble-based method and illustrates
its application with three examples. In §III, we demonstrate the approach on
several dynamical systems problems, starting in §III.1 with the estimation of
correlation dimension for two canonical examples. We show that our method can
accurately detect spurious scaling regions due to noise and other effects, and
that it can be useful in systematically estimating the embedding dimension for
delay reconstruction. We then demonstrate that this method is also useful for
calculating other dynamical invariants, such as the Lyapunov exponent
(§III.2). We present the conclusions in §IV.

##  II Characterizing scaling regions

A scaling law “…describe[s] the functional relationship between two physical
quantities that scale with each other over a significant
interval.”***https://www.nature.com/subjects/scaling-laws Such a scaling law
manifests as a straight segment on a two-dimensional graph of these physical
quantities, known as a scaling region. In practical situations, the data will
typically not lie exactly on a line or may do so only over a limited interval,
see e.g., the synthetic example in Fig. 1(a).

(a) $y(x)=s(x)+d(x)$

(b) Weighted slope histogram with the KD estimator

(c) Weighted distributions of interval endpoints

Figure 1: Extracting a scaling region using an ensemble-based approach. (a) A
synthetic example with a linear range that is bounded below by oscillations
and above by a flat interval. (b) Distribution of slopes (orange histogram)
generated from an ensemble of fits in different intervals of this plot,
weighted by the inverse square of the fitting error. A kernel density
estimation for the probability density is shown in blue, with its mode marked
by the black line. The horizontal axis is clipped to [4,6]. (c) Probability
distributions of the interval endpoints $[x_{l},x_{r}]$ from the ensemble,
with the same weighting as (b).

These data are obtained from a function $y(x)=s(x)+d(x)$, where $s(x)$ is
piecewise linear, with slope $5$ in the region $1\leq x\leq 9$:

$s(x)=\begin{cases}0,&x<1\\\ 5x-5,&1\leq x\leq 9\\\ 40,&9<x\end{cases}.$

This function exhibits a scaling region in the range $[1,9]$ that is bounded
by intervals with zero slope: a shape that is similar to what is often seen in
computations of dynamical invariants even for well-sampled, noise-free data.
To make the problem less trivial, we add a damped oscillation,

$d(x)=10\ e^{-0.5x}\sin(5x),$

to $s(x)$. The data set for this example is taken to be $y(x)$ on $[0,10)$
with a sampling interval $\Delta x=0.1$.

To find a scaling region in a diagram like Fig. 1(a), one would normally
choose two points—markers $x_{l}$ and $x_{r}$ for the bounds of the linear
region—then fit a line to the interval $[x_{l},x_{r}]$. A primary challenge is
to identify the appropriate interval. For Fig. 1(a), the range $7\leq x\leq 9$
appears to be linear, but should this range extend to smaller $x$? One could
certainly argue that the lower boundary could be $x_{l}=6$ instead of $7$, but
it would be harder to defend a choice of $x_{l}\leq 4$ because of the growing
oscillations.

Several approaches have been proposed for formalizing this procedure in the
context of the largest Lyapunov exponent (LLE) problem. These include what
Ref. [Zhou and Wang, 2018] calls the “small data method”—used, for example, in
Ref. [Rosenstein, Collins, and De Luca, 1993]—which is efficient but suffers
from the subjectivity problem described above; an object searching method of
the “maximum nonfluctuant," which often reports a local optimal solution
Yongfeng  _et al._ (2012); and a fuzzy C-means clustering approach Shuang  _et
al._ (2016), which also has limitations since the number of clusters must be
pre-specified by the user.

Our approach formalizes the selection of a scaling interval by first choosing
all possible “sensible” combinations of left- and right-hand-side marker
pairs. Specifically, given a sequence of data points at
$\\{x_{j},\;j=1..N\\}$, we choose $x_{l}$ and $x_{r}$ such that

$r-l>n,$ (1)

so that the left marker $x_{l}$ must be below the right one ($x_{r}$) and
there must also be at least $n+1$ data points on which to perform the linear
fit. For each pair $[x_{l},x_{r}]$, we perform a least-squares fit using the
data bounded by the pair to compute both the slope and the least-squares fit
error. To suppress the influence of intervals with poor fits, we construct a
histogram of their slopes, weighting the frequency of each directly by the
length of the fit and inversely by the square of the associated fitting error.
Finally, we use a kernel density estimator to fit a smooth probability density
function (PDF) to the histogram. The details are presented in Appendix A,
along with a discussion of the choices for $n$ and the exponents associated
with the fit width and the fitting error (which are set to 1 and 2,
respectively, in the experiments reported here).

A central claim in this paper is that the resulting distribution contains
useful information about the existence, extent, and slope of the scaling
region. The slopes for which the PDF is large represent estimates that not
only occur for more intervals in the graph, but also that are longer and have
lower fit errors. In particular, we conjecture that the mode of this
distribution gives the best estimate of the slope of the scaling region: i.e.,
a value with low error that persists across multiple choices for the position
of the interval. More formally, most intervals $[x_{l},x_{r}]$ that result in
slopes near the mode of the PDF (within full width half maximum Weik (2001))
correspond to high-quality fits of regions in the graph whose bounds lie
within the interval of the scaling region. If the graph has a single scaling
region, the PDF will be unimodal; a narrower peak indicates that the majority
of high quality fits lie closer to the mode. Broad peaks and multimodal
distributions might arise for multiple reasons, including noise; we address
these issues in §III.

For the example in Fig. 1, there are $N=100$ data points and we set $n=10$ so
the ensemble of endpoints has $4005$ elements. The resulting weighted slope
distribution is unimodal, as shown in panel (b), with mode $4.97$—not too far
from the expected slope of 5.0. To provide a confidence interval around this
estimate, we compute three quantities:

  1. (a)

the full width half maximum ($FWHM$) of the PDF around the mode;

  2. (b)

the fraction, $p_{FWHM}$, of ensemble members that are within the $FWHM$;

  3. (c)

the standard deviation, $\sigma$, for the slopes of the ensemble members
within the $FWHM$.

For the example in Fig. 1, we obtain

$FWHM=0.36,\quad p_{FWHM}=0.67,\quad\mbox{and }\sigma=0.11.$

That is, we estimate the slope of the scaling region as $4.97$ $\pm$ $0.11$,
noting that fully $67\%$ of the estimated values are within $\pm 0.18$. These
error estimates quantify the shape of the peak, and give assurance that the
diagram does indeed contain a scaling region.

The method also can estimate the boundaries of the scaling region. To do this,
we compute secondary distributions using the positions of the left-hand and
right-hand side markers. Since we choose all possible pairs, subject to the
constraint (1), the probability densities of the resulting histograms will
decrease linearly for $x_{l}$ and grow linearly for $x_{r}$ from left to
right. To penalize poor fits, we again scale the frequency for each point in
these secondary histograms directly by the length and inversely by the square
of the error of the corresponding fit. The resulting PDFs for the synthetic
example of Fig. 1(a) are shown in panel (c). The modes of the LHS and RHS
histograms fall at $x_{l}=6.45$ and $x_{r}=8.92$, respectively. These align
neatly with what an observer might choose as the endpoints of the putative
scaling region in Fig. 1(a). Note that the $x_{l}$ distribution is wider with
a relatively lower peak, whereas the $x_{r}$ distribution is narrower and
taller. This indicates that we can be more certain about the position of the
right endpoint than we are about the left one. There is more flexibility in
choosing $x_{l}$ than $x_{r}$.

(a) $y(x)=x^{2}$

(b) Weighted slope histogram with the KD estimator

(c) Weighted distributions of interval endpoints

Figure 2: Estimating the scaling region of a quadratic function.

To test our ensemble method on a dataset without an obvious scaling region,
consider the quadratic function

$y(x)=x^{2}\\\ .$

We choose the range $x\in[0,2]$ and a sampling interval $\Delta x=0.05$. A
graph of resulting function, Fig. 2(a), of course has no straight regions;
though, one might imagine that the curve becomes nearly linear near the top of
the range even though the actual slope grows linearly with $x$. The slope
distribution, shown in Fig. 2(b), has mode

$2.51\pm 0.63\quad(FWHM=2.18\mbox{ and }p_{FWHM}=0.89),$

and as the endpoint distributions in panel (c) show, this corresponds to the
interval $[x_{l},x_{r}]=[0.90,1.70]$. However, the large width of the slope
distribution and the correspondingly high magnitude of $\sigma$ convey little
confidence as to the existence of a scaling region. This shows that the error
statistics and the $FWHM$ can be useful indicators as to whether or not a
scaling region is present. A high error estimate and/or a high FWHM
contraindicates a scaling region.

In the case of the two simple examples above, scaling regions imply simple
linear relationships between the independent and dependent variables. Our
technique can just as well be applied to plots with suitably transformed
axes—like the log-linear or log-log axes involved in the computation of many
dynamical invariants—as we demonstrate next by computing the dimension of a
fractal.

The box-counting dimension, $d_{box}$, is the growth rate of the number
$N(\epsilon)$ of boxes of size $\epsilon$ that cover a set:

$\displaystyle N(\epsilon)\propto\epsilon^{-d_{box}}.$ (2)

In practice, $d_{box}$ is estimated by computing the slope of the graph of
$\ln\ N(\epsilon)$ versus $\ln\tfrac{1}{\epsilon}$.

(a) 100,000 points on the Sierpinski Triangle

(b) Log-log plot of $N(\epsilon)$

(c) Weighted slope distribution

Figure 3: Estimating the box-counting dimension of the Sierpinski Triangle.

As an example, we consider the equilateral Sierpinski triangle of side one,
generating a set of $10^{5}$ points using an iterated function system Peitgen,
Juergens, and Saupe (1992) to obtain Fig. 3(a). A log-log plot of
$N(\epsilon)$ for this data set is shown in Fig. 3(b) and the resulting slope
distribution is shown in Fig. 3(c). The mode of this distribution falls at
$1.585$ $\pm$ $0.020$, where—as before—the error is the estimate $\sigma$.
This is close to the true dimension $\frac{\ln 3}{\ln 2}$ $\approx$ $1.585$.
This peak is narrow, with $FWHM=0.07$, though $p_{FWHM}=0.68$, similar to the
first example. Thus about $70\%$ of the weighted fits lie within $\pm 0.035$
of the estimated slope, strengthening confidence in the estimate. The small
periodic oscillations in the curve in Fig. 3(b) are due to the self-similarity
of the fractal. Note that the linear growth saturates when
$\ln\tfrac{1}{\epsilon}<0.5$, the point at which the box size becomes
comparable to the diameter of the set. Neither of these effects significantly
influences the mode of the slope distribution. For this example, the LHS and
RHS distributions in panel (b) have similar widths, with modes $x_{l}=1.64$
and $x_{r}=5.59$, respectively.

The examples of this section show that the ensemble-based method effectively
selects an appropriate scaling region—if one exists—and is able to exclude
artifacts near the edges of the data.

##  III Applications to Dynamical Systems

In this section, we apply this scaling region identification and
characterization method to calculations of two dynamical invariants—the
correlation dimension (§III.1) and the Lyapunov exponent (§III.2)—for two
well-known examples, the Lorenz system and the driven, damped pendulum. We
also explore the effects of noise (§III.1.3) and the selection of embedding
parameters for delay reconstructions (§III.1.4).

###  III.1 Correlation Dimension

Correlation dimension, one of the most common and useful dynamical invariants,
can be calculated using the Grassberger-Procaccia algorithm Grassberger and
Procaccia (2004, 1983). Given a set of points that sample an object, such as
an attractor of a dynamical system, this algorithm estimates the average
number of points in an $\epsilon$-neighborhood of a given point,
$C(\epsilon)$. For an attractor of correlation dimension $d_{2}$, this scales
with $\epsilon$ as

$C(\epsilon)\propto\epsilon^{d_{2}}.$ (3)

Estimating $d_{2}$ is therefore equivalent to finding a scaling region in the
log-log plot of $C(\epsilon)$ versus $\epsilon$ Grassberger and Procaccia
(1983).

In this section, we use d2, the TISEAN Kantz and Schreiber (1997)
implementation, which outputs estimates of $C(\epsilon)$ for a range of
$\epsilon$ values. We apply our ensemble-based method to a graph of log
$C(\epsilon)$ versus log $\epsilon$ and identify the mode of the resulting
slope distribution to obtain an estimate of the correlation dimension of the
data set. When the slope distribution is multi-modal, our method can also
reveal the existence of potential scaling regions for different ranges of
$\epsilon$, some of which may not be obvious on first observation. This
analysis, then, not only yields values for the slope and extent of the scaling
region; it also provides insight into the geometry of the dynamics, the
details of the Grassberger-Procaccia algorithm, and the interactions between
the two.

####  III.1.1 Lorenz

We start with the canonical Lorenz system:

$\displaystyle\dot{x}$ $\displaystyle=\sigma(y-x),$ $\displaystyle\dot{y}$
$\displaystyle=x(\rho-z)-y,$ $\displaystyle\dot{z}$ $\displaystyle=xy-\beta
z,$

with $\sigma=10$, $\beta=\frac{8}{3}$, and $\rho=28$ Lorenz (1963). We
generate a trajectory from the initial condition $(0,1,1.05)$ using a fourth-
order Runge-Kutta algorithm for $10^{5}$ points with a fixed time step $\Delta
t=0.01$, discarding the first $10^{4}$ points to remove the transient, to get
a total trajectory length of 90,000 points. The chosen initial condition lies
in the basin of the chaotic attractor and the discarded length is sufficient
to remove any transient effects. Figure 4(a) shows a log-log plot of
$C(\epsilon)$ versus $\epsilon$ produced by the d2 algorithm.

(a) Correlation sum calculations

(b) Weighted slope distribution

(c) Weighted distributions of interval endpoints

(d) Heatmap of the endpoint ensemble

Figure 4: Estimating the correlation dimension of the Lorenz attractor

To the eye, this graph appears to contain a scaling region in the approximate
range $[\ln\epsilon_{l},\ln\epsilon_{r}]=[-3,2]$. As in the box-counting
dimension example in §II, the curve flattens when $\epsilon$ is larger than
the diameter of the attractor, since $C(\epsilon)$ will not change when
$\epsilon$ increases beyond this point. Since the diameter of the Lorenz
attractor is $25.5$, this flattening occurs for $\ln\epsilon>\ln(25.5)\approx
3.2$.

Figures 4(b) and (c) show the results of our ensemble-based approach. The mode
of the slope PDF is

$2.09\pm 0.02\quad(FWHM=0.06\mbox{ and }p_{FWHM}=0.69);$

this gives an estimate for the correlation dimension for the Lorenz attractor
that agrees with the accepted range, $2.06-2.16$ Viswanath (2004); Martínez,
Domínguez-Tenreiro, and Roy (1993). The PDFs of the LHS and RHS markers in
Fig. 4(c) have modes $\ln\epsilon_{l}=-2.6$ and $\ln\epsilon_{r}=1.5$,
respectively. Estimates of both endpoints by our algorithm are close to the
informal ones, although slightly on the conversative side. In other words, our
method can indeed effectively formalize the task of identifying both the
existence and the extent of scaling regions.

The endpoint distributions can be broken down further in the heatmap-style
scatter plot of Fig. 4(d) to reveal additional details. Each dot represents a
single element of the ensemble: i.e., a fit for a particular interval
$[\ln\epsilon_{l},\ln\epsilon_{r}]$. Its color encodes the slope and its
radius is scaled by the fitting weight. If samples come from intervals that
are close, the associated dots will be nearby; if these have low error, the
dots will be large and their effective density will be high. Note that the
domain of panel (d) is a triangle since $\epsilon_{l}<\epsilon_{r}$ by (1).

This visualization provides an effective way to identify scaling region ranges
that produce similar slope estimates: long scaling regions will manifest as
large regions of similarly colored points. Fig. 4(d) shows such a triangular
high-density region in blue (slope $\approx$ 2.0), above the diagonal, bounded
from the left by $\ln\epsilon_{l}\approx-4$, and from above by
$\ln\epsilon_{r}\approx 2$. This clearly corresponds to the primary scaling
region in Fig. 4(a) and the corresponding mode in panel (b). This heatmap can
reveal other, shorter scaling regions. Indeed, panel (d) shows other clusters
of similarly colored dots that are somewhat larger than the dots from
neighboring regions: e.g., the green cluster for $\ln\epsilon_{l}\approx 2$
and $\ln\epsilon_{r}\approx 3$, with a slope near $1.5$. Close examination of
the original curve in Fig. 4(a) reveals a small straight segment in the
interval $[2,3]$. This, not immediately apparent, feature stands out quite
clearly in the distribution visualization. Two much smaller scaling regions at
the two ends of panel (a) are also brought out by this representation: the
interval of slope $\approx$ 3.0 for small $\epsilon$ (the dark blue cluster at
the lower left corner of the scatter plot) and the zero slope for large
$\epsilon$ (the red cluster at the upper right corner).

####  III.1.2 Pendulum

As a second example, we study the driven, damped pendulum:

$\begin{split}\dot{\theta}&=\omega,\\\
\dot{\omega}&=-\beta\omega-\nu_{0}^{2}\sin(\theta)+A\cos(\alpha\ t),\\\
\end{split}$ (4)

where the natural frequency is $\nu_{0}=98$, the damping coefficient is
$\beta=2.5$, and the driving force has amplitude $A=91$ and frequency
$\alpha=0.75\nu_{0}$. The coordinates of the three-dimensional extended phase
space $\mathbb{T}^{2}\times\mathbb{R}$ are the angle $\theta$, the time $t$,
and the angular velocity $\omega$. We generate a trajectory for this system
using a fourth-order Runge-Kutta algorithm with initial condition
$(3.14,0,50)$ for $1.1\times 10^{6}$ points with fixed time step $0.001$,
discarding the first $10^{5}$ points to remove the transient, resulting in a
final time series of length one million. To avoid issues with periodicity in
$\theta$ and $t$, we project the time series from
$\mathbb{T}^{2}\times{\mathbb{R}}$ onto $(\sin(\theta),\sin(\alpha
t),\omega)\in\mathbb{{\mathbb{R}}}^{3}$. The resulting trajectory is shown in
Fig. 5. Note that the attractor has the bound $|\omega|\leq 24.88$, but that
the range of the other two variables is $[-1,1]$ due to the sinusoidal
transformation.

Figure 5: A trajectory of the driven, damped pendulum (4).

To the eye, results of a d2 calculation on this trajectory, shown in Fig.
6(a), exhibit two apparent scaling regions above and below $\ln\epsilon\approx
1$. The slope distributions produced by our method not only confirm, but also
formalize, this observation. The larger peak of the distribution in panel (b)
of the figure falls at $d_{2}=2.19\pm 0.09$ ($FWHM$ = 0.32 and
$p_{FWHM}=0.53$), which is equivalent to the correlation dimension of the
attractor as reported in Ref. [Deshmukh  _et al._ , 2020]. Note that, to
account for the fact that the distribution is clipped on the right before the
density falls below half the peak value, the computation of $FWHM$ uses this
as the upper bound. The interval endpoint distributions in panel (c) indicate
that the extent of this scaling region is $-2.6<\ln\epsilon<0$. To the eye,
$\epsilon_{r}\sim 0.8$ would seem a more-appropriate choice for the RHS
endpoint of the scaling region; however, minor oscillations in the interval
$\epsilon\in[0,1]$ cause the ensemble-based algorithm to be more conservative
and choose $\epsilon_{r}=0$. The presence of a second scaling region in panel
(a) for $\ln\epsilon>1$ gives a second mode at the slope $0.94\pm 0.10$
($FWHM=0.35$ and $p_{FWHM}=0.26$). The secondary peaks in the endpoint
distributions in panel (c) suggest an extent of $0.8<\ln\epsilon<2.8$ for this
second scaling region, which is consistent with visual inspection of panel
(a).

This second scaling region is an artifact of the large aspect ratio of the
attractor: only one of the variables, $\omega$, has a range larger than $1$,
so when $\epsilon>1$, the effective dimension is one. To test this hypothesis,
we rescaled the components of the phase space vectors so that the attractor
has equal extent of $[0,1]$ in all three directions (using the -E flag in the
TISEAN d2 command) and repeated the d2 calculations. This leads to rescaling
of the axis of the $d_{2}$ plot—$\ln\epsilon$ now has the range $[-7,0]$—and
causes the second scaling region in the previous results to disappear, leaving
a single scaling region $-7<\ln\epsilon<-2$ with slope

$d_{2}=2.22\pm 0.02\quad(FWHM=0.05,\mbox{ and }p_{FWHM}=0.70).$

Note that, in addition to resolving the artifact of the second scaling region,
this rescaling also reduces the width of the mode and gives a much tighter
error bound.

(a) Correlation sum calculations

(b) Weighted slope distribution

(c) Weighted distributions of interval endpoints with the correlation sum
graph from (a) superimposed

Figure 6: Estimating the correlation dimension of the driven, damped pendulum
trajectory shown in Fig. 5.

By revealing the existence, the endpoints, and the slopes of different scaling
regions in the data, our ensemble-based approach has not only produced an
objective estimate for the correlation dimension, but also yielded some
insight into the interaction between the d2 algorithm and the attractor
geometry.

####  III.1.3 Noise

Noise is a key challenge for any practical nonlinear time-series analysis
method. We explore the effects of measurement noise on our method using the
Lorenz trajectory of §III.1.1 by adding uniform noise post facto to each of
the three state variables. The noise amplitude in each dimension is
proportional to the radius of the attractor in that dimension. The correlation
sum $C(\epsilon)$ for a noise amplitude $\delta=0.01$—i.e., $1\%$ of the
radius of the attractor in each dimension—is shown in Fig. 7(a).

(a) $\delta=0.01$

(b) Weighted slope distribution for (a)

(c) $\delta=0.1$

(d) Weighted slope distribution for (c)

(e) $\delta=0.1$

(f) Weighted slope distribution for (e)

Figure 7: Correlation dimension of the Lorenz attractor with added noise
$\delta=0.01$ for (a)-(b), and $\delta=0.1$ for (c)-(f). Panels (e)-(f)
exclude $\ln\epsilon<-1$. The red markers in (a), (c) and (e) denote scaling
regions corresponding to noise while the green markers indicate that for the
dynamics.

There appear to be two distinct scaling regions in this plot, above and below
a slight knee at $\ln\epsilon\ \approx\ -1$. The slope distribution produced
by our method is shown in Fig. 7(b). As in the pendulum example, this
distribution is bimodal, indicating the presence of two scaling regions with
slopes

$\displaystyle 2.94\pm 0.09$ $\displaystyle\quad(FWHM=0.33\mbox{ and
}p_{FWHM}=0.42),\mbox{ and }$ $\displaystyle 2.08\pm 0.15$
$\displaystyle\quad(FWHM=0.51\mbox{ and }p_{FWHM}=0.30),$

respectively. We claim that these results are consistent with the geometry of
the noise and of the dynamics, respectively. The taller peak corresponds to
the interval $\ln\epsilon\in[-4.60,-1]$, delineated by red markers in panel
(a). Note that the right endpoint of this region is close to
$\epsilon_{r}\approx 0.01\times 26$, the maximum extent of the noise. The
computed slope of $2.94$ in this interval captures the geometry of a noise
cloud in three-dimensional state space. The smaller, secondary peak at $2.08$
reflects the dimension of the attractor for scales larger than the noise, the
interval $\ln\epsilon\in[-1,2.2]$ that is bounded by green markers in panel
(a).

As the noise level increases, the geometry of the attractor is increasingly
obscured (see Fig. 10 in Appendix B). Figures 7(c) and (d) show results for
$\delta=0.1$. As one would expect, the right-hand boundary of the lower
scaling region is increased, shown as the red markers in panel (c). We observe
that $\ln\ \epsilon_{r}\approx 0.16$, near the the noise cloud radius of
$1.3$. With this noise level, the slope distribution is nearly unimodal with
mode

$3.02\pm 0.05\quad(FWHM=0.18\mbox{ and }p_{FWHM}=0.65),$

again reflecting the geometry of the noise. While there appears to be a
secondary peak in the slope distribution in panel (d), it is nearly obscured.

Interestingly, the identification of the scaling region due to the noise can
be used to refine the calculation and retrieve some information about the
dynamical scaling: one simply repeats the slope-distribution calculations
using only the data for larger $\epsilon$ values; that is, discarding the
values below the noise. Restricting to the interval $[-1.0,4.0]$, gives the
curve and the slope distribution shown in panels (e) and (f). In addition to
the noise peak near $3.0$, this reveals two additional peaks. The first,

$2.59\pm 0.20\quad(FWHM=0.68\mbox{ and }p_{FWHM}=0.67),$

was hinted at by the tail of the distribution in panel (d). This, we
conjecture, corresponds to the correlation dimension of the noisy dynamics.
The final peak, near $1.8$, corresponds to the smaller scaling region
$[2.0,3.0]$ that can be seen in a heatmap-style scatter plot (not shown). The
validity of this region can be easily confirmed by further limiting the
ensemble to the interval $\ln\epsilon\in[1.0,4.0]$ which then gives a single
mode at $1.83$. Note that this region corresponds to a smaller scaling region
that was also seen in the noiseless case: the small green cluster in Fig.
4(d).

The ability of the ensemble-based slope distribution method to reveal
secondary scaling regions and suggesting refinements allows one to identify
noise scales and even extract information that might be hidden by the noise.

####  III.1.4 Time Series Reconstructions

The previous examples assumed that all of the state variables are known. This
is rarely the case in practice; rather, one often has data only from a single
sensor. Delay-coordinate embedding Takens (1981); Packard  _et al._ (1980),
the foundation of nonlinear time-series analysis, allows one to reconstruct a
diffeomorphic representation of the actual dynamics from observations of a
scalar time series $x(t)$ provided that a few requirements are met: $x(t)$
must represent a smooth, generic observation function of the dynamical
variablesTakens (1981); Whitney (1936) and the measurements should be evenly
spaced in time. This method embeds the time series $x(t)$ into
$\mathbb{R}^{m}$ as a set of delay vectors of the form
$[x(t),x(t-\tau),\dots,x(t-(m-1)\tau)]$ given a time delay $\tau$.
Theoretically the only requirements are $\tau>0$Takens (1981) and
$m>2d_{cap}$, where $d_{cap}$ is the capacity dimension of the attractor on
which the orbit is dense Sauer, Yorke, and Casdagli (1991). Many heuristics
have been developed for estimating these two free parametersOlbrich and Kantz
(1997); Garland, James, and Bradley (2016); Fraser and Swinney (1986); Kantz
and Schreiber (1997); Buzug and Pfister (1992a); Liebert, Pawelzik, and
Schuster (1991); Gibson  _et al._ (1992); Buzug and Pfister (1992b); Liebert
and Schuster (1989); Garland and Bradley (2015); Rosenstein, Collins, and De
Luca (1994); Cao (1997); Kugiumtzis (1996); Kennel, Brown, and Abarbanel
(1992); Hegger, Kantz, and Schreiber (1999); Deshmukh  _et al._ (2020); Kennel
and Isabelle (1992), notably the first minimum of the average mutual
information for selecting $\tau$ Fraser and Swinney (1986) and the false near
neighbors algorithm for selecting $m$Kennel and Isabelle (1992). See Refs.
[Bradley and Kantz, 2015; Hegger, Kantz, and Schreiber, 2007] for more
details.

Since the correlation dimension is preserved under a diffeomorphism, it can be
computed with the d2 algorithm on properly reconstructed dynamics. Moreover,
the calculations of the correlation dimension for different values of $m$ can
help diagnose the correctness of an embedding. To explore how our method can
contribute in this context, we embed the $x(t)$ time series of the Lorenz
trajectory from §III.1.1 using $\tau=18$, which was chosen using the curvature
heuristic described in Refs. [Deshmukh  _et al._ , 2020] and corroborated
using the first minimum of the average mutual information. Using TISEAN, we
then create a series of embeddings for a range of $m$ values and apply the
Grassberger-Procaccia algorithm to each one. The results are shown in Fig.
8(a) for $m\in[1,10]$.

(a) Correlation sum calculations

(b) Weighted slope distributions

(c) Wasserstein distance

Figure 8: Correlation dimension from time-delay reconstructions of the Lorenz
system for $x(t)$, $\tau=18$ and $m\in[1,10]$. Each $m$ from (a) gives a slope
distribution (b), and these are pairwise compared in (c) using the Wasserstein
distance.

As is well known, we can assess the convergence of $d_{2}$ by reasoning about
such a sequence of curves. When $m$ is too small, the reconstructed attractor
is improperly unfolded, giving an incorrect dimension, while for large enough
$m$, $d_{2}$ typically converges—modulo data limitation issues—to the
nominally correct value. Since our ensemble methodology automates the
calculation of scaling regions, it can simplify this calculation for multiple
curves. To determine the value of $m$ for which the slopes converge, one
simply computes the modes of each slope distribution and looks for the $m$
value at which the positions of those modes stabilize. For Fig. 8(a), the mode
for $m=1$ is, of course, $d_{2}=1$, but when $m$ reaches $3$, the
distributions begin to overlap and for $3\leq m\leq 5$ the modes are $2.09\pm
0.01$, essentially the estimate from the full dynamics in §III.1.1.

To formalize this procedure, we use the Wasserstein metric $M_{W}$ to compare
the two distributionsVallender (1974). Informally, this metric—called the
“earth mover’s distance”—treats each distribution as a pile of soil, with the
distance between the two distributions defined as the minimum “effort”
required to turn one pile into another. Thus $M_{W}=0$ only when the
distributions are identical. We use the python
scipy.stats.wasserstein_distance package to compute this
metric.†††https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html

Given slope distributions $P_{m}$ and $P_{m-1}$ for successive embedding
dimensions from Fig. 8(b), the distance $M_{W}(P_{m},P_{m-1})$ is shown in
panel (c) of the figure. Note that $M_{W}(P_{m},P_{m-1})$ initially decreases
rapidly with increasing $m$, approaching zero at $m=4$. Up to fluctuations,
$M_{W}(P_{m},P_{m-1})$ remains approximately constant thereafter. Through
experiments across multiple examples, we heuristically propose a threshold
$M_{W}(P_{m},P_{m-1})\sim 0.1$ to estimate convergence of $d_{2}$ as $m$
grows. Such a threshold has two benefits, giving both an estimate of the
embedding dimension and of $d_{2}$.

To assess the generality of these ideas, we applied the ensemble-based method
to four additional scalar time series: the noisy Lorenz data of §III.1.3 with
$\delta=0.1$ using $x(t)$, the pendulum trajectory of §III.1.2 using the
$\omega(t)$, a shorter, $200,000$ point segment of this orbit, and finally, a
trajectory from the Lorenz-96 model Lorenz (2006). In each case, we again used
the curvature-based heuristic of Ref. [Deshmukh  _et al._ , 2020] to estimate
$\tau$. This set of examples represents a range of problems that can arise in
practice: measurements affected by noise, and a trajectory that is too short
to fully cover an attractor, and a high-dimensional attractor.

System | Full Dynamics | Reconstructed
---|---|---
| $d_{2}$ |  $d_{2}$ | $\tau$ | $m$
Lorenz (III.1.1) | $2.09\pm 0.02$ | $2.09\pm 0.01$ | 18 | 4
Noisy Lorenz ($\delta=0.1$) | $2.59\pm 0.20$ | $2.27\pm 0.14$ | 21 | 7
Pendulum (III.1.2) | $2.22\pm 0.02$ | $2.16\pm 0.02$ | 120 | 4
Truncated Pendulum | $2.23\pm 0.07$ | $2.20\pm 0.03$ | 120 | 4
Lorenz-96 | $5.63\pm 0.09$ | $5.80\pm 0.07$ | 24 | 10
Table 1: Correlation dimensions calculated using the ensemble-based method
with Wasserstein distance for five cases. The histograms for each of these
examples are shown in Appendix C.

The results, presented in Table 1 show that $d_{2}$ for the reconstructed
dynamics is close to that of the full dynamics as well as that given
previously by manually fitting scaling regions Deshmukh  _et al._ (2020).
However, the embedding dimensions are often significantly smaller that those
suggested by other work. For the first four examples in Table 1, for example,
Ref. [Deshmukh  _et al._ , 2020] used $m=7=2d+1$, the dimension required by
the Takens theorem for a three-dimensional state space Takens (1981); for the
Lorenz-96 example, that paper used $m=12$, obtained from the heuristic
$m>2d_{cap}$ Sauer, Yorke, and Casdagli (1991). Both $m\geq 2d+1$ and
$m>2d_{cap}$ are sufficient conditions, of course. The Wasserstein test used
in Table 1 shows that accurate estimates of $d_{2}$ can often be obtained with
a lower embedding dimension and without prior knowledge of the original
dynamics.

###  III.2 Lyapunov Exponent

Computing a Lyapunov exponent also often involves identification and
characterization of scaling regions. Here we will use the Kantz algorithm
Kantz (1994) to estimate the maximal Lyapunov exponent $\lambda_{1}$ for a
reconstructed trajectories of the Lorenz and chaotic pendulum examples from
§III.1.4. We embed the scalar data using a delay $\tau$ found by standard
heuristics, and experiment with various embedding dimensions. The Kantz
algorithm starts by finding all of the points in a ball of radius
$\epsilon_{s}$ (also called the “scale”) around randomly chosen reference
points on an attractor. By marching through time, the algorithm then computes
the divergence between the forward image of the reference points and the
forward images of the other points in the $\epsilon_{s}$ ball. The average
divergence across all reference points at a given time is the stretching
factor. To estimate $\lambda_{1}$, one identifies a scaling region in the log-
linear plot of the stretching factor as a function of time.

Note that this procedure involves three free parameters: the embedding
parameters $\tau$ and $m$ and the scale of the balls used in the
algorithm.‡‡‡Other parameters—the number of reference points, Theiler window,
and length of the temporal march, etc.—were set to TISEAN’s default values. To
obtain accurate results one is confronted by an onerous multivariate parameter
tuning problem, and an automated approach can be extremely advantageous. If,
for example, one embeds the $x$ coordinate of the Lorenz data from §III.1.1
with $m$ ranging from 2 to 9 and chooses 10 values of $\epsilon_{s}$ (e.g., on
a logarithmic grid with 10 values between 0.038 and 0.381), then there will be
80 different curves, as seen in Fig. 9(a). Manually fitting scaling regions to
each curve would be a demanding task.

(a) Stretching factor calculations (Lorenz)

(b) $\lambda_{1}$ estimate (Lorenz)

(c) Stretching factor calculations (pendulum)

(d) $\lambda_{1}$ estimate (pendulum)

Figure 9: Estimating the largest Lyapunov exponent $\lambda_{1}$ of the Lorenz
and driven, damped pendulum examples. Panels (a) and (c) show stretching-
factor calculations for the eight embedding dimensions $m$ and ten search
scales $\epsilon_{s}$. The estimated exponents for all the dimensions and
scales are shown in grid form in panels (b) and (d). The gray scale in each
grid cell represents the magnitude of the exponent.

The ensemble method gives the automated results shown in panel (b) of the
figure.

Each row of this 2D grid corresponds to fixed $m$ and each column to fixed
$\epsilon_{s}$. The estimated $\lambda_{1}$—the location of the mode in the
slope PDF—is the value shown in the cell. To help detect convergence, each
cell is shaded according to the corresponding $\lambda_{1}$ value. Note that
the majority of the configurations (48 out of 80, colored in intermediate
gray) give an estimate of $\lambda_{1}=0.89\pm 0.02$, which is consistent with
the known value Grigorenko and Grigorenko (2003).

The pattern in the grid make sense. The cells in its center correspond to
midrange combinations of the free parameters: $m\in[5,9]$ and
$\epsilon_{s}\in[0.064,0.177]$. Values outside these ranges create well-known
issues for the Kantz algorithm in the context of finite data. If
$\epsilon_{s}$ is too small, the ball will not contain enough trajectory
points to allow the algorithm to effectively track the divergence; if
$\epsilon_{s}$ is too large, the ball will straddle a large portion of the
attractor, thereby conflating state-space deformation with divergence. For $m$
too small, the attractor is insufficiently unfolded, so the values in the top
few rows of the grid are suspect. In the upper right corner, these two effects
combine: the attractor is both insufficiently unfolded and also spanned by the
$\epsilon_{s}$ ball. Also, the zero slope portion on the right ends of the
stretching factor curves becomes dominant (since since the $\epsilon_{s}$ ball
reaches the edge of a flattened attractor more quickly) thus causing the
ensemble-based algorithm to return a slope close to zero. Finally, we see a
slightly higher estimate of $\lambda_{1}=1.17$ for $m=9$ and
$\epsilon_{s}=0.038$. On close inspection of the stretching factor plot for
this configuration, we see significant oscillations distorting most of the
curve. A relatively straighter part of these oscillations is chosen by the
algorithm as a scaling region, leading to that specific estimate. Seeing how
visually distorted this curve is from the curves generated for other parameter
combinations, it is clear that this parameter combination should be avoided
when computing the $\lambda_{1}$ estimate.

The process is repeated for the driven damped pendulum example in Fig. 9(c,d).
Here, we use a slightly different trajectory than in section III.1.2:
$500,000$ points (after removing a transient of $50,000$ points) with a time
step of $0.01$. For this trajectory, the heuristic of Ref. [Deshmukh  _et al._
, 2020] suggests an embedding delay of $\tau=21$. From the grid, we observe
similar patterns as in the Lorenz example: 46 out of 80 combinations give
$\lambda=0.93\pm 0.04$ (for $m\in[4,9]$ and $\epsilon_{s}\in[0.083,0.3]$). We
observe significantly lower estimates when the embedding dimension is low and
$\epsilon_{s}$ is high. Additionally, an anomalous behavior similar to the
Lorenz example is observed for higher embedding dimensions ($m\in[8,9]$) and
smaller scales ($\epsilon_{s}\in[0.05,0.065]$). Here, we see lower estimates
($0.60\pm 0.02$) for these parameter configurations. The underlying cause in
this case is low frequency oscillations in the stretching factor plots, with a
relatively straight region towards the right end of the plots where they start
to flatten. As for the Lorenz example, these parameter combinations should be
avoided for the $\lambda_{1}$ estimate.

The ensemble approach allows us to easily automate computation of scaling
regions for various values of the hyperparameters, thus sparing significant
manual labor. Moreover, the grid visualization allows us to find the “sweet
spot” in the free-parameter space. While this is only a relatively crude
convergence criterion, one could instead use a more rigorous test such as the
Wasserstein distance of § III.1.4.

##  IV Conclusions

The technique described above is intended to formalize a subjective task that
arises in the calculation of dynamical invariants, the identification and
characterization of a scaling region: an interval in which one quantity grows
linearly with some other quantity. By manipulating the axes, one can extend
this to detect exponential or power-law relationships: e.g., using log-log
plots for fractal dimensions and log-linear plots for Lyapunov exponents, as
shown in Sections III.1 and III.2. Moreover, linearity is not a limitation in
the applicability of our approach: it can be used to identify regions in a
dataset with any hypothesized functional relationship (e.g. higher order
polynomial, logarthmic, exponential, etc.). One would simply compute the least
squares fit over the data using the hypothesized function instead.

A strength of our method is that the scaling region is chosen automatically as
the optimal family of intervals over which the fits have similar slopes with
small errors. To do this, we create an ensemble of fits over different
intervals to build PDFs, to identify the boundaries and slopes of any scaling
regions, and to obtain error estimates. When the scaling region is clear, our
ensemble-based approach computes a slope similar to that which would be
manually determined. We showed that a convergence test on the slope
distributions helps choose the delay reconstruction dimension. This is a
challenging task: even though there are powerful theoretical guidelines, the
information needed to apply them is not known in advance from a given scalar
time-series.

This kind of objective, formalized method is useful not only because it
reveals scaling regions that may not be visible, but also because it provides
confidence in its existence by providing error estimates. Moreover, since
computing a dynamical invariant, such as the Lyapunov exponent, can involve
finding scaling regions in many curves, an automated method provides a clear
advantage over hand selection of scaling regions. Our method could be
potentially useful in areas outside dynamical systems as well: e.g.,
estimating the scaling exponent for data drawn from a power-law distribution.
This is an active area of research Alstott, Bullmore, and Plenz (2014);
Clauset, Shalizi, and Newman (2009). Standard techniques involve inferring
which power-law distribution (if any) the data is most likely drawn from,
and—if so—which exponent is most likely. Our method could potentially help
narrow down the range of exponents to begin such a search.

There are a number of interesting avenues for further work, beginning with how
to further refine the confidence intervals for slope estimates. We have used
the standard deviation of the slopes within the $FWHM$ around the mode of the
distribution. Alternatively, one could use the average of the least squares
fit error across all the samples within the $FWHM$ as the confidence interval.
Another possibility is to first extract a single scaling region (by
determining the modes for the left- and right-hand side markers) and computing
the error statistics over all the sub-intervals that are sampled from this
scaling region (e.g., standard deviation of the slopes or the average of the
fit errors). We used the Wasserstein distance to assess the convergence of a
sequence of slope distributions to help choose an embedding dimension. This
may be too strict since it requires that the _entire_ distributions are close.
One could instead target the positions of the modes, quantifying closeness
using some relevant summary statistic like their confidence interval. Of
course, in cases of multimodal distributions, the Wasserstein distance test is
more appropriate since we cannot choose a single mode for computing the
intervals.

###### Acknowledgements.

The authors acknowledge support from NSF grants EAGER-1807478 (JG),
AGS-2001670 (VD, EB, and JDM) and DMS-1812481 (JDM). JG was also supported by
Omidyar and Applied Complexity Fellowships at the Santa Fe Institute. Helpful
conversations with Aaron Clauset and Holger Kantz are gratefully acknowledged.

## Data Availability

The data that supports the findings of this study are available within the
article.

##  Appendix A Algorithm

The pseudocode for the proposed ensemble-based approach is described in
Algorithm 1. There are four important choices in this algorithm. The first is
the parameter $n$, the minimum spacing between the left-hand and right-hand
side markers $LHS$ and $RHS$. This choice limits what is deemed to be a
“significant interval” for a scaling region. We set $n=10$ in this paper,
which we chose using a persistence analysis, that is varying $n$ and seeing if
the results change. One could argue for increasing $n$ for a data set with
more points; however, pathologies in the data set could violate that argument.
The second is the choice of a kernel density estimator (KDE) for the
histograms. We used the Gaussian KDE with the python implementation
scipy.stats.gaussian_kdeVirtanen  _et al._ (2020) of Scott’s methodScott
(1979) to automatically select the kernel bandwidth. Alternatively, one might
choose other bandwidth selection methods, such as Silverman’s methodSilverman
(1986), or simply manually specify the bandwidth. Scott’s method, which is the
default method used by the package, worked well for all our examples.

Algorithm 1 Ensemble approach for estimating the slope of scaling regions

1:Assume a plot $(x[0:length-1],y[0:length-1])$ that potentially contains a
scaling region.

2:Initialize empty arrays $slopes$, $lengths$, $errors$, $x_{l}$ and $x_{r}$.

3:for $lhs=0,1,2,\ldots,length-1$ do

4: for $rhs=0,1,2,\ldots,length-1$ do

5: if $rhs-lhs>n$ then

6: Fit a line using least squares to data $x[lhs:rhs],y[lhs:rhs]$.

7: Obtain an estimate for slope $m$ and intercept $c$.

8: Compute the least squares fit length and fitting error,
$\displaystyle\text{fit length}=\sqrt{1+m^{2}}\left|x[rhs]-x[lhs]\right|,$ (5)
$\displaystyle\text{fit
error}=\sqrt{\frac{\sum_{i=lhs}^{rhs}(y[i]-mx[i]-c)^{2}}{rhs-lhs}}.$ (6)

9: Append the slope, length and error to the arrays $slopes$, $lengths$ and
$errors$, respectively.

10: Append $x[lhs]$ and $x[rhs]$ to endpoint arrays $x_{l}$ and $x_{r}$,
respectively.

11:Generate a histogram $H$ from the $slopes$ array, weighting each point by
$\text{weight[i]}\propto\frac{\text{(fit length[i])}^{p}}{\text{(fit
error[i])}^{q}},$ (7) for suitable powers $p$ and $q$.

12:Using a kernel density estimator, generate a probability distribution
function $P$, from $H$ (see e.g. Fig. 1(b)).

13:Compute the mode(s) of $P$ and the error estimates as described in Sec. II.

14:Generate histogram $H_{l}$ and $H_{r}$ for $x_{l}$ and $x_{r}$ from Step
10, weighting the frequency with (7). Generate PDFs and evaluate them on the
$x$ interval, as in Step 12, to generate distributions $P_{l}$ and $P_{r}$
(see e.g. Fig. 1(c)).

The two other important parameters are the powers $p$ and $q$ used for the
weights of the length and the error of the fit in (7). We explored ranges for
these parameters, constraining $p$ and $q$ to nonnegative integers, and
determined suitable values that worked well for all of the examples considered
here. Some interesting patterns emerged in these explorations. Firstly, we
found that $p>0$ helps to reduce the error estimate for the slope $\sigma$ and
improve $p_{FWHM}$. Note that $p>0$ penalizes the shorter fits near the edges
of the $FWHM$, suppressing their influence on the histogram. The $FWHM$
therefore narrows, in turn reducing the error estimate $\sigma$. Secondly,
setting $p\leq q$ was found to be advantageous in all cases, ensuring that the
algorithm does not prioritize unnecessarily longer fits of poor quality.
Enforcing these two conditions, we experimented with a few choices for $p$ and
$q$ in the context of the first example of Fig. 1, and found that lower $p,q$
values generally work better. Higher powers tend to magnify effects of small
errors and small lengths, making the algorithm very conservative in terms of
what constitutes a good fit. With this in mind, we settled on $p=1$ and $q=2$
for the correlation dimension estimation, finding that those values
generalized well across all the examples. For the Lyapunov exponent examples,
on the other hand, we found that $p=1$ and $q=1$ does better, generating
tighter confidence interval bounds around the estimate across the various
parameter combinations.

This algorithm always uses the full data set to compute the ensemble but note
that it does not require even spacing of the data. Given a much larger data
set, it might make sense to downsample to speed up the algorithm. In its
current implementation, the run-time complexity of the algorithm is
$O(N^{2})$, for a relationship curve with $N$ points. In the future, it might
be useful to develop faster algorithms for sampling and generating the slope
distributions.

##  Appendix B Effects of noise on the Lorenz attractor

(a) Noiseless Lorenz trajectory

(b) Noisy Lorenz trajectory ($\delta=0.01$)

(c) Noisy Lorenz trajectory ($\delta=0.1$)

Figure 10: Effect of noise on the Lorenz attractor with different magnitudes
of noise level $\delta$.

Fig. 10 shows the effect of noise on the Lorenz attractor, the resulting noisy
attractor used in Section III.1.3. As we increase the noise level (as a
fraction of the attractor radius) $\delta$ from 0.0 to 0.1, we see the
dynamics of the attractor getting progressively distorted. As discussed in
Section III.1.3, the effects of this distortion on the correlation dimension
plots and hence the slope distributions are clearly visible.

##  Appendix C Time Series Reconstruction: Additional Plots

(a) Weighted slope distributions ($\tau=21$)

(b) Wasserstein distance ($\tau=21$)

(c) Weighted slope distributions ($\tau=60$)

(d) Wasserstein distance ($\tau=60$)

Figure 11: Estimating the correlation dimension for the reconstructed noisy
Lorenz example using two different estimates of $\tau$. The top row shows the
slope distributions and the Wasserstein distance profiles for a value of
$\tau=21$, estimated using the heuristic of Ref. [Deshmukh  _et al._ , 2020],
while the bottom row presents the same plots for $\tau=60$, the estimate
produced using the method of average mutual information Fraser and Swinney
(1986).

(a) Weighted slope distributions

(b) Wasserstein distance

Figure 12: Estimating the correlation dimension for the reconstructed pendulum
example (1M points) using a time delay estimate of $\tau=120$.

(a) Weighted slope distributions

(b) Wasserstein distance

Figure 13: Estimating the correlation dimension for the reconstructed
truncated pendulum example (200,000 points) using a time delay estimate of
$\tau=120$.

(a) Weighted slope distributions ($\tau=24$)

(b) Wasserstein distance ($\tau=24$)

Figure 14: Estimating the correlation dimension for the reconstructed
Lorenz-96 example using a time delay estimate of $\tau=24$.

Here we present figures to support the results of Table 1. For the noisy
Lorenz data, Fig. 11 shows the slope distributions and Wasserstein distance
between slope distributions of consecutive embedding dimensions,
$M_{W}(P_{m},P_{m-1})$. Panels (a) and (b) use $\tau=21$ as in Table 1. Note
that $M_{W}(P_{m},P_{m-1})$ converges more slowly for the noisy data than it
did for the deterministic case in Fig. 8. The distance
$M_{W}(P_{m},P_{m-1})<0.1$ at $m=7$, giving $d_{2}=2.27\pm 0.14$.

Figures 11(c) and (d) show the reconstruction results for $\tau=60$, the
embedding delay given by the average mutual information (AMI) method. Here
$M_{W}(P_{m},P_{m-1})$ does not reach $0.1$ for $m\leq 10$. The implication is
that the embedding dimension should be larger than $10$, contrary to the
sufficient theoretical requirement of $m=7$. Indeed, panel (c) shows that the
mode grows monotonically with $m$, reaching values much higher than the
expected $d_{2}=2.59$ for the full dynamics.

Figures 12 and 13 show the results for the the pendulum trajectories of length
$10^{6}$ and $2(10)^{5}$, respectively, generated using a delay of $\tau=120$.
In both cases, the Wasserstein distance threshold is reached at $m=4$,
resulting in the values in Table 1.

Finally, Fig. 14 shows the delay reconstruction of the Lorenz-96 trajectory
using $\tau=24$. Convergence occurs at $m=10$, giving $d_{2}=5.80\pm 0.07$.

## References

  * Kantz and Schreiber (1997) H. Kantz and T. Schreiber,  _Nonlinear Time Series Analysis_ (Cambridge University Press, Cambridge, 1997).
  * Bradley and Kantz (2015) E. Bradley and H. Kantz, Chaos 25, 097610 (2015).
  * Ji, Zhu, and Jiang (2011) C. Ji, H. Zhu, and W. Jiang, Chinese Science Bulletin 56, 925 (2011).
  * Chen, Liu, and Zhou (2019) Z. Chen, Y. Liu, and P. Zhou, Fractals 27, 1950011 (2019).
  * Clauset, Shalizi, and Newman (2009) A. Clauset, C. Shalizi, and M. Newman, SIAM Review 51, 661 (2009).
  * Zhou and Wang (2018) S. Zhou and X. Wang, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 123118 (2018).
  * Rosenstein, Collins, and De Luca (1993) M. Rosenstein, J. Collins, and C. J. De Luca, Physica D: Nonlinear Phenomena 65, 117 (1993).
  * Yongfeng  _et al._ (2012) Y. Yongfeng, W. Minjuan, G. Zhe, W. Yafeng, and R. Xingmin, Vibration, Testing and Diagnosis 32, 371 (2012).
  * Shuang  _et al._ (2016) Z. Shuang, F. Yong, W. Wenyuan, and W. Weihua, Acta Phys 65, 020502 (2016).
  * Weik (2001) M. H. Weik, “full-width at half-maximum,” in  _Computer Science and Communications Dictionary_ (Springer US, Boston, MA, 2001) pp. 661–661.
  * Peitgen, Juergens, and Saupe (1992) H. Peitgen, H. Juergens, and D. Saupe,  _Chaos and Fractals_ (Spinger-Verlag, New York, 1992).
  * Grassberger and Procaccia (2004) P. Grassberger and I. Procaccia, in  _The Theory of Chaotic Attractors_ (Springer, 2004) pp. 170–189.
  * Grassberger and Procaccia (1983) P. Grassberger and I. Procaccia, Physical Review Letters 50, 346 (1983).
  * Grassberger and Procaccia (1983) P. Grassberger and I. Procaccia, Physica D 9, 189 (1983).
  * Lorenz (1963) E. Lorenz, J. Atmos. Sci. 20, 130 (1963).
  * Viswanath (2004) D. Viswanath, Physica D 190, 115 (2004).
  * Martínez, Domínguez-Tenreiro, and Roy (1993) V. Martínez, R. Domínguez-Tenreiro, and L. Roy, Phys. Rev. E 4, 735 (1993).
  * Deshmukh  _et al._ (2020) V. Deshmukh, E. Bradley, J. Garland, and J. D. Meiss, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 063143 (2020).
  * Takens (1981) F. Takens, in  _Dynamical systems and turbulence_ (Springer, Berlin, 1981) pp. 366–381.
  * Packard  _et al._ (1980) N. Packard, J. P. Crutchfield, J. Farmer, and R. Shaw, Phys. Rev. Lett. 45, 712 (1980).
  * Whitney (1936) H. Whitney, Ann. Math. 37, 645 (1936).
  * Sauer, Yorke, and Casdagli (1991) T. Sauer, J. Yorke, and M. Casdagli, J. Statistical Physics 65, 579 (1991).
  * Olbrich and Kantz (1997) E. Olbrich and H. Kantz, Phys. Lett. A 232, 63 (1997).
  * Garland, James, and Bradley (2016) J. Garland, R. G. James, and E. Bradley, Phys. Rev. E 93, 022221 (2016).
  * Fraser and Swinney (1986) A. Fraser and H. Swinney, Phys. Rev. A 33, 1134 (1986).
  * Buzug and Pfister (1992a) T. Buzug and G. Pfister, Physica D 58, 127 (1992a).
  * Liebert, Pawelzik, and Schuster (1991) W. Liebert, K. Pawelzik, and H. Schuster, Europhysics Letters 14, 521 (1991).
  * Gibson  _et al._ (1992) J. Gibson, J. Farmer, M. Casdagli, and S. Eubank, Physica D 57, 1 (1992).
  * Buzug and Pfister (1992b) T. Buzug and G. Pfister, Phys. Rev. A 45, 7073 (1992b).
  * Liebert and Schuster (1989) W. Liebert and H. Schuster, Physics Letters A 142, 107 (1989).
  * Garland and Bradley (2015) J. Garland and E. Bradley, Chaos 25, 123108 (2015).
  * Rosenstein, Collins, and De Luca (1994) M. Rosenstein, J. Collins, and C. De Luca, Physica D 73, 82 (1994).
  * Cao (1997) L. Cao, Physica D 110, 43 (1997).
  * Kugiumtzis (1996) D. Kugiumtzis, Physica D 95, 13 (1996).
  * Kennel, Brown, and Abarbanel (1992) M. Kennel, R. Brown, and H. Abarbanel, Phys. Rev. A 45, 3403 (1992).
  * Hegger, Kantz, and Schreiber (1999) R. Hegger, H. Kantz, and T. Schreiber, Chaos 9, 413 (1999).
  * Kennel and Isabelle (1992) M. Kennel and S. Isabelle, Phys. Rev. A 46, 3111 (1992).
  * Hegger, Kantz, and Schreiber (2007) R. Hegger, H. Kantz, and T. Schreiber, “Tisean 3.0.1 nonlinear time series analysis,”  (2007).
  * Vallender (1974) S. S. Vallender, Theory of Probability & Its Applications 18, 784 (1974).
  * Lorenz (2006) E. Lorenz, in  _Predictability of Weather and Climate_ (Cambridge University Press, 2006) pp. 40–58.
  * Kantz (1994) H. Kantz, Physics Letters A 185, 77 (1994).
  * Grigorenko and Grigorenko (2003) I. Grigorenko and E. Grigorenko, Phys. Rev. Lett. 91, 034101 (2003).
  * Alstott, Bullmore, and Plenz (2014) J. Alstott, E. Bullmore, and D. Plenz, PloS one 9, e85777 (2014).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Scott (1979) D. W. Scott, Biometrika 66, 605 (1979).
  * Silverman (1986) B. W. Silverman,  _Density Estimation for Statistics and Data Analysis_ (Chapman & Hall, London, 1986).


# First-principles prediction of structural, magnetic properties of Cr-
substituted strontium hexaferrite, and its site preference

Binod Regmi, 1,2 Dinesh Thapa, 1,2 Bipin Lamichhane, 1,2 Seong-Gon Kim, 1,2,a)
1 Department of Physics and Astronomy, Mississippi State University,
Mississippi State, Mississippi 39762, USA
2 Center for Computational Sciences, Mississippi State University, Mississippi
State, Mississippi 39762, USA  Author to whom correspondence should be
addressed. Electronic mail<EMAIL_ADDRESS>
###### Abstract

ABSTRACT
To investigate the structural and magnetic properties of Cr-doped M-type
strontium hexaferrite (SrFe12O19) with x = (0.0, 0.5, 1.0), we perform first-
principles total-energy calculations relied on density functional theory.
Based on the calculation of the substitution energy of Cr in strontium
hexaferrite and formation probability analysis, we conclude that the doped Cr
atoms prefer to occupy the 2a, 12k, and 4f2 sites which is in good agreement
with the experimental findings. Due to Cr3+ ion moment, 3 $\mu_{B}$, smaller
than that of Fe3+ ion, 5 $\mu_{B}$, saturation magnetization (Ms) reduce
rapidly as the concentration of Cr increases in strontium hexaferrite. The
magnetic anisotropic field $\left(H_{a}\right)$ rises with an increasing
fraction of Cr despite a significant reduction of magnetization and a slight
increase of magnetocrystalline anisotropy $\left(K_{1}\right)$.The cause for
the rise in magnetic anisotropy field $\left(H_{a}\right)$ with an increasing
fraction of Cr is further emphasized by our formation probability study. Cr3+
ions prefer to occupy the 2a sites at lower temperatures, but as the
temperature rises, it is more likely that they will occupy the 12k site. Cr3+
ions are more likely to occupy the 12k site than the 2a site at a specific
annealing temperature (>700°C).

##  I  Introduction

Hexaferrites, also known as hexagonal ferrites or hexagonal ferrimagnets, are
a class of magnetic materials that have been of great interest to researchers
since their discovery in the $1950$’s. These hexaferrites are found in
numerous types such as M, Y, Z, W, X, U-type commonly doped with zinc,
strontium, nickel, aluminum, and magnesium. The most common properties to all
hexaferrites include that all are ferrimagnetic, their properties of magnetism
are based on the crystal structure, and they take different amounts of energy
to magnetize in a specific direction within the crystal because of spin-orbit
interaction [1].
Particularly, we are interested in M-type strontium hexaferrite (SrFe12O19,
SFO) that falls to space group $P63/mmc$ which has a crystal structure of
hexagonal magnetoplumbite. The unit cell of SFO having two formula units is
presented in Fig(1).

Figure 1: (a) One double formula unit cell of
$\mathrm{SrFe}_{12}\mathrm{O}_{19}$. Small maroon spheres are $\mathrm{O}$
atoms, and two huge gray spheres are Sr atoms. $\mathrm{Fe}^{3+}$ ions are
represented by colored spheres encircled by polyhedra made of $\mathrm{O}$
atoms in a variety of inequivalent sites: $2a$ (red), $2b$ (blue), $4f_{1}$
(green), $4f_{2}$ (yellow), and $12k$ (purple).(b) A schematic representation
of the $\mathrm{Fe}^{3+}$ ions of $\mathrm{SrFe}_{12}\mathrm{O}_{19}$ in their
lowest energy spin configuration. The local magnetic moment at each atomic
location is represented by the arrows.

The iron ions in this structure are coupled in a tetrahedral, trigonal
bipyramidal, and octahedral manner by oxygen ions.The magnetic property is
retained in SFO mainly due to the occupancy of Fe+3 ions in five inequivalent
sites (namely 2a, 2b, 4f1, 4f2, and 12k), three octahedral sites (namely 2a,
12k, and 4f2), one trigonal bipyramidal site (2b), and one tetrahedral site
(4f1). However, the range (degree) of magnetic properties will be influenced
by the shape and size of material particles mainly in the context of thin
films, nanoparticles [2, 3, 4, 5].
In SFO, there is an involvement of interactions between the moments or moments
with lattice ions, tend to contribute anisotropic energy which is termed as
magnetocrystalline anisotropy (MA). In a more explicit way, it is the
dependence of the magnetic properties on the applied field direction relative
to the crystal lattice. Magnetocrystalline anisotropy energy (MAE), an
integral property of a ferromagnetic crystal, is the energy difference to
magnetize a crystal along easy and hard direction of magnetization [2, 6]. The
primary source of MA is spin-orbit coupling (SOC). Because of this coupling,
orbitals of electrons are coupled with the spin of electrons and follow the
spin direction no matter how the magnetization changes its direction in space
[6]. The anisotropy that arises on a crystal are mainly due to the shape of a
magnetic particle at the quantum scale, atomic diffusion at sufficiently high
temperature, and the interaction between a ferromagnetic and an
antiferromagnetic materials [6, 7, 8, 9].
SFO is one of the best candidates among hexaferrite groups owing to its
industrial and electronic implications. In the beginning years, Technophiles
were motivated about SFO to make permanent magnets, recording media, electric
motors because of its high saturation magnetization, large coercivity, optimum
Curie temperature, supreme magnetocrystalline anisotropy, and more chemical
stability. Nowadays, due to technological advancement, it is equally growing
interest in the development of nano fibres, electronic components for mobile
and wireless communications. Recently, researchers have been characterized Y,
M, U, and Z ferrites as multiferroics even at room temperature. These
multiferroics have a wide range of practical implications like multi-state
memory elements, memory media, and novel functional sensors [2, 10, 11].
Several successful investigations have been done on SFO to understand its
electronic structure by computational and experimental approach. To uplift the
strength of magnetic and electric properties, researchers are substituted ions
or pair of ions in a different concentration mainly on Fe sites of SFO.
Majority of researchers have followed the non-magnetic ions substitution into
Fe sites to enhance further the value of saturation magnetization (Ms). In
case of Zr-Cd substitution (SrFe12-2x(ZrCd)xO19), the value of Ms augmented in
the limit of concentration $x=0.2$, whereas the value of coercivity declined
with increasing concentration of Zr-Cd [12]. Substitution of Er-Ni pair in SFO
showed the continuous rise in the value of Ms and coercivity in accordance
with the concentration [13]. However, the substitution of certain pairs like
Zn-Nb,[14] Zn-Sn, [15, 16, 17] and Sn-Mg [18, 19] showed the increasing trend
of Ms and decreasing pattern of coercivity.
In this study, we performed the first-principles total-energy calculations to
analyze the link between site occupation and magnetic properties of
substituted strontium hexaferrite, SrFe12-xCrxO19 with $x=0.5$ and $x=1.0$.
Every configuration of substituted SFO appears with a particular probability.
To determine the formation probabilities of its various configurations at a
typical annealing temperature ($1000K$), we used the Boltzmann distribution
function. We show that our calculation predicts a decrease of saturation
magnetization ($M_{s}$) as well as a decrease in magnetic anisotropy energy
(MAE) of SrFe12-x(Cr)xO19 at $x=0.5$ and $1.0$ compared to the pure M-type
SFO. This result is also in good agreement with the experimental observation
as observed in Ghasemi et al. (2009) [19].

##  II Computational details

We pursued the first-principles total-energy calculations for configuration
SrFe12-x(Cr)xO19 at $x=0.5$ and $1.0$. In our calculation, a unit cell of two
formula units of SFO is used. The structural optimization calculation along
with total energies and forces were carried out using density functional
theory by projector augmented wave (PAW) potential as executed in VASP.
Depending on the ground state ferrimagnetic spin ordering of Fe, our
calculations were solely focused on spin-polarized [3, 20]. The expansion of
the wave function was in the form of plane waves with a $520eV$ energy cut-off
used for pristine SFO, Cr-substituted SFO. A $7\times 7\times 1$ Monkhorst-
Pack k-mesh was used to sample a Brillouin zone with a Fermi-level smearing of
$0.2eV$ applied through the Methfessel-Paxton method [21, 22]. The electronic
relaxation was accomplished till the change in free energy and the band
structure energy less than $10^{-7}eV$. In addition, we fully optimized the
structure by relaxing the positions of ions and cell shape till the change in
total energies between two ionic steps less than $10^{-4}eV$. We used the
Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) to
describe fully the electron exchange-correlation effect [23]. Furthermore, we
implemented the $GGA+U$ method in the simplified rotationally invariant
approach described by Dudarev et al. to address the localized $3d$ electrons
in Fe [24]. The necessity of Ueff for Fe was fulfilled by setting $3.7eV$
based on previous study. we set Ueff to be zero for all other elements [25].
To evaluate the magnetocrystalline anisotropy energy, we first carried out an
accurate collinear calculation in the ground state, and then we followed the
spin-orbit coupling calculations in two different spin orientations within
non-collinear setup. The substitution of foreign atoms in five
crystallographic inequivalent $\mathrm{Fe}$ sites can change the magnetic
characteristics of SFO. When foreign atoms are substituted in a SFO unit cell,
there are a variety of energetically distinct configurations. The magnetism of
substituted SFO is highly dependent on the on-site preferences of the replaced
atoms since SFO is ferrimagnetic. Understanding the site preference of
substituted atoms is crucial in order to research how substitution affects
magnetic characteristics. The substitution energy can be calculated to find
the substituted atom’s preferred site. The substitution energy
$E_{\text{sub}}[i]$ for configuration $i$ at $0\mathrm{~{}K}$ is given by

$E_{\mathrm{sub}}[i]=E_{\mathrm{CSFO}}[i]-E_{\mathrm{SFO}}-\sum_{\beta}n_{\beta}\epsilon_{\beta}$
(1)

where $E_{\mathrm{CSFO}}[i]$ is the total energy per unit cell of Cr-
substituted SFO in configuration $i$, whereas $E_{\mathrm{SFO}}$ is the total
energy per unit cell of pure $\mathrm{SFO}$ and $\epsilon_{\beta}$ is the
total energy per atom for element $\beta$ ($\beta$= $\mathrm{Cr}$ and
$\mathrm{Fe})$ in its most stable crystal structure. $n_{\beta}$ is the number
of atoms of type $\beta$ added or removed; if one atom is added,
$n_{\beta}=+1$, and when one atom is withdrawn, $n_{\beta}=-1$. The
calculation of Magnetic Anisotropy Energy (MAE) is important for understanding
the preferred magnetization directions in a material. Mathematically, MAE is
defined as the difference between the two total energies where the spin
quantization axes are oriented along two distinct directions: [26]

$E_{a}=E_{(100)}-E_{(001)}$ (2)

where $E_{(100)}$ is the total energy with the spin quantization axis in the
magnetically hard axis and $E_{(001)}$ is the total energy with the spin
quantization axis in the magnetically easy axis. The total energies in Eq.(2)
are computed by the non-self-consistent calculations, where the spin densities
are kept constant. With the help of MAE, the uniaxial magnetic anisotropy
constant, $K_{1}$, can be computed as [27, 28]

$K_{1}=\frac{E_{a}}{V\sin^{2}\theta}$ (3)

where $V$ is the equilibrium volume of the unit cell and $\theta$ is the angle
between the two spin quantization axis orientations $\left(90^{\circ}\right.$
in the present scenario). The anisotropy field, $H_{a}$, which is related to
the coercivity can be expressed as [29]

$H_{a}=\frac{2K_{1}}{M_{s}}$ (4)

where $K_{1}$ is the magnetocrystalline anisotropy constant and $M_{s}$ is the
saturation magnetization.

When the difference in substitution energies $E_{\text{sub }}$ between
different configurations is relatively small compared to the thermal energy at
high annealing temperatures ($\gtrsim 1000\mathrm{~{}K})$, the site preference
of substituted atoms in hexaferrite can change. This change in site occupation
preference can be described using the Maxwell-Boltzmann distribution, which
relates to the formation probability. The site occupation probability or the
formation probability $P_{i}(T)$ of configuration $i$ at temperature $T$ is
given by

$\displaystyle P_{i}(T)$ $\displaystyle=\frac{g_{i}\exp\left(-\Delta
G_{i}/k_{B}T\right)}{\sum_{j}g_{j}\exp\left(-\Delta G_{j}/k_{B}T\right)},$ (5)
$\Delta G_{i}=\Delta E_{i}+P\Delta V_{i}-T\Delta S_{i},\\\ $ (6) $\Delta
S_{i}=k_{B}\ln\left(g_{i}\right)-k_{B}\ln\left(g_{0}\right),$ (7)

where $\Delta G_{i},\Delta E_{i},\Delta V_{i}$, and $\Delta S_{i}$ are the
change in free energy, substitution energy, unit cell volume, and entropy of
the configuration $i$ relative to the ground state configuration. $P$,
$k_{\mathrm{B}}$, and $g_{i}$ are the pressure, Boltzmann constant, and
multiplicity of configuration i. $g_{0}$ is the multiplicity of the ground
state configuration. we considered $\Delta S_{i}$ to be the same for all
configurations based on prior literature [2]. Eq.(7) enhances the model
through the explicit computation of the entropy change concerning the most
stable configuration [30, 31].

Hence, when the probability of higher energy configurations becomes
significant at the annealing temperature, it can be inferred that in a
substituted SFO sample, multiple configurations exist rather than a single
one. Consequently, any physical quantity of the SFO sample will be a weighted
average of the corresponding properties in these different configurations.

$\langle Q\rangle=\sum_{i}P_{1000\mathrm{~{}K}}(i)\cdot Q_{i}$ (8)

where $P_{1000\mathrm{~{}K}}(i)$ and $Q_{i}$ are the formation probability at
$1000\mathrm{~{}K}$ and a physical quantity $Q$ of the configuration $i$. The
weighted average calculated by Eq.(8) is the material’s low-temperature
property even though $1000\mathrm{~{}K}$ is used for computation because the
crystalline configurations of CSFO will be distributed according to this value
during the annealing process.

##  III  Results and Discussion

In order to visualize the doping effect of the $Cr^{+3}$ ion on the structural
and magnetic properties of SFO, we replaced the $Fe^{+3}$ from various lattice
locations. We found that it has a significant effect on the structural and
magnetic properties of this system. We fully relaxed the volume, ionic
positions, and shape of SFO and Cr-doped SFO. The crystal structure remains
hexagonal under all circumstances. We further carried out our calculations
when we match the optimized lattice parameter of pure SFO (a = 5.928 Å, c =
23.195 Å) with experimental lattice constants (a = 5.890 Å, c = 23.182 Å);
less than 1$\%$ difference between the calculated and experimental values. For
x = 1.0, the calculated lattice parameters (a = 5.930 Å, c = 23.076 Å) were
found to be very consistent with experimental lattice parameters (a = 5.902 Å,
c = 23.024 Å) [32, 33]. Fig.2 shows the variation of lattice parameters from
theory and experiment explicitly.

Figure 2: Comparision of predicted and experimental lattice constant of unit
cell as a function of fraction of Cr (x).

To compensate for the unavailability of the experimental lattice parameter at
x = 0.5, we are comparing the experimental values at x = 0.6 with our
calculated values at x = 0.5 for the comparison in Fig. 2. The substitution of
Cr in SFO has little effect on the lattice parameters or unit cell volume.
Because the radius of Cr+3(0.630 Å) is similar to that of Fe+3(0.645 Å), this
is an expected outcome.

In this paper, we measured the various physical quantities by varying the
concentration of Cr in the unit cell of SFO. For the x = 0.5, one Cr atom was
substituted at one of the 24 Fe sites of the unit cell. Many of these Fe sites
are equivalent when crystallographic symmetry operations are applied, leaving
only five inequivalent structures. We label these inequivalent configurations
[2a], [2b], [$4f_{1}$], [$4f_{2}$], and [12k] using the crystallographic name
of the Fe site. These structures were found by fully optimizing the unit cell
shape, volume, and ionic positions. We estimated the substitution energy,
$E_{sub}$, to comprehend the site preference of the substituted Cr atom.

Table 1: Physical properties of inequivalent configurations of $SrFe_{12-x}$CrxO19 with $x=0$ and $0.5$: Doped amount(x), multiplicity($g$), substitution energy($E_{sub}$), total magnetic moment ($M_{tot}$), volume of the unit cell ($V$), saturation magnetization ($M_{s}$), magnetocrystalline anisotropy energy ($E_{a}$), uniaxial magnetic anisotropy constant ($K_{1}$), anisotropy field ($H_{a}$), and the formation probability at $1000$ K ($P_{1000K}$). All values are for a double formula unit cell containing $64$ atoms. x | config. | g | $E_{sub}(eV)$ | $m_{tot}{(\mu_{B})}$ | Volume(Å3) | $M_{s}(emu/g)$ |  $E_{a}$(meV) | $K_{1}(KJ/m^{3})$ | $H_{a}(koe)$ | $P_{1000K}$
---|---|---|---|---|---|---|---|---|---|---
0.0 | SFO | - | - | 39.99 | 706.18 | 105.19 | 0.87 | 197.12 | 7.50 | 1.000
0.5 | $[2a]$ | 2 | -3.41 | 38.00 | 704.70 | 100.00 | 1.20 | 272.45 | 10.89 | 0.076
| $[12k]$ | 12 | -3.33 | 38.00 | 705.04 | 100.00 | 0.90 | 204.25 | 8.17 | 0.884
| $[4f_{2}]$ | 4 | -3.27 | 42.00 | 703.70 | 111.00 | 0.90 | 204.63 | 7.39 | 0.049
| $[2b]$ | 2 | -2.96 | 38.00 | 703.17 | 100.00 | 0.60 | 136.52 | 5.44 | 0.000
| $[4f_{1}]$ | 4 | -2.33 | 42.00 | 705.49 | 111.00 | 0.20 | 45.36 | 1.64 | 0.000

Table I displays the results of our calculation for each of the five
inequivalent configurations in ascending order of substitution
energy($E_{sub}$). The configuration [2a] has the lowest $E_{sub}$ followed by
[12k], and [$4f_{2}$] which is consistent with the experimental outcomes[32,
34]. We can conclude that the [2a] site is the most preferred site for the Cr
atom at 0 K. We used Eq.(5) to calculate the probability of forming each
configuration as a function of temperature. Because the change in volume
between different configurations is so small (less than 0.3 Å3), we may
discard the $P\Delta V$ term as negligible (in the order of
$10^{-7}\mathrm{eV}$ at a standard pressure of 1 atm) compared to the $\Delta
E_{\text{sub}}$ (i) term in Eq.(6). The entropy change $\Delta S$ has two
components: configurational, $\Delta S_{c}$, and vibrational, $\Delta
S_{\mathrm{vib}\cdot}{}$ [31]. $\Delta S_{\mathrm{vib}}$ is around 0.1-0.2
$k_{\mathrm{B}}$ /atom for binary substitutional alloys like the present
system, and $\Delta S_{c}$ is $0.1732k_{\mathrm{B}}$ /atom. As a result, we
assign $\Delta S=0.3732k_{\mathrm{B}}/$ atom.

Figure 3: Temperature dependence of the formation probability of different
configurations of $SrFe_{12-x}$CrxO19 with $x=0.5$. In the plot, the
configurations with trivial probabilities are not displayed.

Fig.(3) shows the formation probability of different configurations of
$\mathrm{SrFe}_{12-x}\mathrm{Cr}_{x}\mathrm{O}_{19}$ with $x=0.5$ at different
temperatures. The foreign $\mathrm{Cr}^{3+}$ ions prefer to occupy host
$\mathrm{Fe}^{3+}$ ions from the $2a$ and the $12k$ sites. The formation
probabilities of $[2b],\left[4f_{1}\right]$, and $\left[4f_{2}\right]$ are
trivial and not displayed in Fig.(3). A $\mathrm{Cr}^{3+}$ ion has a 100$\%$
probability of occupying the [$2a$] site at 0K and it’s value sharply declines
as the temperature rises, while the probability of occupancy of
$\mathrm{Cr}^{3+}$ at the [$12k$] site is maximum (88 $\%$) at 1000K. At a
typical annealing temperature of $1000\mathrm{~{}K}$ for $\mathrm{CSFO}$, the
site occupation probability of the site $[2a]$ and $[12k]$ are $7\%$ and
$88.4\%$ respectively. In CSFO, the doped $\mathrm{Cr}^{3+}$ ions are more
likely to occupy $\mathrm{Fe}^{3+}$ ions at the [$12k$] site than the [$2a$]
site because of the higher multiplicity of the [12k] site.

Table 2: Physical properties of inequivalent configurations of $SrFe_{12-x}$(Cr)xO19 with $x=1.0$: Doped amount(x), substitution energy($E_{sub}$), total magnetic moment ($m_{tot}$),volume of the unit cell ($V$), saturation magnetization ($M_{s}$), magnetocrystalline anisotropy energy ($E_{a}$), uniaxial magnetic anisotropy constant ($K_{1}$), anisotropy field ($H_{a}$), and the formation probability at $1000$ K ($P_{1000K}$). All values are for a double formula unit cell containing $64$ atoms. x | config. | $E_{sub}(eV)$ | $m_{tot}{(\mu_{B})}$ | Volume(Å3) | $M_{s}(emu/g)$ |  $E_{a}$(meV) | $K_{1}(KJ/m^{3})$ | $H_{a}(koe)$ | $P_{1000K}$
---|---|---|---|---|---|---|---|---|---
1.0 | $[2a,2a]$ | -6.84 | 36.00 | 703.73 | 95.00 | 1.50 | 341.39 | 14.36 | 0.001
| $[12k,2a]$ | -6.75 | 30.00 | 703.03 | 79.20 | 1.10 | 250.34 | 12.64 | 0.178
| $[12k,12k]$ | -6.69 | 36.00 | 703.24 | 95.00 | 0.80 | 182.01 | 7.66 | 0.664
| $[2a,4f_{2}]$ | -6.68 | 40.00 | 701.63 | 106.00 | 1.30 | 296.45 | 11.20 | 0.009
| $[12k,4f_{2}]$ | -6.62 | 34.00 | 701.89 | 89.70 | 0.90 | 205.16 | 9.13 | 0.145
| $[4f_{2},4f_{2}]$ | -6.56 | 44.00 | 700.33 | 116.00 | 1.00 | 228.46 | 7.84 | 0.001
| $[12k,2b]$ | -6.30 | 30.00 | 701.93 | 79.20 | 0.50 | 113.97 | 5.75 | 0.001
| $[2b,4f_{2}]$ | -6.15 | 34.00 | 698.58 | 89.70 | 0.70 | 160.33 | 7.10 | 0.000
| $[2a,2b]$ | -6.09 | 30.00 | 701.90 | 79.20 | 0.40 | 91.18 | 4.60 | 0.000
| $[2b,2b]$ | -5.99 | 36.00 | 699.37 | 95.00 | 2.00 | 457.55 | 19.15 | 0.000
| $[12k,4f_{1}]$ | -5.66 | 40.00 | 704.40 | 106.00 | 0.20 | 45.43 | 1.72 | 0.000
| $[4f_{1},4f_{2}]$ | -5.61 | 44.00 | 703.98 | 116.00 | 0.20 | 45.46 | 1.57 | 0.000
| $[2a,4f_{1}]$ | -5.47 | 40.00 | 704.37 | 106.00 | 0.90 | 204.44 | 7.76 | 0.000
| $[2b,4f_{1}]$ | -5.15 | 38.91 | 702.14 | 103.00 | 0.10 | 22.79 | 0.88 | 0.000
| $[4f_{1},4f_{1}]$ | -4.67 | 44.00 | 705.71 | 116.00 | 0.60 | 136.03 | 4.70 | 0.000
Table 3: Weighted averages of physical properties of pure and $\mathrm{Cr}$-doped strontium hexaferrite: volume of the unit cell $(\mathrm{V})$, total magnetic moment $\left(m_{\text{tot }}\right)$, saturation magnetization $\left(M_{s}\right)$, magnetocrystalline anisotropy energy ($\left.E_{a}\right)$), uniaxial magnetic anisotropy constant $\left(K_{1}\right)$, and anisotropy field $\left(H_{a}\right)$. The column labeled $\left(M_{s}\right)$ displays the calculated values at 0K along with the bracketed experimental values from K. Praveena et al. [35] for the comparison. Material |  $V$ (Å3) |  $m_{\text{tot}}$ ($\mu_{\mathrm{B}}$) |  $M_{s}$ (emu/g) |  $E_{a}$ (meV) |  $K_{1}$ (kJ/m3) |  $H_{a}$ (kOe)
---|---|---|---|---|---|---
$\mathrm{SrFe}_{12}\mathrm{O}_{19}$ | $706.18$ | $39.99$ | $105.19(59.33)$ | $0.87$ | $197.12$ | $7.50$
$\mathrm{SrFe}_{11.5}\mathrm{Cr}_{0.5}\mathrm{O}_{19}$ | $704.39$ | $38.16$ | $100.46(36.01)$ | $0.91$ | $208.58$ | $8.30$
$\mathrm{SrFe}_{11.0}\mathrm{Cr}_{1.0}\mathrm{O}_{19}$ | $702.07$ | $34.63$ | $91.40(30.04)$ | $0.87$ | $198.44$ | $8.79$

For the x = 1.0, two Cr atoms were substituted at two of the 24 Fe sites of
the unit cell. Many of these Fe sites are equivalent when crystallographic
symmetry operations are applied, leaving only 15 inequivalent structures.
These structures were found by fully optimizing the unit cell shape, volume,
and ionic positions. To understand the site preference of a substituted Cr
atom, we estimated the substitution energy, $E_{sub}$.

Figure 4: Temperature dependence of the formation probability of different
configurations of $SrFe_{12-x}$CrxO19 with $x=1.0$. In the plot, the
configurations with trivial probabilities are not displayed.

Table II provides the results of our calculation for each of the fifteen
inequivalent configurations in ascending order of substitution
energy($E_{sub}$). The configuration [2a, 2a] has the lowest $E_{sub}$
followed by [12k, 2a], and [12k, 12k]. We also used Eq.(5) to calculate the
probability of forming each configuration as a function of temperature.
Because the change in volume between different configurations is so small
(less than 0.7Å3), we may discard the $P\Delta V$ term as negligible (in the
order of $10^{-7}\mathrm{eV}$ at a standard pressure of 1 atm) compared to the
$\Delta E_{\text{sub }}$ (i) term in Eq.(6).

Fig. (4) shows the formation probability of different configurations of
$\mathrm{SrFe}_{12-x}\mathrm{Cr}_{x}\mathrm{O}_{19}$ with $x=1.0$ at different
temperatures. The foreign $\mathrm{Cr}^{3+}$ ions prefer to occupy host
$\mathrm{Fe}^{3+}$ ions from the $[2a,2a]$, $[12k,2a]$,and $[12k,12k]$ sites.
The formation probabilities of other sites are trivial and not displayed in
Fig.(4). A $\mathrm{Cr}^{3+}$ ion has a 100$\%$ probability of occupying the
$[2a,2a]$ site at 0K and it’s value sharply declines as the temperature rises,
while the probability of occupancy of $\mathrm{Cr}^{3+}$ at the $[12k,12k]$
site is maximum (66.4 $\%$) at 1000 K. At a typical annealing temperature of
$1000\mathrm{~{}K}$ for $\mathrm{CSFO}$, the site occupation probability of
the site $[12k,2a]$ is $17.8\%$. In CSFO, the doped $\mathrm{Cr}^{3+}$ ions
are more likely to occupy $\mathrm{Fe}^{3+}$ ions at the $[12k,12k]$ site than
the $[2a,2a]$ site because of it’s higher multiplicity. We exclusively utilize
the formation probability at elevated temperatures for computing weighted
averages. This is because the arrangement of CSFO configurations during the
annealing process will be distributed according to these values. Table III
displays the weighted average of corresponding quantities as the concentration
of $\mathrm{Cr}^{3+}$ increases. The volume of CSFO decreases as we increase
the concentration of $\mathrm{Cr}^{3+}$ because of the smaller atomic radius
of $\mathrm{Cr}^{3+}$ ion. The magnetic moment of CSFO is also in the
decreasing trend as the amount of doped $\mathrm{Cr}^{3+}$ increases owing to
its smaller magnetic moment (3$\mu_{B}$) than the $\mathrm{Fe}^{3+}$
(5$\mu_{B}$). Similarly, the saturation magnetization is decreasing
monotonically as we increase the $\mathrm{Cr}^{3+}$ concentration which is
consistent with the K. Praveena et al.[35]. Although the value of
magnetocrystalline anisotropy $\left(K_{1}\right)$ is slightly increased, the
reduction in saturation magnetization $\left(M_{s}\right)$ is much more
significant. So, their resultant effect causes the anisotropy field
$\left(H_{a}\right)$ to increase as the fraction of Cr is raised. In Table IV,
we have provided the atomic contribution from each sublattice to the overall
magnetic moment of CSFO. It can be observed that the total magnetic moment of
the unit cell is slightly distinct from the sum of local magnetic moments.
This disparity arises from the contribution of the interstitial region to the
overall magnetic moment.

Table 4: Contribution of atoms in each sublattice to the total magnetic moment of $\mathrm{Cr}$-substituted strontium hexaferrite structures [ $12\mathrm{k}]$ and $[12k,12k]$ compared with pure SFO. All magnetic moments are in $\mu_{B}$, and $\Delta m$ is measured relative to the values for pure SFO. | SFO |  | x = 0.5, [12 k] |  | x = 1.0, [12 k, 12 k]
---|---|---|---|---|---
Site | Atoms | m |  | Atoms | m | $\Delta{m}$ |  | Atoms | m | $\Delta{m}$
2d | 2 Sr | -0.006 |  | 2 Sr | -0.005 | 0.001 |  | 2 Sr | -0.005 | 0.001
2a | 2 Fe | 8.328 |  | 2 Fe | 8.334 | 0.006 |  | 2 Fe | 8.336 | 0.008
2b | 2 Fe | 8.130 |  | 2 Fe | 8.131 | 0.001 |  | 2 Fe | 8.129 | -0.001
$4f_{1}$ | 4 Fe | -16.184 |  | 4 Fe | -16.159 | 0.025 |  | 4 Fe | -16.141 | 0.043
$4f_{2}$ | 4 Fe | -16.416 |  | 4 Fe | -16.410 | 0.006 |  | 4 Fe | -16.415 | 0.001
12k | 1 Fe | 4.180 |  | 1 Cr | 2.680 | 0.319 |  | 1 Cr | 2.689 | 0.311
| 1 Fe | 4.178 |  | 1 Fe | 4.180 | 0.000 |  | 1 Cr | 2.692 | 0.308
| 10 Fe | 41.794 |  | 10 Fe | 41.792 | -0.002 |  | 10 Fe | 41.805 | 0.011
4e | 4 O | 1.400 |  | 4 O | 1.292 | -0.108 |  | 4 O | 1.165 | -0.235
4f | 4 O | 0.356 |  | 4 O | 0.294 | -0.062 |  | 4 O | 0.239 | -0.117
6h | 6 O | 0.160 |  | 6 O | 0.151 | -0.009 |  | 6 O | 0.142 | -0.018
12k | 24 O | 3.140 |  | 24 O | 2.777 | -0.363 |  | 24 O | 2.422 | -0.718
$\Sigma{m}$ |  | 39.053 |  |  | 37.056 |  |  |  | 35.057 |
$m_{tot}$ |  | 40 |  |  | 38 |  |  |  | 36 |

##  IV  Conclusions

First-principles total-energy calculations based on density functional theory
were used to study Cr-substituted SFO
($\mathrm{SrFe}_{12-x}\mathrm{Cr}_{x}\mathrm{O}_{19}$) with $x=0.0,0.5,1.0$.
The results showed that increasing the fraction of Cr atoms reduced the total
magnetic moment of the SFO unit cell. This reduction in magnetization was
obtained by low magnetic Cr atoms replacing $\mathrm{Fe}^{3+}$ ions at two of
the majority spin sites, 2a and 12k, resulting negative contribution to the
magnetization. Our substitution energy and formation probability analysis
predicts that Cr atoms preferentially occupy the 2a, 12k, and 4f2 sites,
consistent with experimental observations. Increasing the fraction of Cr in
SFO leads to a rise in the magnetic anisotropic field $\left(H_{a}\right)$
despite a decrease in magnetization and a slight increase in
magnetocrystalline anisotropy $\left(K_{1}\right)$. This increase in
anisotropic field$\left(H_{a}\right)$ is supported by a formation probability
study, which shows that at higher temperatures (>700°C), $\mathrm{Cr}^{3+}$
ions are more likely to occupy the 12k site rather than the 2a site because of
its higher multiplicity.

###### Acknowledgements.

This work was supported by the Center for Computational Science (CCS) at
Mississippi State University. Computer time allocation has been provided by
the High-Performance Computing Collaboratory $(HPC^{2})$ at Mississippi State
University.

## References

  * Pullar [2012] R. C. Pullar, “Hexagonal ferrites: a review of the synthesis, properties and applications of hexaferrite ceramics,” Progress in Materials Science 57, 1191–1334 (2012).
  * Dixit  _et al._ [2017] V. Dixit, S.-G. Kim, J. Park, and Y.-K. Hong, “Effect of ionic substitutions on the magnetic properties of strontium hexaferrite: A first principles study,” AIP Advances 7, 115209 (2017).
  * Fang  _et al._ [2003] C. Fang, F. Kools, R. Metselaar, R. De Groot, _et al._ , “Magnetic and electronic properties of strontium hexaferrite srfe12o19 from first-principles calculations,” Journal of Physics: Condensed Matter 15, 6229 (2003).
  * Soria  _et al._ [2019] G. D. Soria, P. Jenus, J. Marco, A. Mandziak, M. Sanchez-Arenillas, F. Moutinho, J. E. Prieto, P. Prieto, J. Cerdá, C. Tejera-Centeno, _et al._ , “Strontium hexaferrite platelets: a comprehensive soft x-ray absorption and mössbauer spectroscopy study,” Scientific reports 9, 1–13 (2019).
  * Obradors  _et al._ [1988] X. Obradors, X. Solans, A. Collomb, D. Samaras, J. Rodriguez, M. Pernet, and M. Font-Altaba, “Crystal structure of strontium hexaferrite srfe12o19,” Journal of Solid State Chemistry 72, 218–224 (1988).
  * Aharoni [1996] A. Aharoni,  _Introduction to the Theory of Ferromagnetism_ (Clarendon Press, New York, 1996).
  * Miyazaki and Jin [2012] T. Miyazaki and H. Jin,  _The physics of ferromagnetism_ , Vol. 158 (Springer Science & Business Media, 2012).
  * Meiklejohn and Bean [1956] W. H. Meiklejohn and C. P. Bean, “New magnetic anisotropy,” Physical review 102, 1413 (1956).
  * Hong  _et al._ [2019] Y. Hong, Z. Qiu, Z. Zheng, G. Wang, H. Yu, D. Chen, G. Han, D. Zeng, and J. Liu, “Influencing mechanisms of atomic diffusion and compositional distribution on the magnetic anisotropy of cr/smco/(cu)/cr thin films,” Acta Materialia 164, 627–635 (2019).
  * Li  _et al._ [2016] Q. Li, J. Song, M. Saura-Múzquiz, F. Besenbacher, M. Christensen, and M. Dong, “Magnetic properties of strontium hexaferrite nanostructures measured with magnetic force microscopy,” Scientific reports 6, 1–7 (2016).
  * Tan and Chen [2013] G. Tan and X. Chen, “Structure and multiferroic properties of barium hexaferrite ceramics,” Journal of Magnetism and Magnetic Materials 327, 87–90 (2013).
  * Ashiq, Iqbal, and Gul [2009] M. N. Ashiq, M. J. Iqbal, and I. H. Gul, “Structural, magnetic and dielectric properties of zr–cd substituted strontium hexaferrite (srfe12o19) nanoparticles,” Journal of Alloys and Compounds 487, 341–345 (2009).
  * Ashiq  _et al._ [2012] M. N. Ashiq, M. J. Iqbal, M. Najam-ul Haq, P. H. Gomez, and A. M. Qureshi, “Synthesis, magnetic and dielectric properties of er–ni doped sr-hexaferrite nanomaterials for applications in high density recording media and microwave devices,” Journal of magnetism and magnetic materials 324, 15–19 (2012).
  * Fang  _et al._ [2004] Q. Fang, H. Bao, D. Fang, and J. Wang, “Temperature dependence of magnetic properties of zinc and niobium doped strontium hexaferrite nanoparticles,” Journal of applied physics 95, 6360–6363 (2004).
  * Ghasemi  _et al._ [2010] A. Ghasemi, V. Šepelák, X. Liu, and A. Morisako, “The role of cations distribution on magnetic and reflection loss properties of ferrimagnetic srfe 12- x (sn 0.5 zn 0.5) x o 19,” Journal of Applied Physics 107, 09A734 (2010).
  * Ghasemi and Šepelák [2011] A. Ghasemi and V. Šepelák, “Correlation between site preference and magnetic properties of substituted strontium ferrite thin films,” Journal of magnetism and magnetic materials 323, 1727–1733 (2011).
  * Dixit  _et al._ [2019] V. Dixit, D. Thapa, B. Lamichhane, C. N. Nandadasa, Y.-K. Hong, and S.-G. Kim, “Site preference and magnetic properties of zn-sn-substituted strontium hexaferrite,” Journal of Applied Physics (2019).
  * Davoodi and Hashemi [2011] A. Davoodi and B. Hashemi, “Magnetic properties of sn–mg substituted strontium hexaferrite nanoparticles synthesized via coprecipitation method,” Journal of Alloys and Compounds 509, 5893–5896 (2011).
  * Ghasemi  _et al._ [2009] A. Ghasemi, V. Sepelak, X. Liu, and A. Morisako, “Microwave absorption properties of mn–co–sn doped barium ferrite nanoparticles,” IEEE Transactions on Magnetics 45, 2456–2459 (2009).
  * Gorter [1957] E. Gorter, “Saturation magnetization of some ferrimagnetic oxides with hexagonal crystal structures,” Proceedings of the IEE-Part B: Radio and Electronic Engineering 104, 255–260 (1957).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, “Special points for brillouin-zone integrations,” Physical review B 13, 5188 (1976).
  * Methfessel and Paxton [1989] M. Methfessel and A. Paxton, “High-precision sampling for brillouin-zone integration in metals,” Physical Review B 40, 3616 (1989).
  * Perdew, Burke, and Ernzerhof [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Physical review letters 77, 3865 (1996).
  * Dudarev  _et al._ [1998] S. Dudarev, G. Botton, S. Savrasov, C. Humphreys, and A. Sutton, “Electron-energy-loss spectra and the structural stability of nickel oxide: An lsda+ u study,” Physical Review B 57, 1505 (1998).
  * Liyanage  _et al._ [2013] L. S. Liyanage, S. Kim, Y.-K. Hong, J.-H. Park, S. C. Erwin, and S.-G. Kim, “Theory of magnetic enhancement in strontium hexaferrite through zn–sn pair substitution,” Journal of magnetism and magnetic materials 348, 75–81 (2013).
  * Ravindran  _et al._ [1999] P. Ravindran, A. Delin, P. James, B. Johansson, J. M. Wills, R. Ahuja, and O. Eriksson, “Magnetic, optical, and magneto-optical properties of mnx (x=as, sb, or bi) from full-potential calculations,” Phys. Rev. B 59, 15680–15693 (1999).
  * Muñoz  _et al._ [2013] F. Muñoz, A. H. Romero, J. Mejía-López, and J. L. Morán-López, “Finite size effects on the magnetocrystalline anisotropy energy in fe magnetic nanowires from first principles,” Journal of Nanoparticle Research 15 (2013).
  * Smit and Wijn [1954] J. Smit and H. Wijn, “Physical properties of ferrites,” in  _Advances in Electronics and Electron Physics_ , Vol. 6 (Elsevier, 1954) pp. 69–136.
  * Kittel [1949] C. Kittel, “Physical theory of ferromagnetic domains,” Rev. Mod. Phys. 21, 541–583 (1949).
  * Gilmore [2014] C. Gilmore,  _Materials science and engineering properties_ (Cengage Learning, 2014).
  * van de Walle and Ceder [2002] A. van de Walle and G. Ceder, “The effect of lattice vibrations on substitutional alloy thermodynamics,” Rev. Mod. Phys. 74, 11–45 (2002).
  * Fang  _et al._ [2005] Q. Fang, H. Cheng, K. Huang, J. Wang, R. Li, and Y. Jiao, “Doping effect on crystal structure and magnetic properties of chromium-substituted strontium hexaferrite nanoparticles,” Journal of Magnetism and Magnetic Materials 294, 281–286 (2005).
  * Kimura  _et al._ [1990] K. Kimura, M. Ohgaki, K. Tanaka, H. Morikawa, and F. Marumo, “Study of the bipyramidal site in magnetoplumbite-like compounds, srm12o19 (m = al, fe, ga),” Journal of Solid State Chemistry 87, 186–194 (1990).
  * Parker [1980] R. Parker, “Ferrite,” Proceedings of ICF-3 375 (1980).
  * Praveena  _et al._ [2015] K. Praveena, M. Bououdina, M. P. Reddy, S. Srinath, R. Sandhya, and S. Katlakunta, “Structural, magnetic, and electrical properties of microwave-sintered cr3+-doped sr hexaferrites,” Journal of Electronic Materials 44, 524–531 (2015).
  * Nethala  _et al._ [2018] G. P. Nethala, R. Tadi, G. R. Gajula, K. N. C. Kumar, and V. K. Veeraiah, “Investigations on the structural, magnetic and mossbauer properties of cerium doped strontium ferrite,” Physica B: Condensed Matter (2018).