# Temperature scaling analysis of the 3D disordered Ising model with power-law
correlated defects

S. Kazmin  $|$ ${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$,
W. Janke  $|$ ${}^{\textsf{\footnotesize{\mbox{?}}}}$
<EMAIL_ADDRESS>author<EMAIL_ADDRESS>leipzig.de.

(Received November 12, 2022, in final form January 17, 2023)

###### Abstract

We consider the three-dimensional site-diluted Ising model with power-law
correlated defects and study the critical behaviour of the second-moment
correlation length and the magnetic susceptibility in the high-temperature
phase. By comparing, for various defect correlation strengths, the extracted
critical exponents $\nu$ and $\gamma$ with the results of our previous finite-
size scaling study, we consolidate the exponent estimates.

Key words: 3D site-diluted Ising model, long-range correlations, Monte Carlo
simulation, temperature scaling, critical exponents

###### Abstract

Ìè ðîçãëÿäàìî òðèâèìðíó ðîçâåäåíó ìîäåëü çíãà ç ñòåïåíåâèìè êîðåëÿöÿìè äåôåêòâ
òà äîñëäæóìî êðèòèíó ïîâåäíêó äðóãîãî ìîìåíòó êîðåëÿöéíî¿ äîâæèíè ìàãíòíî¿
ñïðèéíÿòëèâîñò ó âèñîêîòåìïåðàòóðíé ôàç. Ïîðâíþþè îòðèìàí äëÿ ðçíèõ
íòåíñèâíîñòåé êîðåëÿö¿ äåôåêòâ êðèòèí ïîêàçíèêè $\nu$ òà $\gamma$ ç
ðåçóëüòàòàìè íàøîãî ïîïåðåäíüîãî äîñëäæåííÿ ñêíåííî-âèìðíîãî ñêåéëíãó, ðîáèìî
óçãîäæåí îöíêè öèõ ïîêàçíèêâ.

Ключов слова: òðèâèìðíà ðîçâåäåíà ìîäåëü çíãà, äàëåêîñÿæí êîðåëÿö¿,
ìîäåëþâàííÿ Ìîíòå-Êàðëî, òåìïåðàòóðíèé ñêåéëíã, êðèòèí ïîêàçíèêè

##  1 Introduction

It is well known that under certain conditions quenched disorder can affect
the critical behaviour of a physical system. Most extensively studied is
uncorrelated disorder for which the Harris criterion [1] states that
impurities are irrelevant when the specific-heat exponent $\alpha_{\rm pure}$
of the pure system is negative, whereas for $\alpha_{\rm pure}>0$
renormalization-group arguments suggest a modified critical behaviour. This
prediction has been confirmed in numerous studies of different models and
especially also for the three-dimensional Ising model [2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12] for which $\alpha_{\rm pure}\approx 0.1102$ [13, 14].

In realistic physical systems, however, it is more likely that the impurities
or defects exhibit some kind of spatial correlations. When these correlations
decay sufficiently slowly, e.g., they follow asymptotically, for large
distances $r$, the power law $r^{-a}$ with a correlation exponent $a<d$, where
$d$ is the dimension of the system, one observes a new scenario for long-range
correlated (quenched) disorder: An extension of the Harris criterion by
Weinrib and Halperin [15] and the later considerations [16, 17, 18] predict
that in this case the correlation-length exponent $\nu$ obeys quite generally

$\nu=2/a\,.$ (1.1)

Similar predictions for other critical exponents read
$\alpha=2-d\nu=2(a-d)/a$, $\beta=(2-\epsilon)/a+{\cal O}(\epsilon^{2})$, and
$\gamma=4/a+{\cal O}(\epsilon^{2})$, where $\epsilon=4-d$ together with
$\delta=4-a$ enters the employed $\epsilon$-$\delta$ renormalization-group
expansion [15]. As already speculated in [15], the relation (1.1) plays a
special role and is expected to be valid to all orders of this expansion [16].
More recently this relation was confirmed for the special case of the two-
dimensional Ising model via a mapping to Dirac fermions and applying an
alternative renormalization-group scheme with a double expansion in
$\epsilon^{\prime}=2-d$ and $\delta^{\prime}=2-a$ up to two-loop order,
$\nu=2/a+O(\delta^{\prime 3})$ [19]. On the other hand, when the correlations
decay more rapidly, e.g., exponentially or power-law-like with $a\geq d$, one
falls back into the universality class of uncorrelated disorder.

In [20, 21], we studied the power-law correlated case for the site-diluted
three-dimensional (3D) Ising model with extensive Monte Carlo (MC) computer
simulations in the vicinity of criticality by employing finite-size scaling
(FSS) techniques for the data analyses. Here, we complement these studies by
reporting alternative estimates for the critical exponents $\nu$ and $\gamma$
obtained from the analyses of the temperature scaling of MC data for the
second-moment correlation length and magnetic susceptibility when approaching
the critical point in the high-temperature phase.

The rest of the paper is organized as follows. In section 2 we briefly recall
the employed model and the simulation method. Our results are presented and
discussed in section 3, and in section 4 we conlude the paper with a summary
and brief outlook to future work.

##  2 Model and methods

The three-dimensional Ising model with site disorder is defined by the
Hamiltonian

${\cal H}=-J\sum_{\langle ij\rangle}\epsilon_{i}\epsilon_{j}s_{i}s_{j},$ (2.1)

where the spins $s_{i}$ take on the values $\pm 1$ and the sum runs over all
nearest-neighbor pairs denoted by $\langle ij\rangle$ of a simple-cubic
lattice of size $V=L^{3}$ with periodic boundary conditions. The defect
variables are $\eta_{i}=1-\epsilon_{i}=0$ when a spin is present at site $i$
and $\eta_{i}=1$ when site $i$ is empty, i.e., occupied by a defect. The
coupling constant is set to $J=1$, fixing the unit of energy and, by setting
the Boltzmann constant $k_{\text{B}}=1$, also the temperature scale.

For uncorrelated disorder, the defects are chosen randomly according to the
probability density

$f(\eta)=p\delta_{\eta,0}+p_{d}\delta_{\eta,1},$ (2.2)

where $\delta_{i,j}$ is the Kronecker delta symbol. Here, $p_{d}$ denotes the
concentration of defects and $p=1-p_{d}$ is the concentration of spins.111Note
that in the corresponding definition in [20, 21], $p$ and $p_{d}$ are
inadvertently interchanged. We use the grand-canonical approach where the
desired defect concentration $p_{d}$ is the mean value over all the considered
disorder realizations.

For correlated disorder, we additionally introduce a long-range spatial
correlation between the defects at sites $i$ and $j$ that decays
asymptotically for large distances $r_{ij}$ according to the power law,

$\langle\eta_{i}\eta_{j}\rangle_{c}\equiv\langle(\eta_{i}-\langle\eta_{i}\rangle)(\eta_{j}-\langle\eta_{j}\rangle)\rangle\propto\frac{1}{r_{ij}^{a}}\,,$
(2.3)

where $a>0$ is the correlation exponent. For the numerical generation of the
defect correlation we employed the Fourier filter method described by Makse et
al. [22, 23] in the publicly available C++ implementation222The C++ code is
available at github.com/CQT-Leipzig/correlated_disorder. of [24], which for
technical reasons considers a slightly modified correlation function
$C(r)=(1+r^{2})^{-a/2}$ that agrees asymptotically with (2.3) (see also [25]).
The resulting $1/r^{2}$ corrections in combination with finite-size effects
make the measurement of the actual correlation exponent $\overline{a}\approx
a$ an important analysis step, cf. table 1; for details we refer to [20].

We considered the correlation exponents $a=1.5$, $2.0$, $2.5$, $3.0$, $3.5$,
and $\infty$ (standing symbolically for the uncorrelated case) and studied in
each case the eight defect concentrations $p_{d}=0.05$, $0.1$, $0.15$, $0.2$,
$0.25$, $0.3$, $0.35$, and $0.4$. For each disorder realization, the MC
simulations of this model were performed at various temperatures $T$ with the
Swendsen-Wang multiple-cluster update algorithm [26], collecting $N=10\,000$
measurements after 500 thermalization sweeps. All final results are the
averages over $N_{c}=1000$ randomly chosen disorder realizations. The linear
lattice size was taken for all temperatures $T$ to be $L=256$, the largest
lattice of our FSS studies [20, 21]. Since the correlation length $\xi(T)$ and
hence finite-size effects quickly diminish away from the critical point, one
could in principle adapt the lattice size to satisfy $L\gg\xi(T)$. However, in
the case of correlated defects, the measured correlation exponent
$\overline{a}$ was found to be slightly $L$-dependent [20] so that mixing
different lattice sizes in scaling analyses of the high-temperature data
should be avoided.

We studied two observables, the second-moment correlation length $\xi$
calculated as [27]

$\xi=\frac{1}{2\sin(\piup/L)}\sqrt{\frac{S(\mathbf{0})}{S(\mathbf{1})}-1},$
(2.4)

where $S(\mathbf{k})$ is the discrete Fourier transform of the spatial spin-
spin correlation function $\langle s_{i}s_{j}\rangle$ in the high-temperature
phase evaluated at $\mathbf{0}=(0,0,0)$ and $\mathbf{1}=(2\piup/L,0,0)$, and
the (high-temperature) susceptibility

$\tilde{\chi}=\beta V\langle m^{2}\rangle,$ (2.5)

where $m=(1/V)\sum_{i}\epsilon_{i}s_{i}$ is the magnetization density and
$\beta\equiv 1/T$. Note that $S(\mathbf{0})=\tilde{\chi}/\beta$.

When $T$ approaches the critical temperature $T_{c}$, the expected
temperature-scaling behaviour of the disorder-averaged observables (indicated
by $[\dots]$) reads

$\displaystyle[\xi(T)]$ $\displaystyle=$ $\displaystyle
a|t|^{-\nu}(1+\dots)\,,$ (2.6) $\displaystyle[\tilde{\chi}(T)]$
$\displaystyle=$ $\displaystyle b|t|^{-\gamma}(1+\dots),$ (2.7)

where $t=(1-T/T_{c})$ $(\leqslant 0)$ is the reduced temperature and
$(1+\dots)$ indicates analytical and confluent scaling corrections which
vanish as $T\to T_{c}$.

##  3 Results

Using the sufficiently precise estimates of $T_{c}$ from [21], we performed
linear fits of $\ln[\xi(T)]$ and $\ln[\tilde{\chi}(T)]$ in $\ln|t|$ which
provided us with the estimates for critical exponents $\nu$ and $\gamma$,
respectively. In what follows, we describe the analysis steps for the
observable $\xi$ and the exponent $\nu$ in some detail and then present an
analogous brief discussion for $\tilde{\chi}$ and the exponent $\gamma$ .

(a) $a=\infty$ (b) $a=3.5$ (c) $a=3.0$ (d) $a=2.5$ (e) $a=2.0$ (f) $a=1.5$ (g)

Figure 1: (Colour online) Correlation length $[\xi(T)]$ as function of the
reduced temperature $t$ for all considered correlation exponents $a$ and
defect concentrations $p_{d}$. The definition (2.4) of $\xi$ is valid only in
the high-temperature phase with $t\leqslant 0$, but we extended the curves in
order to see the crossing points at $t\approx 0$ better.

###  3.1 Critical exponent $\nu$

By plotting the disorder averaged correlation length $[\xi(T)]$ for various
defect concentrations $p_{d}$ as a function of $t$, we visually verified the
used $T_{c}$ estimates (that depend on both $a$ and $p_{d}$), as can be seen
in figure 1 where a negative reduced temperature $t\leqslant 0$ corresponds to
the high-temperature phase. All the curves for different $p_{d}$ intersect at
$t\approx 0$. Only for the strongest correlations with $a\leqslant 2.0$,
visually not all of them intersect in one point. Strictly speaking, $\xi$ is
defined only in the high-temperature phase where $t\leqslant 0$, but we
extended it to $t>0$ in order to see the intersections better.

Considering only the high-temperature values with $t\leqslant 0$, we performed
for each correlation exponent $a$ and each defect concentration $p_{d}$ an
individual fit with the ansatz

$\displaystyle\ln[\xi(t)]=A-\nu\ln|t|\,.$ (3.1)

Since the power-law behaviour only starts at a certain distance away from
$t=0$, i.e., once finite-size effects become neglectable, we varied the
smallest $|t|_{\rm min}$ included in the fits from its smallest value near
$t=0$ to a maximum value where only three degrees of freedom remained.
Examples of the fits are presented in figure 2 which shows one main problem
with this procedure. We clearly observe finite-size effects for each $p_{d}$,
since the data points curve down as $|t|\to 0$. Compare, e.g., the plot for
$a=1.5$ where this effect is most pronounced. However, since the statistical
errors (estimated with the Jackknife method [28]) are quite large, linear fits
still provide reasonable $\chi^{2}_{\rm red}$ values per degree of freedom.
The resulting estimates of the exponent $\nu$ together with the $\chi^{2}_{\rm
red}$ values for each $a$ and $p_{d}$ and for all the considered $|t|_{\rm
min}$ are shown in figure 3. The data show a clear dependence on the defect
concentration $p_{d}$ and also on $|t|_{\rm min}$. They do not reach the
plateau value even for the largest $|t|_{\rm min}$ and the estimates for
$p_{d}\leqslant 0.1$ are clearly influenced by the crossover to the pure Ising
model. This reflects the observation in figure 1 that $[\xi]$ exhibits the
strongest curvature for the smallest defect concentration $p_{d}$. Therefore,
for a quantitative comparison, we computed the error weighted mean
$\overline{\nu}_{\rm ts}^{w}$ over all estimates for $p_{d}\geqslant 0.15$
where for each $p_{d}$ we used the fit with the largest possible $|t|_{\rm
min}$ having three degrees of freedom. This may be not an optimal solution but
at least it was closer to possible plateau values than using the fits with the
smallest possible $|t|_{\rm min}$ for which $\chi^{2}_{\rm red}\leqslant 1.0$
was satisfied for the first time. The latter results are way too low and
clearly do not represent the asymptotic behaviour. As mentioned above, this is
due to the relatively large statistical errors which made the simple linear
fits acceptable, even though finite-size effects were still present.

The weighted means $\overline{\nu}_{\rm ts}^{w}$ are compared in the narrow
right-hand panels of figure 3 with the estimates from our FSS analysis:
Weighted means $\overline{\nu}^{w}_{\rm lin}$ of individual linear fits
neglecting the scaling corrections [29] and $\nu^{\rm g}$ from non-linear
“global” fits including corrections-to-scaling [29, 21, 30]. The numerical
values are compiled in table 1 where we additionally include the weighted
means $\overline{\nu}^{w}$ of non-linear FSS fits including the scaling
corrections [29]. Except for $a=1.5$, the estimates
$\overline{\nu}^{w}_{\text{ts}}$ obtained from temperature scaling are closer
to $\overline{\nu}^{w}_{\rm lin}$ than to $\overline{\nu}^{w}$ respectively
$\nu^{\rm g}$ and are slightly larger. The value for $a=1.5$ is possibly
smaller because the estimates for larger $p_{d}$ show very large errors and
hence the smallest included defect concentration $p_{d}=0.15$ dominates the
weighted mean. Hence, the estimate for this value of $a$ should be considered
with some reservation even though the exemplary fits for $a=1.5$ displayed in
figure 2(f) do not look particularly worrying. The biggest deviations can be
seen for the two correlation exponents $a=3.0$ and $3.5$ which is exactly the
same behaviour as for the two types of FSS estimates, i.e.,
$\overline{\nu}^{w}_{\rm lin}$ and $\overline{\nu}^{w}$ or $\nu^{\rm g}$. We
interpret this as a signal for the theoretically expected crossover from the
correlated to effectively uncorrelated behaviour at $a\approx d=3$.

(a) $a=\infty$ (b) $a=3.5$ (c) $a=3.0$ (d) $a=2.5$ (e) $a=2.0$ (f) $a=1.5$ (g)

Figure 2: (Colour online) Examples of the fits of $[\xi(t)]$ with the ansatz
(3.1). For each concentration of defects $p_{d}$ the smallest $|t|_{\rm min}$
for which $\chi^{2}\leqslant 1.0$ was true for the first time is used in the
plots.

Although one should take the estimates $\overline{\nu}^{w}_{\text{ts}}$ with
some care, we nevertheless can qualitatively confirm the FSS results in all
the considered cases. The prediction (1.1) of Weinrib and Halperin [15] that
$\nu=2/a$ is not matched quantitatively. The results lie above this
prediction, but the dependence on $a$ respectively the measured $\overline{a}$
is clearly in accordance with the FSS results which indeed show a $\propto
1/\overline{a}$ behaviour [21]. For a comparison with previous results for
selected cases of $a$ and $p_{d}$ by other groups [31, 32, 33, 34, 35], we
refer to table I and to the discussion in [20, 21]. Let us finally note that
we also have checked the influence of the statistical error of the $T_{c}$
estimates on the results, but it turned out that it can be neglected due to
much larger errors coming from the fits themselves.

###  3.2 Critical exponent $\gamma$

For the fits of the susceptibility $[\tilde{\chi}(T)]$ in the high-temperature
phase, we used the ansatz

$\displaystyle\ln[\tilde{\chi}(t)]=B-\gamma\ln|t|$ (3.2)

and again performed individual fits for all correlation exponents $a$ and
defect concentrations $p_{d}$. As in the case of $\xi(T)$, we varied the
minimal $|t|_{\rm min}$ included in the fits to see the asymptotic behaviour.

(a) $a=\infty$ (b) $a=3.5$ (c) $a=3.0$ (d) $a=2.5$ (e) $a=2.0$ (f) $a=1.5$ (g)

Figure 3: (Colour online) Final fit results of $[\xi(t)]$ using the ansatz
(3.1) for different $|t|_{\rm min}$ and all concentrations of defects $p_{d}$.
The weighted means $\overline{\nu}^{w}_{\text{ts}}$ over all concentrations
with $p_{d}\geqslant 0.15$ and for the largest possible $|t|_{\rm min}$ are
shown together with the results $\overline{\nu}^{w}_{\rm lin}$ and $\nu^{\rm
g}$ from the FSS analysis [21, 29]. The narrow right-hand panels show a
separate comparison between the different estimates for $\nu$ which are
plotted as horizontal lines in the main plots.

The susceptibility as function of $t$ is presented in figure 4. It has the
same qualitative behaviour as $\xi$—for different $p_{d}$ the curves cross
each other at $t\approx 0$. Note that as in the case of $\xi$, the definition
of $\tilde{\chi}$ is valid only in the high-temperature phase, and we only
extended the $T$ values below $T_{c}$ in order to see the crossing points
better. In figure 5 we show examples of the fits for all correlation exponents
$a$ and defect concentrations $p_{d}$. The estimates of the critical exponent
$\gamma$ in dependence on the chosen $|t|_{\rm min}$ are presented in figure
6. The error weighted means of $\gamma$ over all $p_{d}\geqslant 0.15$ with
the largest $|t|_{\rm min}$ are summarized in table 1.

(a) $a=\infty$ (b) $a=3.5$ (c) $a=3.0$ (d) $a=2.5$ (e) $a=2.0$ (f) $a=1.5$ (g)

Figure 4: (Colour online) Susceptibility $[\tilde{\chi}(T)]$ as function of
the reduced temperature $t$. The definition of $\tilde{\chi}$ is valid only in
the high-temperature phase with $t\leqslant 0$, but we extended the curves in
order to see the crossing points better.

(a) $a=\infty$ (b) $a=3.5$ (c) $a=3.0$ (d) $a=2.5$ (e) $a=2.0$ (f) $a=1.5$ (g)

Figure 5: (Colour online) Examples of the fits of $[\tilde{\chi}(t)]$ with the
ansatz (3.2). For each concentration of defects $p_{d}$ the minimum $|t|_{\rm
min}$ for which $\chi^{2}\leq 1.0$ was true for the first time is used in the
plots.

The first observation is the same as in the case of $\xi$: The smallest
concentrations $p_{d}\leqslant 0.1$ show a crossover behaviour and therefore
we excluded them in the weighted mean. Again, the curves do not reach the
asymptotic values even for the largest $|t|_{\rm min}$. The final weighted
mean estimates $\overline{\gamma}^{w}_{\text{ts}}$ lie slightly above the
corrected global fit estimates $\gamma^{\rm g}$ from the FSS analysis except
for the case of $a=1.5$. They match very well for the correlation exponents in
the range $2.0\leqslant a\leqslant 3.0$ but do not agree well in the
uncorrelated case, $a=\infty$. The crossover region with $a\approx 3.0-3.5$
shows the largest deviations between $\overline{\gamma}^{w}_{\text{ts}}$ and
$\gamma^{\rm g}$. In general, however, the qualitative cross-check with
temperature scaling does support our estimates from the FSS analysis [21].
Unfortunately, here we cannot compare with the individual uncorrected fit
ansatz in the FSS case, since we have not performed it for the critical
exponent $\gamma$.

Table 1: Comparison of final estimates for the critical exponents $\nu$ and $\gamma$ from temperature-scaling (labeled with “ts”) and FSS analyses for several (input) correlation exponents $a$ together with their actual numerically measured values $\overline{a}$. The weighted means $\overline{\nu}^{w}_{\text{ts}}$ and $\overline{\gamma}^{w}_{\text{ts}}$ over $p_{d}\geq 0.15$ were calculated over the estimates with the maximal $|t|_{\rm min}$ (and hence three degrees of freedom). The FSS estimates $\nu^{\rm g}$ and $\gamma^{\rm g}$ are taken from [21] and $\overline{\nu}^{w}_{\rm lin}$ and $\overline{\nu}^{w}$ from [29].  $a$ | $\overline{a}$ | $\overline{\nu}^{w}_{\text{ts}}$ | $\overline{\nu}^{w}_{\rm lin}$ | $\overline{\nu}^{w}$ | $\nu^{\rm g}$ | $\overline{\gamma}^{w}_{\text{ts}}$ | $\gamma^{\rm g}$
---|---|---|---|---|---|---|---
$\infty$ | $\infty$ | 0.6928(17) | 0.6913(15) | 0.6843(31) | 0.6831(30) | 1.3430(18) | 1.3324(64)
3.5 | 3.30(18) | 0.7557(25) | 0.7427(25) | 0.7122(49) | 0.7117(49) | 1.4875(33) | 1.451(15)
3.0 | 2.910(96) | 0.7898(34) | 0.7812(35) | 0.7532(53) | 0.7484(52) | 1.5726(50) | 1.566(16)
2.5 | 2.451(26) | 0.8905(82) | 0.8887(61) | 0.8725(96) | 0.8719(96) | 1.787(11) | 1.783(24)
2.0 | 1.979(18) | 1.073(23) | 1.079(14) | 1.067(23) | 1.060(23) | 2.171(27) | 2.149(51)
1.5 | 1.500(30) | 1.348(61) | 1.449(32) | 1.435(56) | 1.421(55) | 2.791(70) | 2.93(14)

(a) $a=\infty$ (b) $a=3.5$ (c) $a=3.0$ (d) $a=2.5$ (e) $a=2.0$ (f) $a=1.5$ (g)

Figure 6: (Colour online) Final results of the fits of $[\tilde{\chi}(t)]$
with the ansatz (3.2) for different $|t|_{\rm min}$ and all defect
concentrations $p_{d}$. The weighted means $\overline{\gamma}^{w}_{\text{ts}}$
over all concentrations with $p_{d}\geqslant 0.15$ and for the largest
possible $|t|_{\rm min}$ are shown together with the result from the FSS
analysis $\gamma^{\rm g}$. The narrow right-hand panels show a separate
comparison between the different estimates for $\gamma$ which are plotted as
horizontal lines in the main plots.

##  4 Conclusions

The main merit of temperature scaling is its conceptual simplicity. It
provides direct estimates of the critical exponents $\nu$, $\gamma$, …,
whereas in the complementary FSS approach they can only be computed from the
fitted exponent ratios $1/\nu$, $\gamma/\nu$, …which requires some care with
statistical error propagation. A drawback of temperature scaling is that the
determination of a suitable fit interval requires care at both ends. If the
included temperatures $T$ are too far away from $T_{c}$, corrections-to-
scaling cannot be neglected. At the other end, the included $T$ values should
not be too close to $T_{c}$, because finite-size effects become important. By
contrast, in the FSS approach, only the lower end of the fit interval needs to
be controlled: The minimal lattice size $L$ must be large enough to avoid
sizeable corrections-to-scaling. In this case, in FSS analyses one generically
deals with simple linear two-parameter fits (since $T_{c}$ only enters
indirectly), whereas in temperature scaling, this is only possible when
$T_{c}$ is known from other sources (otherwise more cumbersome non-linear
three-parameter fits are necessary). Of course, in both approaches, the
situation becomes more complicated when corrections-to-scaling should be
included since this introduces additional parameters and generically requires
a non-linear many-parameter fitting.

Here, we successfully used the temperature-scaling analysis to validate our
FSS results [21] for the critical exponents $\nu$ and $\gamma$. Considering
that mostly different data entered the analysis and also different observables
were used, i.e., the correlation length $\xi$ that was not studied in the FSS
approach and the high-temperature definition of the susceptibility, we can
clearly solidify our FSS results. Additionally, the critical temperatures
estimated in the FSS study [21] were confirmed to be reasonably accurate to be
used in the temperature-scaling analysis.

The available data turned out, however, to be not sufficiently accurate to
perform fits including corrections-to-scaling, and the uncorrected fits show a
clear dependence on the temperature range used. In order to improve the
temperature-scaling analysis, we would need more simulated temperatures and
also change the setup of the entire simulation process by using more disorder
realizations instead of longer measurement time series for each realization,
which is planned for future work.

## Acknowledgements

This paper is dedicated to Professor Bertrand Berche on the occasion of his
60th birthday, with whom the basics of the present work was laid 20 years ago
in a series of joint papers on uncorrelated bond disorder. It has been always
a pleasure to collaborate with Bertrand — and to enjoy the after-work
activities!

We thank the Max Planck Society, the Max Planck Institute for Mathematics in
the Sciences (MIS), and especially the Graduate School IMPRS-MIS for financial
support of SK and providing the computational resources at the Max Planck
Computing and Data Facility in Munich.

We also gratefully acknowledge further support by the Deutsch-Französische
Hochschule (DFH-UFA) through the Doctoral College “$\mathbb{L}^{4}$” under
Grant No. CDFA-02-07 which is co-directed by Bertrand. Many thanks go to him,
Malte Henkel, Dragi Karevski, and Christophe Chatelain for many fruitful
discussions and hosting SK’s visit in Nancy. We also thank Yurko Holovatch,
Viktoria Blavatska, and Mikhail Nalimov for interesting discussions that
contributed to a deeper insight into the topic.

## References

  * [1] Harris A. B., J. Phys. C: Solid State Phys., 1974, 7, 1671–1692, doi:10.1088/0022-3719/7/9/009.
  * [2] Ballesteros H. G., Fernández L. A., Martín-Mayor V., Muñoz Sudupe A., Parisi G., Ruiz-Lorenzo J. J., Phys. Rev. B, 1998, 58, 2740–2747, doi:10.1103/PhysRevB.58.2740.
  * [3] Calabrese P., Martín-Mayor V., Pelissetto A., Vicari E., Phys. Rev. E, 2003, 68, 036136 (17 pages), doi:10.1103/PhysRevE.68.036136.
  * [4] Berche P.-E., Chatelain C., Berche B., Janke W., Comput. Phys. Commun., 2002, 147, 427–430, doi:10.1016/S0010-4655(02)00319-3.
  * [5] Berche P.-E., Chatelain C., Berche B., Janke W., Monte Carlo Studies of Three-Dimensional Bond-Diluted Ferromagnets, In: High Performance Computing in Science and Engineering, Munich 2002, Wagner S., Hanke W., Bode A., Durst F. (Eds.), Springer, Berlin, 2003, 227–238, Preprint arXiv:cond-mat/0212504.
  * [6] Janke W., Berche P.-E., Chatelain C., Berche B., Phase Transitions in Disordered Ferromagnets, In: NIC Symposium 2004, Proceedings, Wolf D., Münster G., Kremer M. (Eds.), John von Neumann Institute for Computing, Jülich, NIC Series, 2003, Vol. 20, 241–250.
  * [7] Janke W., Berche P.-E., Chatelain C., Berche B., Quenched Disorder Distributions in Three-Dimensional Diluted Ferromagnets, In: Computer Simulation Studies in Condensed-Matter Physics XVI, Landau D. P., Lewis S. P., Schüttler H.-B. (Eds.), Springer, Berlin, 2004, 89–94, Preprint arXiv:cond-mat/0304642.
  * [8] Berche P.-E., Chatelain C., Berche B., Janke W., Eur. Phys. J. B, 2004, 38, 463–474,
doi:10.1140/epjb/e2004-00141-x.

  * [9] Berche B., Berche P.-E., Chatelain C., Janke W., Condens. Matter Phys., 2005, 8, 47–58, doi:10.5488/CMP.8.1.47.
  * [10] Janke W., Berche B., Chatelain C., Berche P.-E., Hellmund M., PoS, 2005, LAT2005, 018 (22 pages),
doi:10.22323/1.020.0018.

  * [11] Murtazaev A. K., Kamilov I. K., Babaev A. B., J. Exp. Theor. Phys., 2004, 99, 1201–1206,
doi:10.1134/1.1854807.

  * [12] Hasenbusch M., Toldin F. P., Pelissetto A., Vicari E., J. Stat. Mech., 2007, 2007, No. 02, P02016 (43 pages), doi:10.1088/1742-5468/2007/02/P02016.
  * [13] Kos F., Poland D., Simmons-Duffin D., Vichi A., J. High Energ. Phys., 2016, 2016, No. 8, 36 (15 pages), doi:10.1007/JHEP08(2016)036.
  * [14] Ferrenberg A. M., Xu J., Landau D. P., Phys. Rev. E, 2018, 97, 043301 (12 pages),
doi:10.1103/PhysRevE.97.043301.

  * [15] Weinrib A., Halperin B. I., Phys. Rev. B, 1983, 27, 413–427, doi:10.1103/PhysRevB.27.413.
  * [16] Honkonen J., Nalimov M. Y., J. Phys. A: Math. Gen., 1989, 22, 751–763, doi:10.1088/0305-4470/22/6/024.
  * [17] Korzhenevskii A. L., Luzhkov A. A., Schirmacher W., Phys. Rev. B, 1994, 50, 3661–3666,
doi:10.1103/PhysRevB.50.3661.

  * [18] Korzhenevskii A. L., Luzhkov A. A., Heuer H.-O., Europhys. Lett., 1995, 32, 19–24,
doi:10.1209/0295-5075/32/1/004.

  * [19] Dudka M., Fedorenko A. A., Blavatska V., Holovatch, Yu., Phys. Rev. B, 2016, 93, 224422 (13 pages), doi:10.1103/PhysRevB.93.224422.
  * [20] Kazmin S., Janke W., Phys. Rev. B, 2020, 102, 174206 (13 pages), doi:10.1103/PhysRevB.102.174206.
  * [21] Kazmin S., Janke W., Phys. Rev. B, 2022, 105, 214111 (12 pages), doi:10.1103/PhysRevB.105.214111.
  * [22] Makse H., Havlin S., Stanley H. E., Schwartz M., Chaos, Solitons & Fractals, 1995, 6, 295–303, doi:10.1016/0960-0779(95)80035-F.
  * [23] Makse H. A., Havlin S., Schwartz M., Stanley H. E., Phys. Rev. E, 1996, 53, 5445–5449,
doi:10.1103/PhysRevE.53.5445.

  * [24] Zierenberg J., Fricke N., Marenz M., Spitzner F. P., Blavatska V., Janke W., Phys. Rev. E, 2017, 96,
062125 (11 pages), doi:10.1103/PhysRevE.96.062125.

  * [25] Fricke N., Zierenberg J., Marenz M., Spitzner F. P., Blavatska V., Janke W., Condens. Matter Phys., 2017, 20, 13004 (11 pages), doi:10.5488/CMP.20.13004.
  * [26] Swendsen R. H., Wang J.-S., Phys. Rev. Lett., 1987, 58, 86–88, doi:10.1103/PhysRevLett.58.86.
  * [27] Janke W., In: Computational Many-Particle Physics, Lecture Notes in Physics, Vol. 739, Fehske H., Schneider R., Weiße A. (Eds.), Springer, Berlin Heidelberg, 2008, 79–140, doi:10.1007/978-3-540-74686-7_4.
  * [28] Efron B., Tibshirani R., An Introduction to the Bootstrap, Monographs on Statistics and Applied Probability, Chapman & Hall, New York, 1993.
  * [29] Kazmin S., Ising Model in Three Dimensions with Long-Range Power-Law Correlated Site Disorder: A Monte Carlo Study, PhD Thesis, Universität Leipzig, May 2021,
URL https://nbn-resolving.org/urn:nbn:de:bsz:15-qucosa2-798018.

  * [30] Blavats’ka V., von Ferber C., Holovatch Y., Phys. Rev. E, 2001, 64, 041102 (10 pages),
doi:10.1103/PhysRevE.64.041102.

  * [31] Ballesteros H. G., Parisi G., Phys. Rev. B, 1999, 60, 12912–12917, doi:10.1103/PhysRevB.60.12912.
  * [32] Prudnikov V. V., Prudnikov P. V., Fedorenko A. A., Phys. Rev. B, 2000, 62, 8777–8786,
doi:10.1103/PhysRevB.62.8777.

  * [33] Prudnikov V. V., Prudnikov P. V., Dorofeev S. V., Kolesnikov V. Y., Condens. Matter Phys., 2005, 8, 213–224, doi:10.5488/CMP.8.1.213.
  * [34] Ivaneyko D., Berche B., Holovatch Yu., Ilnytskyi J., Physica A, 2008, 387, 4497–4512,
doi:10.1016/j.physa.2008.03.034.

  * [35] Wang W., Meier H., Lidmar J., Wallin M., Phys. Rev. B, 2019, 100, 144204 (8 pages),
doi:10.1103/PhysRevB.100.144204.

Ukrainian 0

[Ñòåïåíåâà cêîðåëüîâàíà íåâïîðÿäêîâàíà 3D ìîäåëü çíãà]Àíàëç òåìïåðàòóðíîãî
ñêåéëíãó òðèâèìðíî¿ íåâïîðÿäêîâàíî¿ ìîäåë çíãà ç ñòåïåíåâèìè ñêîðåëüîâàíèìè
äåôåêòàìè Ñ. Êàçìí${}^{\textsf{\footnotesize{\mbox{?},\penalty
1000\mbox{?}}}}$, Â. ßíêå${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ íñòèòóò òåîðåòèíî¿ ôçèêè, Óíâåðñèòåò Ëåéïöèãà,
IPF 231101, 04801 Ëåéïöèã, Íìåèíà ${}^{\textsf{\footnotesize 2}}$ Íìåöüêèé
íåêîìåðöéíèé íàóêîâî-äîñëäíèé öåíòð áîëîãíèõ ìàòåðàëâ, Òîðãàóåð øòð. 116,
04347 Ëåéïöèã, Íìåèíà


# Prediction of MoO2 as high capacity electrode material for (Na, K, Ca)$-$ion
batteries11footnotemark: 1

Yong-Chao Rao, Song Yu, Xiao Gu  Xiang-Mei Duan∗ Department of Physics,
Faculty of Science, Ningbo University, Ningbo-315211, P.R. China
<EMAIL_ADDRESS>
###### Abstract

Na, K and Ca$-$ion battery electrode materials with appropriate
electrochemical properties are desirable candidates for replacing
lithium$-$ion batteries (LIBs) because of their natural richness and low cost.
Recently, MoO2 has been reported as the anode material in LIBs, but so far not
received attention in Na and other ion batteries. In this paper, the behaviors
of Na, K, and Ca on MoO2 are investigated by first$-$principles calculations.
These metal atoms strongly absorb on the hexagonal center of MoO2 and the
adsorption results in semiconducting$-$metallic transition. The low diffusion
barrier, 0.13, 0.08, 0.22 eV for Na, K, Ca, respectively, leads to an
ultrahigh diffusivity. Importantly, the maximum metal$-$storage phases of MoO2
monolayer correspond to Na4MoO2, K2.5MoO2 and Ca3MoO2, with considerable
theoretical specific capacities of 837, 523 and 1256 mAh g-1. The electrode
materials exhibit moderate average voltage of 0.30, 0.75, 0.35 V,
respectively. Our findings suggest that MoO2 monolayer can be utilized as a
promising anode material with high capacities and high rate performance for
next generation ion batteries.

###### keywords:

MoO2 monolayer; Ion$-$battery; First principles calculations

††journal: Journal of LaTeX Templates

##  1 Introduction

With the tremendous and rapid development of portable electronic devices and
electric vehicles, rechargeable batteries with large capacity, high rate
capability and good cycle performance have attracted great attention both in
basic research and industry applications.[1] As one of the most important
energy storage devices, lithium$-$ion batteries (LIBs) are widely used in the
commercial society.[2] However, issues such as cost, safety, energy storage
density, charge/discharge rates, and lifetime will continue to spur the
research for new batteries beyond LIBs.[3, 4] Room temperature sodium$-$ion
batteries (SIBs) have gained increasing attention due to the cheaper cost and
relative abundance of sodium.[5, 6] Calcium$-$ion batteries (CIBs) offer low
price, chemical safety and lighter mass$-$to$-$charge ratio.[7] Moreover,
potassium is indeed a competitive metal with less complicated interfacial
reaction [8] and higher ionic conductivity [9] in solution.

The battery$-$performance characteristic strongly depends on the
electrochemical properties of electrode materials. Great challenge of advanced
metal$-$ion batteries is to find better electrode substrates. Currently, a
wide range of compounds have been proposed for SIB anodes, and much progress
has been made .[10, 5, 11] However, capacity for SIBs are still limited owing
to the relatively large ionic radius of Na$-$ion (1.02 Å) compared to Li$-$ion
(0.76 Å). The commercially used anode material in LIBs, graphite, is not
suitable for Na$-$ion batteries for the low capacity and poor rate
capability.[12] With the high specific surface area, 2D materials provide
variable choice in designing high energy and high motility devices.[13, 14]
Much research efforts turn to the development of SIBs, KIBs (potassium
ion$-$battery) and CIBs based on 2D materials, such as graphene[15] and its
dopant,[16, 17] transition metal oxides,[18, 19] TMDs,[20, 21] silicene[22]
and other van der Waals$-$bonded layered materials.[10, 23] Recently, MoO2
sub$-$micro sheet was synthesized by using a novel chemical vapour deposition
method.[24] The dynamical and thermal stabilities MoO2 monolayer have been
confirmed via ab initio molecular dynamics simulation and phonon dispersion
spectrum. Moreover, Li$-$binding on MoO2 sheet possesses superb conductivity,
lower diffusion barrier (75 meV), and astonishing theoretical storage capacity
(2513 mAh g-1).[25]

Due to the electrochemical similarity between LIBs and SIBs, KIBs and CIBs, an
interesting question arises: Can the superior electrochemical property of the
MoO2 sheet, demonstrated in LIBs, extend into Na, K, Ca$-$ion cells? In this
work, DFT calculations are employed to systematically investigate the
electronic, metal adsorption and storage behaviors of MoO2 monolayer. Our
results indicate that Na, K and Ca atoms can strongly absorb on the surface of
MoO2, and the metallic characteristics they presented, are crucial for the
anode materials. The diffusion barrier for Na, K, Ca on MoO2 is only 0.13,
0.08, 0.23 eV, respectively, implying good diffusion mobility during the
charging/discharging process. Furthermore, the theoretical specific capacities
are as high as 837, 523, 1256 mAh g-1 for Na, K, Ca$-$ion batteries. Most
especially, the capacity for Na is 4.5 times higher than that of the oxidized
graphite electrode (184 mAh g-1).[26] Therefore, as electrode materials for
Na, K, Ca$-$ion batteries, the MoO2 sheet possesses great potential and
outstanding performance.

##  2 Computational details

DFT calculations were performed by using the Vienna ab initio Simulation
Package (VASP). The exchange$-$correction interaction were treated within the
generalized gradient approximation (GGA) in the form of the
Perdew$-$Burke$-$Ernzerhof (PBE) functional.[27] The electron wave functions
were expanded using plane waves with a cutoff energy of 500 eV, and the
convergence criteria for the residual force and energy on each atom during
structure relaxation was set to 10-2 eV/Å and 10-5 eV, respectively. $4\times
4\times 1$ supercells were built and complemented by vacuum layers of 20 Å
thickness in the out$-$of$-$plane direction to avoid interaction between
periodic images. The $k-$point sampling of $3\times 3\times 1$ and $7\times
7\times 1$ were chosen for geometry optimization and electronic structure
computation, respectively. For the adatom diffusion on MoO2 sheet, the
climbing image nudge elastic band method was used to perform minimum energy
path profiling and figure out the diffusion barrier.[28]

The strength between the metal atoms and the substrate was described by the
average adsorption energy, $E_{\rm ads}$, which is defined as

$E_{\rm ads}=(E_{{{\rm A}_{x}\rm{MoO}_{\rm 2}}}-E_{\rm{MoO_{\rm 2}}}-x\mu_{\rm
A})/x$ (1)

where $x$ is the total number of adsorbed metal atoms in the primitive cell
and $E_{{{\rm A}_{x}\rm{MoO}_{\rm 2}}}$, $E_{\rm{MoO_{\rm 2}}}$, $\mu_{\rm A}$
is total energies of the metal ion adsorbed on the system, pristine MoO2
monolayer, and the energy per atom in the bulk metal, respectively.

The average voltage of AxMoO2 (A $=$ Na, K, and Ca) in the range of $x_{\rm
1}\leq x\leq x_{\rm 2}$ is given as

$V_{\rm ave}\approx[E_{{{\rm A}_{x_{\rm 1}}\rm{MoO}_{\rm 2}}}-E_{{{\rm
A}_{x_{\rm 2}}\rm{MoO}_{\rm 2}}}+(x_{\rm 2}-x_{\rm 1})\mu_{\rm A}]/(x_{\rm
2}-x_{\rm 1})e$ (2)

According to the Faraday equation, the maximum specific capacity, $C_{\rm A}$,
is obtained by[16]

$C_{\rm A}=x_{\rm max}F/M$ (3)

Where $x_{\rm max}$ represents maximum number of A in the primitive cell, $F$
is the Faraday constant with the value of 26806 mAh mol-1, and $M$ is the mass
of MoO2 in g mol-1.

##  3 Results and discussion

###  3.1 Adsorption of Na, K, Ca on MoO2 surface

As the common transition$-$metal dichalcohenides (TMDs), MoO2 monolayer is
composed of a triple layer (sandwich structure), where atoms are stacked in a
sequence of O$-$Mo$-$O with thickness of 2.47 Å, the optimized lattice
parameter is found to be $a=b=2.83$ Å, which are in good agreement with the
previous calculation values of $2.82$ Å. [29]

To systematically study metal atoms intercalation on MoO2 monolayer, we first
examine four possible high$-$symmetric adsorption sites (H site: hexagonal
center site; B site: over the center of the Mo$-$O bond; T site: the top of a
O atom; T′ site: the top of a Mo atom) for isolated Na, K, Ca atom, as plotted
in Fig. 1 (a). The adsorption energies ($E_{\rm ads}$) of a single Na, K, and
Ca atom adsorbed on MoO2 are illustrated in Fig. 1 (c). In general, the
K$-$adsorbed configurations are the most stable, while the strength of
adsorption of one Ca atom on MoO2 is the weakest with the smallest $E_{\rm
ads}$. The negative adsorption energies for H, B and T′ configurations
indicate the Na, K, and Ca storage is exothermic and spontaneous, which is
foundational and favorable for the application. H site is the most stable
configuration with the largest adsorption energy of $-$1.46, $-$2.01, $-$1.16
eV for Na, K and Ca, respectively. We note that a single Na atom adsorbed on
the MoO2 monolayer is more stable, when compared with other 2D materials, such
as boron$-$doped graphene ($-$0.79 eV)[30] and single$-$layer MoS2 ($-$1.23
eV).[31] High adsorption strength ensures that metal atoms do not form
clusters, which is beneficial to keep good cycling during the
charging/discharging process of batteries. More particularly, the adsorption
energies on T′ site ($-$1.42, $-$1.99, $-$1.07 eV for Na, K and Ca,
respectively) are slightly smaller that that of H site, which implies that the
metal atoms may move from the H site to the T′ site with lower energy barrier.

Conductivity is a key factor in determining the electrochemical performance of
an electrode. Pristine MoO2 is an indirect semiconductor with a band gap of
0.92 eV, which is consistent with previous theoretical work of 0.91 eV,[29]
implying poor electrical conductivity. The valence band maximum (VBM) of MoO2
is contributed by O$-2p$ orbitals and its conduction band minimum (CBM) is
dominated by Mo$-4d$ orbitals. The total density of states (TDOS) and
projected density of states (PDOS) for the adsorption of single Na, K, Ca,
respectively, are shown in Fig. 2 (a), (b) and (c). It can be seen that three
adsorbed systems exhibit metallic characteristic, the CBM of each system is
still dominated by Mo 4$d$ orbital, and the $s$ orbital of metal atoms is
empty. To obtain insights into the nature of metals$-$MoO2 binding, we refer
to the charge density differences of (110) section defined as
$\Delta\rho=\rho_{\rm metal/MoO_{2}}-\rho_{\rm metal}-\rho_{\rm MoO_{2}}$,
which are shown in Fig. 3. Clearly, electrons accumulate around oxygen atoms
and their densities surrounding metal atom decrease, suggesting ionic
character of metals$-$MoO2 binding. By using Bader charge analysis, the
absorbed Na, K and Ca atoms act as electron donor, and transfer 0.9, 0.93 and
1.61 $e$ to the underlying MoO2, respectively.

The spin$-$dependent GGA plus Hubbard U (GGA$+$U) method was employed as a
test. Comparing the $E_{\rm ads}$ and DOS with the results obtained from the
common GGA method, the difference is negligible, implying the weak onsite
Coulomb repulsion between the 4$d$ electrons of molybdenum. Moreover, the
systems with/without metals absorption show non$-$magnetic behavior.

###  3.2 Diffusion of metal adatoms on MoO2 surface

Another factor that plays a central role in the performance of electrode
materials is the diffusion rate of ions, which contributes directly to the
charging/discharging loading rates accessible by a battery. We investigate the
diffusion pathways and energy barriers of Na, K and Ca on MoO2 surface. Only
the diffusion path between two neighboring, most preferred binding sites are
taken into consideration. Due to the symmetry, there are three representative
pathways for metal atom migration between two adjacent H sites: (a) through
the T site (H$-$T$-$H), (b) across the Mo$-$O bond (H$-$B$-$H), and (c)
through the T′ site. As we mentioned above, the almost equivalent adsorption
energy at H and T′ sites motivates us to only consider the pathway (c). The
optimized route is shown in Fig. 4(a). The corresponding energy profiles of
Na, K and Ca along pathway on the MoO2 surface are depicted in Fig. 4(b). The
barrier for Na is 0.13 eV, implying very fast Na diffusion kinetics, which
would lead MoO2 to be an effective anode electrode for Na$-$ion batteries. The
value of the barrier is comparable with other 2D materials such as MoS2 (0.11
eV),[31] graphene (0.13 eV)[16] and boron$-$doped graphene (0.16$-$0.22
eV).[30] Particularly, the barrier in K$-$adsorbed system is as low as 0.08
eV, smaller than corresponding Na system, which can be ascribed to the larger
vertical distance (2.16 Å) between K and substrate, while Na$-$substrate
distance is 1.68 Å. Similar trends have been observed in VS2[32] and MoN2.[33]
For the case of Ca, the relative large mass results in a short
cation$-$substrate distance of 1.50 Å, the diffusion barrier of Ca is 0.22 eV.
Therefore, we conclude that MoO2 monolayer have the ability to rapidly
charge/discharge Na and K ions.

###  3.3 Average open$-$circuit voltage and metals storage capacity of the
MoO2 monolayer

From the practical point of view, the average open$-$circuit voltage and metal
storage capacity are significant features for the electrode materials. Due to
the high symmetry of MoO2 structure, we consider all the adsorption with
double sides to obtain the higher capacity. The typical NaxMoO2 configurations
are illustrated in Fig. 5. To avoid the inherently strong coulomb repulsion as
metals absorb on the top of the MoO2 surface, the adsorption structures appear
to be mixed configurations. As shown in Fig. 5 (a), four Na atoms absorb on
the H site due to its relatively larger $E_{\rm ads}$. The additional four Na
atoms tend to adsorb on the next nearest neighbor T′ site, as shown in Fig.
5(b), to reduce the repulsion as possible. As the concentration increases, the
physical attraction between the further far metal atoms and MoO2 becomes
weaker, so the adsorption prefers to proceed layer by layer, with each layer
containing eight metal atoms [see Na3MoO2 in Fig. 5 (c)].

In Fig. 6, we present the OCVs as function of the concentration by varying $x$
in the AxMoO2 system. With the increase of $x$, the OCVs decrease gradually
because the distances between adjacent metal atoms become shorter and their
repulsive interaction is enhanced. Once the OCV is zero, no more atom can be
adsorbed, and the maximum $x$ in AxMoO2 is obtained. Encouragingly, the large
initial OCVs of 1.27, 1.79, and 1.16 V for Na, K, and Ca systems,
respectively, indicate the potential application of MoO2 in ion$-$batteries.
For the NaxMoO2 system, the electrode potential varies in the range of
1.27$-$0.09 V. The $V_{\rm ave}$, 0.30 V, is in the desirable voltage range
for a sodium anode of 0.0$-$1.0 V. At $x=4$, the corresponding maximum
specific capacity of MoO2 for SIBs is 837 mAh g-1, which is larger than that
of existing SIB anode materials (including hard carbon,[34, 35] Na4Ti5O12,[36]
Na2Ti3O7[37]) and is even larger than some alloy$-$based materials.[38, 39]
Compared with the capacities of 2D materials, such as MoS2/graphene of 338 mAh
g-1,[23] VS2 of 232.91 and 116.45 mAh g-1 for 1H and 1T phases,
respectively[40] and GeS nano$-$sheet of 512 mAh g-1,[41] the capacity of
Na4MoO2 is also competitive. Therefore MoO2 monolayer should be a promising
candidate for an electrode material in SIBs. For the case of K and Ca adsorbed
MoO2, electrode potentials are in the range of 0.02$-$1.79 V and 0.25$-$1.16
V, with the $V_{\rm ave}$ of 0.75 and 0.35 V, respectively. Due to the large
ionic radius of K and Ca, the maximum $x$ is estimated to be 2.5 and 3, and
the corresponding capacity of KIB and CIB is 523 and 1256 mAh g-1.

When evaluating cycling stability, the degree of structural deformation would
cause substantial loss of capacity in several cycles. After the absorption of
maximum Na, K, and Ca, the thickness of MoO2 is 2.50, 2.48, 2.52 Å
respectively, slightly stretched from the original thickness of 2.47 Å. The
bond length of Mo$-$O changes negligibly. During the process of ions
intercalation/deintercalation, the structure remains intact to ensure good
cycle performance and a smooth surface that will undoubtedly minimizes
diffusion resistance during the metal ion transport. So far, the MoO2 anode
for KIBs is predicted to be a perfect material, which simultaneously possesses
a large specific capacity (523 mAh g-1), excellent stability and conductivity
[see Fig. 2(b)], and a modest open circuit voltage (0.75 V).

##  4 Conclusion

In summary, we investigate the properties of 2D MoO2 monolayer as potential
anodes for batteries beyond LIBs. With large adsorption energies, Na, K or Ca
atom is preferentially absorbed over the center of the honeycomb lattice, and
the charge transfer from metal atoms to the substrate results in a transition
of semicondutive MoO2 to metallic AxMoO2. The diffusion barriers of metal
atoms on the surface of MoO2 are relatively low, especially for Na and K,
which are only 0.13 and 0.08 eV. The theoretical specific capacity of Na4MoO2,
K2.5MoO2 and Ca3MoO2 are as high as 837 mAh g-1, 523 mAh g-1 and 1256 mAh g-1,
respectively. When the MoO2 anodes are used for SIBs, KIBs and CIBs, the open
voltage is 0.30 V, 0.75 V and 0.35 V, respectively. Moreover, MoO2 exhibits
high mechanical stability and integrity upon metal insertion. Our results give
insights for further experimental work in exploring the possibility of MoO2
for Na, K and Ca$-$ion batteries.

## Acknowledgments

This research is supported by the Natural Science Foundation of China (grant
No. 11574167), the New Century 151 Talents Project of Zhejiang Province and
the KC Wong Magna Foundation in Ningbo University.

## Conflict of interest

The authors declare they have no conflict of interest

## References

## References

  * [1] B. Dunn, H. Kamath, J. M. Tarascon, Electrical energy storage for the grid: a battery of choices, Science 334 (2011) 928–935.  doi:10.1126/science.1212741.
  * [2] A. Manthiram, Y. Fu, Y. S. Su, In charge of the world: electrochemical energy storage, J. Phys. Chem. Lett. 4 (2013) 1295–1297.  doi:10.1021/jz4006652.
  * [3] J. B. Goodenough, K. S. Park, The Li-ion rechargeable battery: a perspective, J. Am. Chem. Soc. 135 (2013) 1167–1176.  doi:10.1021/ja3091438.
  * [4] J. M. Tarascon, Is lithium the new gold ?, Nat. Chem. 2 (2010) 510.  doi:10.1038/nchem.680.
  * [5] S. W. Kim, D. H. Seo, X. Ma, G. Ceder, K. Kang, Electrode materials for rechargeable sodium-ion batteries: potential alternatives to current lithium-ion batteries, Adv. Energy Mater. 2 (2012) 710–721.  doi:10.1002/aenm.201200026.
  * [6] V. Palomares, P. Serras, I. Villaluenga, K. B. Hueso, J. Carreterogonzaĺez, T. Rojo, Na-ion batteries, recent advances and present challenges to become low cost energy storage systems, Energy Environ. Sci. 5 (2012) 5884–5901.  doi:10.1039/C2EE02781J.
  * [7] M. Hayashi, H. Arai, H. Ohtsuka, Y. Sakurai, Electrochemical characteristics of calcium in organic electrolyte solutions and vanadium oxides as calcium hosts, J. Power Sources 119 (2003) 617–620.  doi:10.1016/S0378-7753(03)00307-0.
  * [8] M. Moshkovich, Y. Gofer, D. Aurbach, M. Moshkovich, Y. Gofer, D. Aurbach, Investigation of the electrochemical windows of aprotic alkali metal (Li, Na, K) salt solutions, J. Electrochem. Soc. 148 (2001) 7282–7289.  doi:10.1149/1.1357316.
  * [9] A. Eftekhari, Z. Jian, X. Ji, Potassium secondary batteries, ACS Appl. Mater. Interfaces 9 (2017) 4404–4419.  doi:10.1021/acsami.6b07989.
  * [10] D. Su, S. Dou, G. Wang, WS2@graphene nanocomposites as anode materials for Na-ion batteries with enhanced electrochemical performances, Chem. Commun. 50 (2014) 4192–4195.  doi:10.1039/C4CC00840E.
  * [11] C. S. Park, H. Kim, R. A. Shakoor, E. Yang, S. Y. Lim, R. Kahraman, Y. Jung, J. W. Choi, Anomalous manganese activation of a pyrophosphate cathode in sodium ion batteries: a combined experimental and theoretical study, J. Am. Chem. Soc. 135 (2013) 2787–2792.  doi:10.1021/ja312044k.
  * [12] M. Winter, J. O. Besenhard, M. E. Spahr, P. Novák, Insertion electrode materials for rechargeable lithium batteries, Adv. Mater. 10 (2010) 725–763.  doi:10.1002/(SICI)1521-4095(199807)10:10<725::AID-ADMA725>3.0.CO;2-Z.
  * [13] E. Yoo, J. Kim, E. Hosono, H. S. Zhou, T. Kudo, I. Honma, Large reversible Li storage of graphene nanosheet families for use in rechargeable lithium ion batteries, Nano Lett. 8 (2008) 2277–2282.  doi:10.1021/nl800957b.
  * [14] S. Yu, Y. C. Rao, S. F. Li, X. M. Duan, Net W monolayer: a high-performance electrode material for Li-ion batteries, Appl. Phys. Lett. 112 (2018) 053903\.  doi:10.1063/1.5013617.
  * [15] O. I. Malyi, K. Sopiha, V. V. Kulish, T. L. Tan, S. Manzhos, C. Persson, A computational study of Na behavior on graphene, Appl. Surf. Sci. 333 (2015) 235–243.  doi:10.1016/j.apsusc.2015.01.236.
  * [16] D. Datta, J. Li, V. B. Shenoy, Defective graphene as a high-capacity anode material for Na- and Ca-ion batteries, ACS Appl. Mater. Interfaces 6 (2014) 1788–1795.  doi:10.1021/am404788e.
  * [17] S. Gong, Q. Wang, Boron-doped graphene as a promising anode material for potassium-ion batteries with a large capacity, high rate performance and good cycling stability, J. Phys. Chem. C 121 (2017) 24418–24424.  doi:10.1021/acs.jpcc.7b07583.
  * [18] A. Caballero, L. Hernán, J. Morales, L. Sánchez, J. Santos Pen̈a, M. A. G. Aranda, Synthesis and characterization of high-temperature hexagonal P2-Na0.6MnO2 and its electrochemical behaviour as cathode in sodium cells, J. Mater. Chem. 12 (2002) 1142–1147.  doi:10.1039/B108830K.
  * [19] H. Liu, H. Zhou, L. Chen, Z. Tang, W. Yang, Electrochemical insertion/deinsertion of sodium on NaV6O15 nanorods as cathode material of rechargeable sodium-based batteries, J. Power. Sources 196 (2011) 814–819.  doi:10.1016/j.jpowsour.2010.07.062.
  * [20] J. Shuai, H. Deog Yoo, Y. Liang, Y. Li, Y. Yao, L. C. Grabow, Density functional theory study of Li, Na, and Mg intercalation and diffusion in MoS2 with controlled interlayer spacing, Mater. Res. Express 3 (2016) 064001.  doi:10.1088/2053-1591/3/6/064001.
  * [21] E. Yang, H. Ji, Y. Jung, Two-dimensional transition metal dichalcogenide monolayers as promising sodium ion battery anodes, J. Phys. Chem. C 119 (2015) 26374–26380.  doi:10.1021/acs.jpcc.5b09935.
  * [22] J. Zhu, U. Schwingenschlögl, Silicene for Na-ion battery applications, 2D Mater. 3 (2016) 035012.  doi:10.1088/2053-1583/3/3/035012.
  * [23] L. David, R. Bhandavat, G. Singh, MoS2/Graphene composite paper for sodium-ion battery electrodes, ACS Nano 8 (2014) 1759–1770.  doi:10.1021/nn406156b.
  * [24] X. Liu, Y. He, S. Wang, Q. Zhang, Preparation of MoO2 sub-micro scale sheets and their optical properties, J. Alloy Comp. 509 (2011) S408–S411.  doi:10.1016/j.jallcom.2011.01.089.
  * [25] Y. Zhou, C. Geng, A MoO2 sheet as promising electrode material: ultrafast Li-diffusion and astonishing Li-storage capacity, Nanotech. 28 (2017) 105402.  doi:10.1088/1361-6528/aa56d0.
  * [26] W. Yang, H. Kai, Y. Zhu, F. Han, Y. Xu, I. Matsuda, Y. Ishii, J. Cumings, C. Wang, Expanded graphite as superior anode for sodium-ion batteries, Nat. Commun. 5 (2014) 4033.  doi:10.1038/ncomms5033.
  * [27] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 78 (1996) 3865–3868.  doi:10.1103/PhysRevLett.77.3865.
  * [28] G. Henkelman, H. Jńsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys. 113 (2000) 9978–9985.  doi:10.1063/1.1323224.
  * [29] F. A. Rasmussen, K. S. Thygesen, Computational 2D materials database: electronic structure of transition-metal dichalcogenides and oxides, J. Phys. Chem. C 119 (2015) 13169–13183.  doi:10.1021/acs.jpcc.5b02950.
  * [30] C. Ling, F. Mizuno, Boron-doped graphene as a promising anode for Na-ion batteries, Phys. Chem. Chem. Phys. 16 (2014) 10419–10424.  doi:10.1039/C4CP01045K.
  * [31] J. Su, Y. Pei, Z. Yang, X. Wang, Ab initio study of graphene-like monolayer molybdenum disulfide as a promising anode material for rechargeable sodium ion batteries, RSC Adv. 4 (2014) 43183–43188.  doi:10.1039/C4RA06557C.
  * [32] M. J. Momeni, E. Targholi, M. Mousavi-Khoshdel, Maleic anhydride as a promising anode material for Na-ion and Li-ion batteries with using a proper substrate: a first principles study, Comp. Mater. Sci. (2016) 166–172doi:10.1016/j.commatsci.2016.07.034.
  * [33] X. Zhang, Z. Yu, S. S. Wang, S. Guan, H. Y. Yang, Y. Yao, S. A. Yang, Theoretical prediction of MoN2 monolayer as a high capacity electrode material for metal ion batteries, J. Mater. Chem. A 4 (2016) 15224–15231.  doi:10.1039/C6TA07065E.
  * [34] M. Dahbi, N. Yabuuchi, K. Kubota, K. Tokiwa, S. Komaba, Negative electrodes for Na-ion batteries, Phys. Chem. Chem. Phys. 16 (2014) 15007–15028.  doi:10.1039/C4CP00826J.
  * [35] D. Kundu, E. Talaie, V. Duffort, L. F. Nazar, The emerging chemistry of sodium ion batteries for electrochemical energy storage, Angew. Chem. Int. Ed. Engl. 46 (2015) 3431–3448.  doi:10.1002/anie.201410376.
  * [36] P. J. P. Naeyaert, M. Avdeev, N. Sharma, H. B. Yahia, C. D. Ling, Synthetic, structural, and electrochemical study of monoclinic Na4Ti5O12 as a sodium-ion battery anode material, Chem. Mater. 26 (2014) 7067–7072.  doi:10.1021/cm5035358.
  * [37] P. Senguttuvan, G. Rousse, V. Seznec, J. M. Tarascon, M. R. Palacin, Na2Ti3O7: lowest voltage ever reported oxide insertion electrode for sodium ion batteries, Chem. Mater. 23 (2011) 4109–4011.  doi:10.1021/cm202076g.
  * [38] S. Komaba, Y. Matsuura, T. Ishikawa, N. Yabuuchi, W. Murata, S. Kuze, Redox reaction of Sn-polyacrylate electrodes in aprotic Na cell, Electrochem. Commun. 21 (2012) 65–68.  doi:10.1016/j.elecom.2012.05.017.
  * [39] J. Qian, Y. Chen, L. Wu, Y. Cao, X. Ai, H. Yang, High capacity Na-storage and superior cyclability of nanocomposite Sb/C anode for Na-ion batteries, Chem. Commun. 48 (2012) 7070–7072.  doi:10.1039/C2CC32730A.
  * [40] D. B. Putungan, S. H. Lin, J. L. Kuo, Metallic VS2 monolayer polytypes as potential sodium-ion battery anode via ab initio random structure searching, ACS Appl. Mater. Interfaces 8 (2016) 18754–18762.  doi:10.1021/acsami.6b03499.
  * [41] F. Li, Y. Qu, M. Zhao, Germanium sulfide nanosheet: a universal anode material for alkali metal ion batteries, J. Mater. Chem. A 4 (2016) 8905–8912.  doi:10.1039/C6TA03210A.

Figure 1: (a) Top view for the adsorption of metal atom at different sites on
the MoO2 sheet. (b) Side view of adsorbate on the H site. (c) The adsorption
energy for Na, K, Ca based on the four distinct binding sites on MoO2. Figure
2: Total and projected density of states for (a) Na@MoO2, (b) K@MoO2, and (c)
Ca@MoO2, respectively. The horizontal dashed lines represent the Fermi level.
Figure 3: Charge density difference of the (110) section for single adatom (a)
Na, (b) K and (c) Ca (being adsorbed on the H site) on the MoO2. Red and blue
colors indicate the electrons accumulation and depletion, respectively. The
color scale is in the units of e/Å3.  Figure 4: (a) the top and side views of
the optimized migration pathways and the corresponding diffusion barrier
profiles of (b) Na (black), K (red), Ca (blue) on the MoO2 sheet,
respectively.  Figure 5: Top and side views of the optimized structures of (a)
Na0.25MoO2, (b) Na1MoO2, and (c) Na3MoO2. The color of yellow, green and
violet represent the first, second and third Na layer, respectively. Figure 6:
The variation of OCVs with the increasing concentration of metal atoms on the
MoO2 surface. The dashed open circles correspond to the maximum $x$ in the
compounds