# Thermophysical Properties of Molten FLiNaK: a Moment Tensor Potential
Approach

Nikita Rybin  Dmitrii Maksimov  Yuriy Zaikov  Alexander Shapeev

###### Abstract

Fluoride salts demonstrate significant potential for applications in next-
generation nuclear reactors, necessitating a comprehensive understanding of
their thermophysical properties for technological advancements. Experimental
measurement of these properties poses challenges, due to factors such as high
temperatures, impurity control, and corrosion. Consequently, precise
computational modeling methods become essential for predicting the
thermophysical properties of molten salts. In this work, we performed
molecular dynamics (MD) simulations of several thermophysical properties of
the eutectic salt mixture LiF-NaF-KF (FLiNaK) melt, including density, self-
diffusion coefficients, viscosity, and thermal conductivity. We demonstrated
the successful application of moment tensor potentials (MTP) as an accurate
model for interatomic interactions in FLiNaK. Our results on thermophysical
properties calculations exhibit strong agreement with experimental data. An
important aspect of our methodology is the incorporation of an active learning
scheme, which enables the generation of a robust and accurate potential, while
maintaining a moderate-sized training dataset.

###### keywords:

FLiNaK , molten salts , molten salt density , molten salt diffusivity , molten
salt thermal conductivity , moment tensor potential , active learning

††journal: Journal of Molecular Liquids

[first]organization=Skolkovo Institute of Science and
Technology,addressline=Bolshoy Boulevard 30, bld. 1, city=Moscow,
postcode=121205, country=Russia

[second]organization=Institute of High Temperature Electrochemistry, Ural
Branch of the Russian Academy of Sciences,addressline=Academicheskaya Str. 20,
city=Yekaterinburg, postcode=620066, country=Russia

##  1 Introduction

Molten salts play a pivotal role in various industrial applications, including
molten salt batteries Ong et al. (2020); Bell et al. (2019); Cui et al.
(2022); Parasotchenko et al. (2023); LeBlanc (2010), molten salt reactor
systems Beneš and Konings (2009); Magnusson et al. (2020), as well as in
pyroprocessing methods for recycling nuclear fuel Zhitkov et al. (2020, 2022);
Mullabaev et al. (2022). Most important characteristics of melts from an
application standpoint are temperature-dependent values of density,
diffusivity, viscosity and thermal conductivity, which are usually measured
experimentally. Despite the availability of numerous experimental measurements
for different molten salts, reported thermophysical properties exhibit
significant scatter and uncertainties, reaching up to 20% Romatoski and Hu
(2017). This variability is attributed to challenges associated with
impurities, high-temperature measurements, and deviations in composition.

As an alternative to experimental methods, accurate theoretical approaches
based on ab initio molecular dynamics (AIMD) simulations have emerged as
valuable tools for evaluating thermophysical properties in well-controlled
conditions Porter et al. (2022). While AIMD has significantly advanced our
understanding of the structure of various molten salts, essential properties
such as diffusion coefficients, viscosity, and thermal conductivity
necessitate larger unit cell sizes and extended simulation times for relevant
statistical analysis Gheribi et al. (2014); Jabbari et al. (2017). Thus, there
is a critical need for an efficient computational method to reliably predict
the thermophysical properties of molten salts.

An alternative computational approach involves MD simulations, where
interatomic interaction potential is fitted to first-principles calculations
or experimental data Salanne et al. (2009); Galashev et al. (2023); Maxwell
(2022). This overcomes limitations related to simulation cell sizes and time
scales, enabling simulations of larger structures for more extended periods
compared to AIMD. Properties derived from classical MD simulations generally
achieve sufficient numerical precision, with errors dominated by the accuracy
of the underlying interatomic potentials Lu et al. (2021); Lee et al. (2021).

In the last decade, significant progress has been made in developing the so-
called machine learning interatomic potentials (MLIPs) Deringer et al. (2019).
MLIPs have demonstrated promise for MD simulations of molten salts modelling
with near ab initio accuracy, on time and length scales comparable to
traditional interatomic potentials Liang et al. (2020); Feng et al. (2022);
Lam et al. (2021); Li et al. (2021); Rodriguez et al. (2021); Attarian et al.
(2022). The moment tensor potential Shapeev (2016) (MTP) was shown to be among
the most efficient (in terms of data utilization) and accurate machine-learned
models for interatomic interactions Zuo et al. (2020). In the case of molten
LiF-BeF2 mixture, it was recently shown that MTP can be used to compute
thermophysical properties with both high precision and low data utilization
Attarian et al. (2022).

In this study, we employ MTPs to calculate thermophysical properties of molten
FLiNaK at eutectic composition (46.5-11.5-42 mol%). FLiNaK is among the most
promising and well-investigated salts Locatelli et al. (2013); Maslennikova et
al. (2023); Lizin et al. (2017), with a lot of experimental and theoretical
data. We fitted an MTP to approximate the potential energy surface of FLiNaK
on the data calculated from first principles in the dispersion-corrected
Density-Functional Theory (DFT-D3) framework. As will be demonstrated, MTP
model used in couple with MD simulations allows to obtain thermophysical
properties of FLiNaK in a good agreement with reported literature data, both
experimental and theoretical.

##  2 Methodology

###  2.1 Moment Tensor Potential

In this work, we used MTP approach implemented in the MLIP-2 package Novikov
et al. (2020) to investigate the thermophysical properties of molten FLiNaK at
finite temperatures. The potential energy of an atomic system as described by
the MTP interatomic potential is defined as a sum of the energies of atomic
environments of the individual atoms:

$E_{\text{MTP}}=\sum_{i=1}^{N}V(n_{i}),$

where the index $i$ label $N$ atoms of the system, and $n_{i}$ describes the
local atomic neighborhood around atom i within a certain cutoff radius
$R_{\text{cut}}$ and the function $V$ is the moment tensor potential:

$V(n_{i})=\sum_{\alpha}\xi_{\alpha}B_{\alpha}(n_{i}),$

where $\xi_{\alpha}$ are the fitting parameters and $B_{\alpha}(n_{i})$ are
the basis functions that will be defined below. Moment tensors descriptors are
used as representations of atomic environments and defined as:

$M_{\mu,\nu}\left({n}_{i}\right)=\sum_{j}f_{\mu}\left(\left|r_{ij}\right|,z_{i},z_{j}\right)\underbrace{r_{ij}\otimes\ldots\otimes
r_{ij}}_{\nu\text{ times }},$

where the index $j$ goes through all the neighbors of atom $i$. The symbol
“$\otimes$” stands for the outer product of vectors, thus
${r}_{ij}\otimes\cdots\otimes{r}_{ij}$ is the tensor of rank $\nu$ encoding
the angular part which itself resembles moments of inertia. The function
$f_{\mu}$ represents the radial component of the moment tensor:

$f_{\mu}\left(\left|r_{ij}\right|,z_{i},z_{j}\right)=\sum_{k}c_{\mu,z_{i},z_{j}}^{(k)}Q^{(k)}(r),$

where $z_{i}$ and $z_{j}$ denote the atomic species of atoms $i$ and $j$,
respectively, $r_{ij}$ describes the positioning of atom $j$ relative to atom
$i$, $c_{\mu,z_{i},z_{j}}^{(k)}$ are the fitting parameters and

$Q^{(k)}(r):=T_{k}(r)\left(R_{\text{cut }}-r\right)^{2}$

are the radial functions consisting of the Chebyshev polynomials $T_{k}(r)$ on
the interval $[R_{\text{min}},R_{\text{cut}}]$ with the term
($R_{\text{cut}}-r)^{2}$ that is introduced to ensure a smooth cut-off to
zero. The descriptors $M_{\mu,\nu}$ taking $\nu$ equal to $0,1,2,\ldots$ are
tensors of different ranks that allow to define basis functions as all
possible contractions of these tensors to a scalar, for instance:

$\displaystyle B_{0}\left({n}_{i}\right)=M_{0,0}\left({n}_{i}\right),$
$\displaystyle B_{1}\left({n}_{i}\right)=M_{0,1}\left({n}_{i}\right)\cdot
M_{0,1}\left({n}_{i}\right),$ $\displaystyle
B_{2}\left({n}_{i}\right)=M_{0,0}\left({n}_{i}\right)\left(M_{0,2}\left({n}_{i}\right):M_{0,2}\left({n}_{i}\right)\right).$

Therefore the level of $M_{\mu,\nu}$ is defined by ${\rm lev}M_{\mu,\nu}$ =
2$\mu$ \+ $\nu$ and if $B_{\alpha}$ is obtained from $M_{\mu_{1},\nu_{1}}$ ,
$M_{\mu_{2},\nu_{2}}$ , $\dots$, then ${\rm lev}B_{\alpha}$ =
$(2\mu_{1}+\nu_{1})$ \+ $(2\mu_{2}+\nu_{2})$ \+ $\dots$ . By including all
basis functions such that ${\rm lev}B_{\alpha}<d$ we obtain the moment tensor
potential of level $d$, which we denote as MTPd.

###  2.2 Dataset Generation and Potential Fitting

The MTP utilized in this study is trained using data computed within the
Density-Functional Theory (DFT) framework. All DFT calculations were carried
out using VASP (Vienna ab initio simulations package) Kresse and Furthmüller
(1996) with the projector augmented wave method Kresse and Joubert (1999). The
Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) Perdew et
al. (1996) was employed for the exchange–correlation functional, and the
DFT-D3 method Grimme et al. (2010) was utilized to account for dispersion
forces.

The initial dataset was generated from four independent AIMD trajectories (at
temperatures $T=800$ K, $1000$ K, $1200$ K, and $1400$ K), which were
conducted in the isothermal-isobaric ensemble (NPT) using the Nosé-Hoover
thermostat Nosé (1984). Each AIMD trajectory was simulated for 2 ps with a 1
fs time step and for supercells containing 56 atoms (28 F atoms, 13 Li atoms,
3 Na atoms, and 12 K atoms) to closely mimic eutectic composition. The plane-
wave basis set had an energy cutoff of 550 eV, and a single gamma point was
used to sample the Brillouin zone. The initial structures for AIMD were
generated by randomly distributing atoms in the cubic cell, which had the
volume equal to the total volume of spheres, each with a radius matching the
van der Waals radii of the respective atoms. The first picosecond of each AIMD
simulation was discarded, and the second picosecond was subsampled with the 5
fs time intervals, resulting in a set of 800 samples. These samples were
employed to train an initial MTP of level 14.

The root-mean-squared-errors (RMSEs) on this dataset are 2.1 meV/atom for
energies and 48 meV/$\AA$ for forces. Previous studies Feng et al. (2022); Lam
et al. (2021) have indicated that RMSE of energy smaller than 5 meV/atom and
RMSEs on forces smaller than 100 meV/$\AA$ are generally sufficient for
predicting properties such as density, radial distribution functions,
diffusion coefficient, and viscosity of molten salts. However, achieving this
accuracy does not by itself ensure the robustness (ability to run long MD) of
the potential. Generally, two strategies are employed to train a robust MLIP
capable of capturing the entire range of local atomic environments: (i)
including a variety of different systems which requires lengthy AIMD
simulations, and (ii) utilizing an active learning (active sampling, learning
on-the-fly) strategy to selectively add data points on which the MLIP
extrapolates significantly, i.e., the prediction of energies and forces of
these structures is done with high uncertainty.

Table 1: Potential energy and force RMSEs on training and test sets for the fitted Moment Tensor Potential. | Train | Test |  |
---|---|---|---|---
Energy (meV/atom) | 1.83 | 1.93 |  |
Forces (meV/$\AA$) | 39 | 39 |  |

For the purpose of methodological validation, we have explored both
strategies. We extended AIMD simulations by an additional 3 ps at each
temperature, yielding (after subsampling) a dataset of 3200 samples. We
shuffled this dataset and divided it into the training (20% of the total
dataset) and test sets (remaining 80% of the dataset). Fig. 1 presents a
comparative analysis of DFT and MTP calculated energies and forces for both
training and test sets. As presented in Tab. 1, the RMSE of MTP on energies is
just 1.93 meV/atom, and on forces, it is 39 meV/$\AA$. The displayed RMSE
values in Fig. 1 demonstrate that the MTP, trained on a 20% subset of
configurations derived from AIMD simulations, accurately predicts the energies
and forces of the remaining 80% (test set) with comparable errors to those of
the training set. Low errors on the test set, suggests that the MTP is well-
fitted and likely capable of extrapolating to diverse structures and chemical
environments.

However, even with these promising results, an ability to run long MD
simulations without failure remains to be demonstrated. To assess stability of
the fitted potential, we conducted MD simulations using MTP as the model for
interatomic interactions in the micro-canonical (NVE) ensemble within LAMMPS
(Large-scale Atomic/Molecular Massively Parallel Simulator) Thompson et al.
(2022). Unfortunately, these simulations proved to be unstable, exhibiting
instability after several tens of picoseconds. Notably, extending the
potential training to the entire dataset did not alleviate this issue. This
outcome underscores that while low energy and force errors are achieved, they
do not guarantee the robust applicability of a MLIP for practical MD
simulations, even under the same thermodynamic conditions as those employed in
the training data.

We now revisit the active learning strategy, which means interactively
selecting a diverse, but minimal, set of training data in the feature space to
effectively fit the potential. Various active learning schemes are employed in
the development of machine learning (ML) potentials Zhang et al. (2019);
Sivaraman et al. (2020). In our study, we utilized the D-optimality-based
active learning procedure developed in Podryabinkin and Shapeev (2017) and
available in the MLIP-2 package.

Within the active learning algorithm, we initiate MD calculations in LAMMPS
using a pretrained MTP. At each step of the MD simulations, the algorithm
assesses the extrapolation grade $\gamma$ of the atomic configuration based
solely on atomic coordinates. Configurations with $\gamma>2$ are added to the
preselected set. When $\gamma$ exceeds 10, the MD simulation halts, and all
sufficiently different configurations from the preselected set are
incorporated into the training, followed by the refitting of the potential.
This procedure repeats until MD simulations can run without failure for 200
ps. In our case, MTP achieved robustness with a training set comprising 880
samples, indicating that only a small portion of the actively selected
configurations required single-point DFT calculations. The potential trained
in this manner demonstrates the ability to robustly conduct MD simulations for
at least hundreds of picoseconds, showcasing the effectiveness of the active
learning procedure in terms of reducing the necessary DFT data for training
and enhancing potential robustness. The training set generated during the
active learning procedure can later be used to train a MTP with enhanced
accuracy, as will be demonstrated in this study.

Thus, the active learning scheme plays a crucial role in enhancing the
robustness of any MLIP, and MTP, in particular. Its employment not only
ensures the applicability of MTP at large length scales and extended time
scales in the MTP-MD simulations, but also significantly reduces the number of
required DFT calculations. In the subsequent sections of the paper, we present
the results of our calculations on the thermophysical properties of FLiNaK.
Specifically, we assess the MTP’s capability to accurately calculate the
temperature dependencies in density, diffusivity, viscosity, and thermal
conductivity.

(a)

(b)

Figure 1: Parity plots of the energies (a) and forces (b) for training (blue)
and test (orange) sets. The diagonal line in each figure shows the perfect
fit.

##  3 Results

###  3.1 Radial Distribution Functions and Self-Diffusion Coefficients

To assess how well the developed MTP can predict local structural features of
molten FLiNaK, we first calculated radial distribution functions (RDF) and
diffusion coefficients of ions. The smallest unit cell of molten FLiNaK in an
eutectic composition contains 400 atoms: 93 Li atoms, 200 F atoms, 23 Na
atoms, and 84 K atoms, respectively. For MTP-MD simulations, the simulation
cell was replicated in a 2$\times$2$\times$2 arrangement, which resulted in a
supercell with 3200 atoms. The usage of larger supercells yields better
statistics of calculations and MTP-MD allows to perform long simulations,
which is inaccessible by AIMD. The first RDF peaks for F-Li, F-Na, and F-K are
shown in Fig. 2. RDF for Li-F shows a sharp peak at 1.842 $\AA$ and then
rapidly decays. This suggests a strong bonding between Li and F within FLiNaK.
On the other hand, the RDF for K-F has a wider first peak at 2.595 $\AA$ that
does not decay to zero. This result suggests that the first nearest neighbor
shell for K and F is more diffusive and not as well defined as for Li and F,
which further indicated that the bonding between K and F in FLiNaK is not as
strong as between Li and F. The values of the calculated first peak distances
obtained using MTP-MD for F-Li, F-Na, and F-K are close to the experimental
values Igarashi et al. (1988) as shown in an inset of Fig. 2. This indicates
that MTP can accurately capture the local environment of ions in molten
FLiNaK.

Figure 2: Radial distribution functions, extracted from MTP-MD at temperature
T=793 K for F-Li, F-Na, and F-K ion pairs. Experimental nearest-neighbor peak
distance for aforementioned ion pairs was determined from X-ray scattering
experiment Igarashi et al. (1988).

Self-diffusion of atoms in a melt is another fundamental property of liquid
dynamics that provides important structural information Rollet et al. (2009).
Fig. 3(a, b, c, d) show the temperature-dependence of self-diffusion
coefficients of F, Li, Na, and K on a logarithmic scale, obtained from MTP-MD
and experiment Umesaki et al. (1981). Although the experimental data is
limited to a relatively short temperature intervals, for all four atomic
species the calculated values and the slope of temperature dependencies are
both close to the experimental results. Such results further validates our
approach and demonstrated benefits of using atomistic simulations to accompany
the experiment.

(a)

(b)

(c)

(d)

Figure 3: Temperature dependence of the self-diffusion coefficients for (a)
Li, (b) K, (c) Na, (d) F. Experimental data is taken from Umesaki et al.
(1981). Error bars are calculated from the standard deviation among 5
separately trained MTPs.

###  3.2 Density

We next study the temperature dependence of the density as one of the most
technologically important thermophysical characteristics of the molten salt.
While experimental measurements are feasible, computational modeling enables
fast and inexpensive exploration of various salt mixtures with diverse
components and compositions. The density of molten FLiNaK has been
experimentally measured at different temperatures Romatoski and Hu (2017); An
et al. (2017); W. D. Powers and Greene (1963), and recent efforts have
employed ML potentials based on neural networks to calculate it Lee et al.
(2021).

We perform density calculations of molten FLiNaK, using the same simulation
cell as for the diffusion calculations (comprising 3200 atoms at eutectic
composition). First, we investigated the influence of different levels of MTP
(i.e., varying numbers of parameters) in predicting the density of melt at
T=1000K. Fig. 4 illustrates that beyond the level 16, increasing the
complexity of MTP does not substantially enhance density values towards
experimental results, while computational costs grow. For all subsequent
calculations in this study, we employed the potential with level 16. Notably,
such a clear convergence of the potential complexity and the property of
interest (density in this case) as in Fig. 4 is, however, not always observed.
We also note that fitting of the MTPs is based on the Broyden-Fletcher-
Goldfarb-Shanno method, and the optimized parameters depended on the
initialization of MTP parameters, resulting in varying values and errors for
the targeted property. To address this variability, we conducted 5
optimization sessions for each potential level, with error bars depicted in
Fig. 4. Our findings indicate that the level of as low as 10 can yield
reasonably accurate results for molten FLiNaK density. The errors on energies
and forces exhibit nearly identical patterns across different potential levels
(from 10 to 20), aligning with observations in a recent study of molten LiF-
BeF2, where MTP was also used Attarian et al. (2022).

Next, we proceed toward the comparison of our results with previous
theoretical calculations. Recently, DeepMD potential was used to study molten
FLiNaK Lee et al. (2021). We performed density calculations for temperatures
T=600 K, 800 K, 1000 K, 1200 K. Fig. 5 reveals an overall agreement between
our results and those reported in Lee et al. (2021). Notably, the
configurations used for potential training in Lee et al. (2021) were sampled
at a fixed experimental density in the NVT ensemble, whereas we performed AIMD
simulations with the NPT ensemble.

Given that FLiNaK is an ionic compound, accounting for long-range interactions
becomes crucial in atomic dynamics calculations. Although the functional form
of MTP generated by the MLIP-2 package currently lacks explicit terms for
dispersion corrections, our tests of the convergence of energies and forces
with respect to the cutoff radius of interactions showed no substantial
improvements beyond a cutoff radius of 5.5 $\AA$. This distance is deemed
sufficient to include multiple neighbor shells in the liquid, with negligible
electrostatic interactions expected beyond this distance for near-equilibrium
FLiNaK configurations. While directly incorporating long-range interactions
treatment into potential construction could potentially enhance simulation
accuracy, such functionality has not been developed yet.

To the best of our knowledge, the literature lacks an in-depth exploration of
the impact of long-range interactions on thermophysical properties
calculations of molten salts. A study Frandsen et al. (2020) demonstrated that
including vdW-corrections increased the error in density determination by 4%.
In our case, Fig. 5 illustrates how a potential trained on data with the
DFT-D3 correction leads to slightly different results (about a 10% difference)
compared to a potential trained without the the D3 correction. In our case,
including long-range interactions during the dataset generation step improved
agreement with experimental results. This test was performed by generating a
dataset using the active learning methodology described earlier.

Finally, Fig. 6 presents a comparison of temperature-dependent density of
molten FliNaK obtained using MTP-MD with several experimental results An et
al. (2017); W. D. Powers and Greene (1963); Romatoski and Hu (2017). Our
findings align well with the experimental results, despite not fixing
experimentally-known density during the preparation of the training set. Given
that experimental density may not be available for other salts, we emphasize
this as an appropriate approach for conducting every MLIP benchmarks.

Figure 4: Density of FLiNaK calculated at T=1000K using MTP with different
levels of complexity (number of parameters). The experimental data is taken
from An et al. (2017); W. D. Powers and Greene (1963); Romatoski and Hu
(2017). Error bars represent one standard deviation calculated over 5
separately trained MTPs. Figure 5: Temperature dependence of the density of
FLiNaK calculated using MLIP. Olive disks represent results of MTP calculated
on a dataset in which D3 dispersion corrections was not taken into account,
while red disks represent MTP, which was trained on dataset with the D3
correction Grimme et al. (2010). Density calculated using DeepMD (orange
triangles) is taken from Lee et al. (2021). Error bars represent one standard
deviation calculated over 5 separately trained MTPs. Figure 6: Temperature
dependence of the density of FLiNaK. The experimental data is taken from An et
al. (2017); W. D. Powers and Greene (1963); Romatoski and Hu (2017). Error
bars represent one standard deviation calculated over 5 separately trained
MTPs.

In the following subsections we will demonstrate how MTP can be used in
conjunction with the Green-Kubo method Kubo et al. (1957); Green (2004) to
calculate viscosity of molten salt and with the Müller-Plathe non-equilibrium
method Müller-Plathe (1997) for thermal conductivity calculations. It should
be noted that due to the high computational cost of conducting viscosity and
thermal conductivity calculations the reported values are presented without
error bars.

###  3.3 Viscosity

We calculated viscosity using the Green-Kubo (GK) approach Kubo et al. (1957);
Green (2004), employing the implementation available in the LAMMPS package.
The GK approach involves the calculation of viscosity through the integral of
the auto-correlation function of the off-diagonal elements of the stress
tensor, expressed by the following relation:

$\eta=\frac{V}{k_{B}T}\int_{0}^{\infty}\braket{P_{\alpha\beta}(t)P_{\alpha\beta}(0)}dt,$

where $\eta$ is viscosity, $k_{B}$ is the Boltzmann constant, and
$P_{\alpha\beta}$ are the off-diagonal elements of the stress tensor.

In our viscosity calculations we again utilized a simulation cell containing
3200 atoms, with a simulation time step of 1 fs. For all temperatures, an
auto-correlation time of 10 ps was selected ensuring the convergence of the
auto-correlation function of the diagonal stress components to zero. After the
initial equilibration at each temperature, the simulation was extended for 5
ns, employing the NVE ensemble.

Fig. 7 illustrates the temperature-dependent viscosity of FLiNaK calculated
here and experimentally J. Ambrosek and Allen (2009); Vriesema (1979); Rudenko
et al. (2022a). Our results exhibit excellent agreement with the experimental
data. The relationship between viscosity and temperature established in our
work is described by the following equation:

$\eta=0.00543\cdot\exp\left(\frac{5422}{T}\right)~{}[{\rm mPa}\cdot{\rm s}].$

Figure 7: Temperature dependence of the viscosity of molten FLiNaK (eutectic
composition). Experimental data is taken from J. Ambrosek and Allen (2009);
Vriesema (1979); Rudenko et al. (2022a). Figure 8: Temperature dependence of
the thermal conductivity of molten FLiNaK (eutectic composition). The
experimental results are taken from Smirnov et al. (1987); Sohal et al.
(2010); Gallagher et al. (2022); Robertson et al. (2022); Rudenko et al.
(2022b). In Smirnov et al. (1987) $\kappa$ was measured with a correction of
$\pm$0.012 W/mK, introducing an uncertainty of 1-5%.

###  3.4 Thermal Conductivity

Before delving into the discussion of our heat transport calculations, it is
essential to provide an overview of the available experimental data on the
temperature dependence of a thermal conductivity (labeled as $\kappa$ from now
on). The earliest works provide $\kappa$ values ranging from 1.3 W/mK to 4.5
W/mK, Vriesema (1979); Grele and Gedeon (1953); Hoffman and Lones (1955). Some
experimental works Smirnov et al. (1987); Sohal et al. (2010) indicate a
slight increase in thermal conductivity with temperature, as depicted in Fig.
8. Work Williams (2006) indicated that the thermal conductivity at 973 K will
be in the range 0.6–1.0 W/mK. Most recent experiments Gallagher et al. (2022);
Robertson et al. (2022); Rudenko et al. (2022b) agree in suggesting a slight
decrease of $\kappa$ with temperature, indicating that the thermal
conductivity value should be below 1 W/mK.

We determined thermal conductivity of molten FLiNaK using the Müller-Plathe
non-equilibrium method Müller-Plathe (1997). It was demonstrated in Pan et al.
(2021) that this approach allows computing heat transport in molten salts in a
better agreement with experimental values, comparably to the Green-Kubo
method. We used a time step of 0.5 fs and the supercell with 4000 atoms (2.6
nm $\times$ 2.9 nm $\times$ 23 nm dimensionality) and a kinetic energy swap
rate of 1 in every 1000 steps. For each temperature, after the initial
equilibration in the NPT ensemble for 10 ps, the simulations were done for 4
ns in the NVE ensemble. Fig. 8 shows results of our calculations and mentioned
above experimental data. The linear fit of our results lead to the following
equation for the thermal conductivity:

$\kappa=0.880-0.116\cdot 10^{-3}T\text{ [W/m}\cdot\text{K]}$

As demonstrated in Fig. 8, although the experimental data is slightly
scattered, our calculated values agree well with the most recent experimental
works Gallagher et al. (2022); Robertson et al. (2022); Rudenko et al.
(2022b). Taken into account the magnitude of thermal conductivity changes with
temperature, we are of the opinion that it is reasonable to state that the
thermal conductivity in FLiNaK stays nearly constant with temperature grows
and its value is below 1 W/mK.

##  4 Summary and conclusions

In this work, we utilized a MTP and assessed its performance for the
calculations of thermophysical properties of molten FLiNaK in a range of
temperatures. Employing an active learning approach enabled the rapid training
of a robust potential capable of predicting density, diffusion, viscosity, and
thermal conductivity with near ab initio accuracy. While previous
computational works demonstrated good agreement with experimental data, the
MTP, coupled with active learning, stands out for its ability to deliver
results much faster, thanks to its efficient utilization of data. This
significantly reduces the computational cost compared to previous studies
employing alternative MLIPs. This breakthrough paves the way for swift yet
precise exploration of thermophysical properties across various salt systems.

Our analysis, encompassing radial distribution function and diffusivity
calculations, demonstrates the MTP’s accuracy in predicting local structures
of molten FLiNaK. The extrapolation of the potential to larger systems,
achieved using a smaller simulation cell during the ML potential training,
underscores its versatility in handling diverse local configurations. While
the MTP has a medium-range cutoff distance (up to 5.5 $\AA$ in this study) and
does not explicitly consider long-range interactions, our results indicate
that it performs well in ionic systems. The absence of explicit charge
treatment in the MTP, justified by local ionic solvation shells limiting long-
range interactions, makes it a useful tool for fast and accurate computations
of molten salt’s thermophysical properties. However, certain conditions must
be met for optimal results. Our findings suggest that accounting for van der
Waals dispersion during the generation of training data improves predictions
of the temperature dependence of density. While the potential currently lacks
an implicit treatment of long-range interactions, recent research Ko et al.
(2021); Gao and Remsing (2022) highlights the possible benefits of their
explicit inclusion in charge systems. Such features should be included in the
next generation of MTPs to improve modeling accuracy.

We also note that utilization of different DFT approximations might improve
the results for molten salts calculations with MLIPs. In particular, it was
recently shown in Tisi et al. (2021) that training of a MLIP on the dataset
generated at the DFT-PBE level of accuracy leads to the predicted values of
the thermal conductivity of water 60% larger than in experiments. Retraining
on the dataset generated with the strongly constrained and appropriately
normed semilocal density functional (SCAN) Sun et al. (2015, 2016) decreases
errors by two times. Proper benchmarking of the exchange-correlation density
functional influence on the evaluated properties of molten salts is, to the
best of our knowledge, currently missing in literature and have to be explored
in future works.

## Acknowledgements

N.R. and A.S. acknowledge funding from the Russian Science Foundation (Project
No. 23-13-00332).

## References

  * An et al. (2017) An, X.H., Cheng, J.H., Su, T., Zhang, P., 2017.  Determination of thermal physical properties of alkali fluoride/carbonate eutectic molten salt.  AIP Conference Proceedings 1850, 070001.  URL: https://doi.org/10.1063/1.4984415, doi:10.1063/1.4984415.
  * Attarian et al. (2022) Attarian, S., Morgan, D., Szlufarska, I., 2022.  Thermophysical properties of flibe using moment tensor potentials.  Journal of Molecular Liquids 368, 120803.  URL: https://www.sciencedirect.com/science/article/pii/S016773222202342X, doi:https://doi.org/10.1016/j.molliq.2022.120803.
  * Bell et al. (2019) Bell, S., Steinberg, T., Will, G., 2019.  Corrosion mechanisms in molten salt thermal energy storage for concentrating solar power.  Renewable and Sustainable Energy Reviews 114, 109328.  URL: https://www.sciencedirect.com/science/article/pii/S1364032119305362, doi:https://doi.org/10.1016/j.rser.2019.109328.
  * Beneš and Konings (2009) Beneš, O., Konings, R., 2009.  Thermodynamic properties and phase diagrams of fluoride salts for nuclear applications.  Journal of Fluorine Chemistry 130, 22–29.  URL: https://www.sciencedirect.com/science/article/pii/S0022113908002030, doi:https://doi.org/10.1016/j.jfluchem.2008.07.014. fluorine and Nuclear Energy.
  * Cui et al. (2022) Cui, K., Zhao, W., Li, S., Zhou, D., Liu, C., Qu, X., Li, P., 2022.  Low-temperature and high-energy-density li-based liquid metal batteries based on licl–kcl molten salt electrolyte.  ACS Sustainable Chemistry & Engineering 10, 1871–1879.  URL: https://doi.org/10.1021/acssuschemeng.1c07560, doi:10.1021/acssuschemeng.1c07560.
  * Deringer et al. (2019) Deringer, V.L., Caro, M.A., Csányi, G., 2019.  Machine learning interatomic potentials as emerging tools for materials science.  Advanced Materials 31, 1902765.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.201902765, doi:https://doi.org/10.1002/adma.201902765.
  * Feng et al. (2022) Feng, T., Zhao, J., Liang, W., Lu, G., 2022.  Molecular dynamics simulations of lanthanum chloride by deep learning potential.  Computational Materials Science 210, 111014.  URL: https://www.sciencedirect.com/science/article/pii/S0927025621007011, doi:https://doi.org/10.1016/j.commatsci.2021.111014.
  * Frandsen et al. (2020) Frandsen, B.A., Nickerson, S.D., Clark, A.D., Solano, A., Baral, R., Williams, J., Neuefeind, J., Memmott, M., 2020.  The structure of molten flinak.  Journal of Nuclear Materials 537, 152219.  URL: https://www.sciencedirect.com/science/article/pii/S0022311520304736, doi:https://doi.org/10.1016/j.jnucmat.2020.152219.
  * Galashev et al. (2023) Galashev, A.Y., Rakhmanova, O.R., Abramova, K.A., Katin, K.P., Maslov, M.M., Tkacheva, O.Y., Rudenko, A.V., Kataev, A.A., Zaikov, Y.P., 2023.  Molecular dynamics and experimental study of the effect of cef3 and ndf3 additives on the physical properties of flinak.  The Journal of Physical Chemistry B 127, 1197–1208.  URL: https://doi.org/10.1021/acs.jpcb.2c06915, doi:10.1021/acs.jpcb.2c06915. pMID: 36696698.
  * Gallagher et al. (2022) Gallagher, R.C., Birri, A., Russell, N.G., Phan, A.T., Gheribi, A.E., 2022.  Investigation of the thermal conductivity of molten lif-naf-kf with experiments, theory, and equilibrium molecular dynamics.  Journal of Molecular Liquids 361, 119151.  URL: https://www.sciencedirect.com/science/article/pii/S0167732222006894, doi:https://doi.org/10.1016/j.molliq.2022.119151.
  * Gao and Remsing (2022) Gao, A., Remsing, R.C., 2022.  Self-consistent determination of long-range electrostatics in neural network potentials.  Nature Communications 13, 1572.  URL: https://doi.org/10.1038/s41467-022-29243-2, doi:10.1038/s41467-022-29243-2.
  * Gheribi et al. (2014) Gheribi, A.E., Corradini, D., Dewan, L., Chartrand, P., Simon, C., Madden, P.A., Salanne, M., 2014.  Prediction of the thermophysical properties of molten salt fast reactor fuel from first-principles.  Molecular Physics 112, 1305–1312.  URL: https://doi.org/10.1080/00268976.2014.897396, doi:10.1080/00268976.2014.897396.
  * Green (2004) Green, M.S., 2004.  Markoff Random Processes and the Statistical Mechanics of Time‐Dependent Phenomena. II. Irreversible Processes in Fluids.  The Journal of Chemical Physics 22, 398–413.  URL: https://doi.org/10.1063/1.1740082, doi:10.1063/1.1740082.
  * Grele and Gedeon (1953) Grele, M.D., Gedeon, L., 1953.  Forced-convection heat-transfer characteristics of molten sodium hydroxide.
  * Grimme et al. (2010) Grimme, S., Antony, J., Ehrlich, S., Krieg, H., 2010.  A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu.  The Journal of Chemical Physics 132, 154104.  URL: https://doi.org/10.1063/1.3382344, doi:10.1063/1.3382344.
  * Hoffman and Lones (1955) Hoffman, H., Lones, J., 1955.  Fused salt heat transfer. part ii. forced convection heat transfer in circular tubes containing naf-kf-lif eutectic.
  * Igarashi et al. (1988) Igarashi, K., Okamoto, Y., Mochinaga, J., Ohno, H., 1988.  X-ray diffraction study of molten eutectic lif–naf–kf mixture.  J. Chem. Soc., Faraday Trans. 1 84, 4407–4415.  URL: http://dx.doi.org/10.1039/F19888404407, doi:10.1039/F19888404407.
  * J. Ambrosek and Allen (2009) J. Ambrosek, M. Anderson, K.S., Allen, T., 2009.  Current status of knowledge of the fluoride salt (flinak) heat transfer.  Nuclear Technology 165, 166–173.  URL: https://doi.org/10.13182/NT165-166, doi:10.13182/NT165-166.
  * Jabbari et al. (2017) Jabbari, F., Rajabpour, A., Saedodin, S., 2017.  Thermal conductivity and viscosity of nanofluids: A review of recent molecular dynamics studies.  Chemical Engineering Science 174, 67–81.  URL: https://www.sciencedirect.com/science/article/pii/S0009250917305432, doi:https://doi.org/10.1016/j.ces.2017.08.034.
  * Ko et al. (2021) Ko, T.W., Finkler, J.A., Goedecker, S., Behler, J., 2021.  A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer.  Nature Communications 12, 398.  URL: https://doi.org/10.1038/s41467-020-20427-2, doi:10.1038/s41467-020-20427-2.
  * Kresse and Furthmüller (1996) Kresse, G., Furthmüller, J., 1996.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  Phys. Rev. B 54, 11169–11186.  URL: https://link.aps.org/doi/10.1103/PhysRevB.54.11169, doi:10.1103/PhysRevB.54.11169.
  * Kresse and Joubert (1999) Kresse, G., Joubert, D., 1999.  From ultrasoft pseudopotentials to the projector augmented-wave method.  Phys. Rev. B 59, 1758–1775.  URL: https://link.aps.org/doi/10.1103/PhysRevB.59.1758, doi:10.1103/PhysRevB.59.1758.
  * Kubo et al. (1957) Kubo, R., Yokota, M., Nakajima, S., 1957.  Statistical-mechanical theory of irreversible processes. ii. response to thermal disturbance.  Journal of the Physical Society of Japan 12, 1203–1211.
  * Lam et al. (2021) Lam, S.T., Li, Q.J., Ballinger, R., Forsberg, C., Li, J., 2021.  Modeling lif and flibe molten salts with robust neural network interatomic potential.  ACS Applied Materials & Interfaces 13, 24582–24592.  URL: https://doi.org/10.1021/acsami.1c00604, doi:10.1021/acsami.1c00604. pMID: 34019760.
  * LeBlanc (2010) LeBlanc, D., 2010.  Molten salt reactors: A new beginning for an old idea.  Nuclear Engineering and Design 240, 1644–1656.  URL: https://www.sciencedirect.com/science/article/pii/S0029549310000191, doi:https://doi.org/10.1016/j.nucengdes.2009.12.033.
  * Lee et al. (2021) Lee, S.C., Zhai, Y., Li, Z., Walter, N.P., Rose, M., Heuser, B.J., Z, Y., 2021.  Comparative studies of the structural and transport properties of molten salt flinak using the machine-learned neural network and reparametrized classical forcefields.  The Journal of Physical Chemistry B 125, 10562–10570.  URL: https://doi.org/10.1021/acs.jpcb.1c05608, doi:10.1021/acs.jpcb.1c05608. pMID: 34496565.
  * Li et al. (2021) Li, Q.J., Küçükbenli, E., Lam, S., Khaykovich, B., Kaxiras, E., Li, J., 2021.  Development of robust neural-network interatomic potential for molten salt.  Cell Reports Physical Science 2, 100359.  URL: https://www.sciencedirect.com/science/article/pii/S2666386421000448, doi:https://doi.org/10.1016/j.xcrp.2021.100359.
  * Liang et al. (2020) Liang, W., Lu, G., Yu, J., 2020.  Molecular dynamics simulations of molten magnesium chloride using machine-learning-based deep potential.  Advanced Theory and Simulations 3, 2000180.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adts.202000180, doi:https://doi.org/10.1002/adts.202000180.
  * Lizin et al. (2017) Lizin, A., Tomilin, S., Ponomarev, L.I., Fedorov, Y.S., Hirose, Y., 2017.  12 - fast-spectrum, liquid-fueled reactors, in: Dolan, T.J. (Ed.), Molten Salt Reactors and Thorium Energy. Woodhead Publishing, pp. 375–433.  URL: https://www.sciencedirect.com/science/article/pii/B9780081011263000129, doi:https://doi.org/10.1016/B978-0-08-101126-3.00012-9.
  * Locatelli et al. (2013) Locatelli, G., Mancini, M., Todeschini, N., 2013.  Generation iv nuclear reactors: Current status and future prospects.  Energy Policy 61, 1503–1520.  URL: https://www.sciencedirect.com/science/article/pii/S0301421513006083, doi:https://doi.org/10.1016/j.enpol.2013.06.101.
  * Lu et al. (2021) Lu, J., Yang, S., Pan, G., Ding, J., Liu, S., Wang, W., 2021.  Thermal and transport properties of molten chloride salts with polarization effect on microstructure.  Energies 14.  URL: https://www.mdpi.com/1996-1073/14/3/746, doi:10.3390/en14030746.
  * Magnusson et al. (2020) Magnusson, J., Memmott, M., Munro, T., 2020.  Review of thermophysical property methods applied to fueled and un-fueled molten salts.  Annals of Nuclear Energy 146, 107608.  URL: https://www.sciencedirect.com/science/article/pii/S0306454920303066, doi:https://doi.org/10.1016/j.anucene.2020.107608.
  * Maslennikova et al. (2023) Maslennikova, A.A., Mushnikov, P.N., Dub, A.V., Tkacheva, O.Y., Zaikov, Y.P., Liu, Y.L., Shi, W.Q., 2023.  Determination of the oxygen content in the lif-naf-kf melt.  Materials 16.  URL: https://www.mdpi.com/1996-1944/16/11/4197, doi:10.3390/ma16114197.
  * Maxwell (2022) Maxwell, C.I., 2022.  Molecular dynamics study of fission gas behaviour and solubility in molten flinak salt.  Journal of Nuclear Materials 563, 153633.  URL: https://www.sciencedirect.com/science/article/pii/S0022311522001295, doi:https://doi.org/10.1016/j.jnucmat.2022.153633.
  * Mullabaev et al. (2022) Mullabaev, A., Kovrov, V., Kholkina, A., Zaikov, Y., 2022.  Anode processes on pt and ceramic anodes in chloride and oxide-chloride melts.  Nuclear Engineering and Technology 54, 965–974.  URL: https://www.sciencedirect.com/science/article/pii/S1738573321005374, doi:https://doi.org/10.1016/j.net.2021.08.034.
  * Müller-Plathe (1997) Müller-Plathe, F., 1997.  A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity.  The Journal of Chemical Physics 106, 6082–6085.  URL: https://doi.org/10.1063/1.473271, doi:10.1063/1.473271.
  * Nosé (1984) Nosé, S., 1984.  A unified formulation of the constant temperature molecular dynamics methods.  The Journal of chemical physics 81, 511–519.
  * Novikov et al. (2020) Novikov, I.S., Gubaev, K., Podryabinkin, E.V., Shapeev, A.V., 2020.  The mlip package: moment tensor potentials with mpi and active learning.  Machine Learning: Science and Technology 2, 025002.  URL: https://dx.doi.org/10.1088/2632-2153/abc9fe, doi:10.1088/2632-2153/abc9fe.
  * Ong et al. (2020) Ong, T.C., Sarvghad, M., Lippiatt, K., Griggs, L., Ryan, H., Will, G., Steinberg, T.A., 2020.  Review of the solubility, monitoring, and purification of impurities in molten salts for energy storage in concentrated solar power plants.  Renewable and Sustainable Energy Reviews 131, 110006.  URL: https://www.sciencedirect.com/science/article/pii/S1364032120302975, doi:https://doi.org/10.1016/j.rser.2020.110006.
  * Pan et al. (2021) Pan, G., Ding, J., Chen, P., Yan, H., Du, Y., Lee, D.J., Lu, Y., 2021.  Finite-size effects on thermal property predictions of molten salts.  Solar Energy Materials and Solar Cells 221, 110884.  URL: https://www.sciencedirect.com/science/article/pii/S0927024820304827, doi:https://doi.org/10.1016/j.solmat.2020.110884.
  * Parasotchenko et al. (2023) Parasotchenko, Y., Gevel, T., Pavlenko, O.B., Gorshkov, L.V., Leonova, N.M., Suzdaltsev, A.V., Zaikov, Y.P., 2023.  Choice of the composition of the chloride melts for the electrochemical synthesis of silicon.  Silicon URL: https://api.semanticscholar.org/CorpusID:264982487.
  * Perdew et al. (1996) Perdew, J.P., Burke, K., Ernzerhof, M., 1996.  Generalized gradient approximation made simple.  Phys. Rev. Lett. 77, 3865–3868.  URL: https://link.aps.org/doi/10.1103/PhysRevLett.77.3865, doi:10.1103/PhysRevLett.77.3865.
  * Podryabinkin and Shapeev (2017) Podryabinkin, E.V., Shapeev, A.V., 2017.  Active learning of linearly parametrized interatomic potentials.  Computational Materials Science 140, 171–180.  URL: https://www.sciencedirect.com/science/article/pii/S0927025617304536, doi:https://doi.org/10.1016/j.commatsci.2017.08.031.
  * Porter et al. (2022) Porter, T., Vaka, M.M., Steenblik, P., Della Corte, D., 2022.  Computational methods to simulate molten salt thermophysical properties.  Communications Chemistry 5, 69.  URL: https://doi.org/10.1038/s42004-022-00684-6, doi:10.1038/s42004-022-00684-6.
  * Robertson et al. (2022) Robertson, S.G., Wiser, R., Yang, W., Kang, D., Choi, S., Baglietto, E., Short, M.P., 2022.  The curious temperature dependence of fluoride molten salt thermal conductivity.  Journal of Applied Physics 131.
  * Rodriguez et al. (2021) Rodriguez, A., Lam, S., Hu, M., 2021.  Thermodynamic and transport properties of lif and flibe molten salts with deep learning potentials.  ACS Applied Materials & Interfaces 13, 55367–55379.  URL: https://doi.org/10.1021/acsami.1c17942, doi:10.1021/acsami.1c17942. pMID: 34767334.
  * Rollet et al. (2009) Rollet, A.L., Sarou-Kanian, V., Bessada, C., 2009.  Measuring self-diffusion coefficients up to 1500 k: A powerful tool to investigate the dynamics and the local structure of inorganic melts.  Inorganic Chemistry 48, 10972–10975.  URL: https://doi.org/10.1021/ic9010086, doi:10.1021/ic9010086. pMID: 19863113.
  * Romatoski and Hu (2017) Romatoski, R., Hu, L., 2017.  Fluoride salt coolant properties for nuclear reactor applications: A review.  Annals of Nuclear Energy 109, 635–647.  URL: https://www.sciencedirect.com/science/article/pii/S0306454917301391, doi:https://doi.org/10.1016/j.anucene.2017.05.036.
  * Rudenko et al. (2022a) Rudenko, A., Kataev, A., Tkacheva, O., 2022a.  Dynamic viscosity of the naf-kf-ndf3 molten system.  Materials 15.  URL: https://www.mdpi.com/1996-1944/15/14/4884, doi:10.3390/ma15144884.
  * Rudenko et al. (2022b) Rudenko, A., Redkin, A., Il’ina, E., Pershina, S., Mushnikov, P., Zaikov, Y., Kumkov, S., Liu, Y., Shi, W., 2022b.  Thermal conductivity of flinak in a molten state.  Materials 15.  URL: https://www.mdpi.com/1996-1944/15/16/5603, doi:10.3390/ma15165603.
  * Salanne et al. (2009) Salanne, M., Simon, C., Turq, P., Madden, P.A., 2009.  Heat-transport properties of molten fluorides: Determination from first-principles.  Journal of Fluorine Chemistry 130, 38–44.  URL: https://www.sciencedirect.com/science/article/pii/S0022113908002078, doi:https://doi.org/10.1016/j.jfluchem.2008.07.013. fluorine and Nuclear Energy.
  * Shapeev (2016) Shapeev, A.V., 2016.  Moment tensor potentials: A class of systematically improvable interatomic potentials.  Multiscale Modeling & Simulation 14, 1153–1173.  URL: https://doi.org/10.1137/15M1054183, doi:10.1137/15M1054183.
  * Sivaraman et al. (2020) Sivaraman, G., Krishnamoorthy, A.N., Baur, M., Holm, C., Stan, M., Csányi, G., Benmore, C., Vázquez-Mayagoitia, Á., 2020.  Machine-learned interatomic potentials by active learning: amorphous and liquid hafnium dioxide.  npj Computational Materials 6, 104.  URL: https://doi.org/10.1038/s41524-020-00367-7, doi:10.1038/s41524-020-00367-7.
  * Smirnov et al. (1987) Smirnov, M., Khokhlov, V., Filatov, E., 1987.  Thermal conductivity of molten alkali halides and their mixtures.  Electrochimica Acta 32, 1019–1026.  URL: https://www.sciencedirect.com/science/article/pii/0013468687900272, doi:https://doi.org/10.1016/0013-4686(87)90027-2.
  * Sohal et al. (2010) Sohal, M., Ebner, M., Sabharwall, P., Sharpe, P., 2010.  Engineering database of liquid salt thermophysical and thermochemical properties.  Idaho Natl. Lab. Idaho Falls CrossRef doi:10.2172/980801.
  * Sun et al. (2016) Sun, J., Remsing, R.C., Zhang, Y., Sun, Z., Ruzsinszky, A., Peng, H., Yang, Z., Paul, A., Waghmare, U., Wu, X., Klein, M.L., Perdew, J.P., 2016.  Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional.  Nature Chemistry 8, 831–836.  URL: https://doi.org/10.1038/nchem.2535, doi:10.1038/nchem.2535.
  * Sun et al. (2015) Sun, J., Ruzsinszky, A., Perdew, J.P., 2015.  Strongly constrained and appropriately normed semilocal density functional.  Phys. Rev. Lett. 115, 036402.  URL: https://link.aps.org/doi/10.1103/PhysRevLett.115.036402, doi:10.1103/PhysRevLett.115.036402.
  * Thompson et al. (2022) Thompson, A.P., Aktulga, H.M., Berger, R., Bolintineanu, D.S., Brown, W.M., Crozier, P.S., in ’t Veld, P.J., Kohlmeyer, A., Moore, S.G., Nguyen, T.D., Shan, R., Stevens, M.J., Tranchida, J., Trott, C., Plimpton, S.J., 2022.  Lammps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Computer Physics Communications 271, 108171.  URL: https://www.sciencedirect.com/science/article/pii/S0010465521002836, doi:https://doi.org/10.1016/j.cpc.2021.108171.
  * Tisi et al. (2021) Tisi, D., Zhang, L., Bertossa, R., Wang, H., Car, R., Baroni, S., 2021.  Heat transport in liquid water from first-principles and deep neural network simulations.  Physical Review B 104, 224202.
  * Umesaki et al. (1981) Umesaki, N., Iwamoto, N., Tsunawaki, Y., Ohno, H., Furukawa, K., 1981.  Self-diffusion of lithium, sodium, potassium and fluorine in a molten lif + naf + kf eutectic mixture.  J. Chem. Soc., Faraday Trans. 1 77, 169–175.  URL: http://dx.doi.org/10.1039/F19817700169, doi:10.1039/F19817700169.
  * Vriesema (1979) Vriesema, B., 1979.  Aspects of molten fluorides as heat transfer agents for power generation.
  * W. D. Powers and Greene (1963) W. D. Powers, S.I.C., Greene, N.D., 1963.  Physical properties of molten reactor fuels and coolants.  Nuclear Science and Engineering 17, 200–211.  URL: https://doi.org/10.13182/NSE63-5, doi:10.13182/NSE63-5.
  * Williams (2006) Williams, D., 2006.  Assessment of candidate molten salt coolants for the ngnp/nhi heat-transfer loop.
  * Zhang et al. (2019) Zhang, L., Lin, D.Y., Wang, H., Car, R., E, W., 2019.  Active learning of uniformly accurate interatomic potentials for materials simulation.  Phys. Rev. Mater. 3, 023804.  URL: https://link.aps.org/doi/10.1103/PhysRevMaterials.3.023804, doi:10.1103/PhysRevMaterials.3.023804.
  * Zhitkov et al. (2022) Zhitkov, A., Potapov, A., Karimov, K., Kholkina, A., Shishkin, V., Dedyukhin, A., Zaykov, Y., 2022.  Interaction between un and cdcl2 in molten licl–kcl eutectic. ii. experiment at 1023 k.  Nuclear Engineering and Technology 54, 653–660.  URL: https://www.sciencedirect.com/science/article/pii/S1738573321005192, doi:https://doi.org/10.1016/j.net.2021.08.030.
  * Zhitkov et al. (2020) Zhitkov, A., Potapov, A., Karimov, K., Shishkin, V., Dedyukhin, A., Zaykov, Y., 2020.  Interaction between un and cdcl2 in molten licl-kcl eutectic. i. experiment at 773 k.  Nuclear Engineering and Technology 52, 123–134.  URL: https://www.sciencedirect.com/science/article/pii/S1738573319301780, doi:https://doi.org/10.1016/j.net.2019.07.006.
  * Zuo et al. (2020) Zuo, Y., Chen, C., Li, X., Deng, Z., Chen, Y., Behler, J., Csányi, G., Shapeev, A.V., Thompson, A.P., Wood, M.A., Ong, S.P., 2020.  Performance and cost assessment of machine learning interatomic potentials.  The Journal of Physical Chemistry A 124, 731–745.  URL: https://doi.org/10.1021/acs.jpca.9b08723, doi:10.1021/acs.jpca.9b08723. pMID: 31916773.


††thanks: Corresponding author<EMAIL_ADDRESS>
# Local electronic structure of dilute hydrogen in $\beta$-MnO2

H. Okabe  Muon Science Laboratory and Condensed Matter Research Center,
Institute of Materials Structure Science, High Energy Accelerator Research
Organization (KEK-IMSS), Tsukuba, Ibaraki 305-0801, Japan  M. Hiraishi  Muon
Science Laboratory and Condensed Matter Research Center, Institute of
Materials Structure Science, High Energy Accelerator Research Organization
(KEK-IMSS), Tsukuba, Ibaraki 305-0801, Japan  A. Koda  Muon Science Laboratory
and Condensed Matter Research Center, Institute of Materials Structure
Science, High Energy Accelerator Research Organization (KEK-IMSS), Tsukuba,
Ibaraki 305-0801, Japan  Department of Materials Structure Science, The
Graduate University for Advanced Studies (Sokendai), Tsukuba, Ibaraki
305-0801, Japan  S. Takeshita  Muon Science Laboratory and Condensed Matter
Research Center, Institute of Materials Structure Science, High Energy
Accelerator Research Organization (KEK-IMSS), Tsukuba, Ibaraki 305-0801, Japan
Department of Materials Structure Science, The Graduate University for
Advanced Studies (Sokendai), Tsukuba, Ibaraki 305-0801, Japan  K. M. Kojima
Centre for Molecular and Materials Science TRIUMF, Vancouver, B.C. V6T 2A3,
Canada  I. Yamauchi  Graduate School of Science and Engineering, Saga
University, Saga 840-8502, Japan  T. Ohsawa  National Institute for Materials
Science (NIMS), Tsukuba, Ibaraki 305-0044, Japan  N. Ohashi  National
Institute for Materials Science (NIMS), Tsukuba, Ibaraki 305-0044, Japan  H.
Sato  Department of Physics, Chuo University, Tokyo 112-8551, Japan  R. Kadono
Muon Science Laboratory and Condensed Matter Research Center, Institute of
Materials Structure Science, High Energy Accelerator Research Organization
(KEK-IMSS), Tsukuba, Ibaraki 305-0801, Japan  Department of Materials
Structure Science, The Graduate University for Advanced Studies (Sokendai),
Tsukuba, Ibaraki 305-0801, Japan

###### Abstract

The electronic and magnetic states of $\beta$-MnO2 in terms of hydrogen
impurities have been investigated by muon spin rotation ($\mu$SR) technique
combined with density-functional theory (DFT) calculations for muon as pseudo-
hydrogen. We found that 85 % of implanted muons are localized in the oxygen
channels of the rutile structure and behave as interstitial protons
(Mu${}^{+})$ except those (7.6%) forming a charge-neutral state (Mu${}^{0})$
at 2.3 K, which indicates that interstitial hydrogen acts as a shallow donor
within less than 0.1 meV of ionization energy. The residual 15% of muons are
attributed to those related to lattice imperfection as Mn vacancies. Detailed
analyses combined with DFT approach suggested that the muon is localized at
the center of the oxygen channel due to its large zero-point vibration energy.

##  I Introduction

Manganese dioxide (MnO${}_{2})$ is well known as a catalyst for a variety of
chemical reactions such as the deoxidization process of hydrogen peroxide,
with which we are familiar in elementary school chemistry, or most recently,
it is attracting a lot of attention as a highly efficient heterogeneous
catalyst for the oxidation of various types of substrates including bio-mass-
derived compounds cat . The most known applications of MnO2 are in batteries;
MnO2 synthesized by the electrolytic method ($\gamma$-MnO${}_{2})$ is used
worldwide as the cathode material in lithium or zinc primary battery bat . The
artificially synthesized MnO2 generally contains hydrogen in the crystal
structure, which are supposed to influence various physical, chemical, and
electrochemical properties.

There are two types of hydrogens (protons) existing in $\gamma$-MnO2 proposed
by Ruetschi and Giovanoli Rucy1 ; Rucy2 ; Colm . One is that protons
associated with cation vacancies, termed ”Ruetschi” protons, which compensate
for the charge instead of Mn4+ cations. The other is that protons associated
with Mn3+ cations, termed ”Coleman” protons, which are located in the tunnels
of the structure and more mobile than Ruetschi protons. The existence of the
two types of protons has been confirmed experimentally by neutron diffraction
and NMR studies Filax ; JTES ; Pitel ; JMS ; NMR . However, these protons are
mainly distinguished by differences in O-H bond lengths, therefore its
specific positions and electronic states remain unclear. The main reason for
the problem is the difficulty of detecting dilute amounts of hydrogen.

Since hydrogen is the lightest and smallest of all the elements, and
ubiquitous presence in nature, hydrogen easily penetrates into materials to
cause embrittlement of metals, or serves as possible sources for unintentional
carrier doping in semiconductors Van1 ; Van2 . Nevertheless, it is difficult
to find out the behavior of dilute hydrogen in matter due to its own weak
signals and/or false detections of background-originated signals from
extraneous hydrogen in various diffraction analyses. It is therefore desirable
to develop an effective method to ascertain the behavior of hydrogen in host
materials.

To obtain microscopic information on the local electronic structure of
isolated hydrogen in matter, positive muon studies have been making
significant contributions. While positive muon as pseudo-hydrogen exhibits a
remarkable isotope effect for kinetics of muon/proton in matter, its local
electronic structure as an atom determined by muon-electron interaction
(resulting in various charged states designated as Mu+, Mu0, or Mu${}^{-})$ is
almost identical with that of hydrogen. Thus, implanted muon can be regarded
as a simulator for the electronic structure of interstitial hydrogen. We
previously identified the electronic structure of muons in oxide
semiconductors using muon spin rotation ($\mu$SR) technique, and found
evidence for hydrogen impurities as a possible origin of unintentional
$n$-type conductivity Gil ; Cox ; Shimo1 ; Salman ; Shimo2 ; Vilao ; C12A7 ;
FeS2 ; IGZO .

As a part of our continued effort to elucidate the behavior of hydrogen in
oxide semiconductors, we focused on $\beta$-MnO2 as our next research target.
MnO2 is known to exhibit a variety of crystalline polymorphs, which results in
various tunnel and layered structures polym , where $\beta$-MnO2 (pyrolusite)
is the simplest and thermodynamically stable phase. The crystal structure of
$\beta$-MnO2 consists of MnO6 octahedra shared by corners or edges with 1
$\times$ 1 tunnels (oxygen channels, see Fig. 1), which is also part of the
basic structure of $\gamma$-MnO2. It would be interesting to note that the
structure of $\beta$-MnO2 is identical with rutile (TiO${}_{2})$, which we
have intensively studied by $\mu$SR Shimo2 ; Vilao .

$\beta-$MnO2 is regarded as a narrow gap semiconductor and it is known to
exhibt a $n$-type conductivity Noda ; Sato . The $n$-type nature of
$\beta$-MnO2 may stem from oxygen deficiency and/or lesser impurity
interstitials like hydrogen. If oxygen vacancy concentration is negligibly
small, hydrogen could be a major source of carriers. However it is often
difficult to separate all these different contributions in the bulk property
measurements. On the other hand, $\mu$SR can provide useful information to
extract intrinsic electronic properties out of those extrinsic contributions
because muon monitors local magnetic states in the atomic scale. Even more
conveniently, $\beta$-MnO2 is known as a classical example of magnetic
materials with a well-defined screw magnetic order below $T_{N}=$ 92 K SatoX ;
Regls . Since muon is a direct probe of local magnetic fields, we can
determine the position and electronic state of hydrogen in $\beta$-MnO2 by
$\mu$SR, using the magnetic structure as a guide.

In this work, we used muon as a microscopic simulator for extrinsic
hydrogen/proton to identify the local electronic and magnetic structures
stemmed from interstitial hydrogen in $\beta$-MnO2, by probing the local
fields from surrounding Mn ions that exhibit the well-defined magnetic order.
The paper proceeds as follows. Section II describes the sample
characterization, $\mu$SR experimental, and computational details. Section III
presents the results of transverse-field (TF) $\mu$SR, zero-field (ZF) $\mu$SR
measurements, and density functional theory (DFT) calculations. Then, we
discuss the positions of hydrogen and the effect on the electronic and
magnetic state of $\beta$-MnO2. Finally, Sec. IV summarizes this paper.

Figure 1: Crystal structure and electrostatic potential maps for $\beta$-MnO2
calculated by using reported lattice parameters Baur . The Hartree potential
calculated using VASP are represented by the color contour from blue to red in
the lattice planes (002) and (100). The solid lines indicate the unit cell.

##  II Experimental Details

Powder sample of $\beta$-MnO2 (3N, Soekawa Co., Ltd.) was used in the present
study. The crystal structure and phase purity were checked using a Rigaku
MultiFlex powder X-ray diffractometer under ambient conditions with Cu
K$\alpha$ radiation. The XRD data indicate that the sample crystallized into
$\beta$-MnO2 single phase and the absence of secondary phase (see Supplemental
Material Fig. S1).

The hydrogen content in the sample was observed by thermal desorption
spectrometry (TDS) placed in the National Institute for Material Science
(NIMS). The sample was linearly heated from room temperature to 1073 K in 2
hours and kept at 1073 K for 1 hour. During the heating, the amount of
desorbed hydrogen was measured by a quadrupole mass spectrometer (QMS). The
QMS signal intensity of hydrogen molecule from sample increase with increasing
temperature, and 2.918 $\times$ 10${}^{19\,}$cm-3 of hydrogen was detected
(Fig. S2).

The magnetic susceptibility $\chi$ was measured using a Quantum Design SQUID
magnetometer from 5 to 400 K. The temperature dependence of $\chi$ shows a
kink at 92 K, corresponding to the Néel temperature of $\beta$-MnO2 (Fig. S3).
The peak anomaly around 40 K is probably due to hydrogen impurities, since it
was suppressed after heat treatment at 723 K in air. $\chi$ also shows the
Curie-Weiss behavior with the Curie constant $C=$ 2.51 emu k/mol and the Weiss
temperature $\Theta=-$788 K, which are in good agreement with the earlier
report Sato  (Fig. S3, inset).

Conventional $\mu$SR experiments were performed using the ARTEMIS spectrometer
installed in the S1 area at Muon Science Establishment (MUSE), Japan Proton
Accelerator Research Complex (J-PARC). Time dependent $\mu$-e decay asymmetry
A($t)$ was measured to monitor the local field under zero field (ZF),
longitudinal field (LF), and weak transverse field (TF) conditions. Additional
$\mu$SR experiments were conducted to resolve the detailed local field
distribution using two instruments, i.e., the LAMPF spectrometer on the M20
beamline and NuTime spectrometer on the M15 beamline at TRIUMF, Canada, for
the magnetic and paramagnetic phases of $\beta$-MnO2, respectively. The
$\mu$SR spectra have been analyzed using the MUSRFIT software package Suter .

The DFT calculations were performed using the projector augmented wave
approach PAW  implemented in the Vienna ab initio simulation package (VASP)
VASP  with the Perdew-Burke-Ernzerhof (PBE) exchange correlation potential PBE
. The cutoff energy for the plane-wave basis set was 500 eV. All atoms were
relaxed until the Hellmann$-$Feynman forces on them were smaller than 0.01
eV/Å. The distribution of the local magnetic field at the muon sites was
calculated using Dipelec program Koji . Crystal structures were visualized
using the VESTA program Izumi .

##  III Results and discussion

###  III.1 Local fields in the paramagnetic phase

At first, we attempted to narrow down the position of hydrogen/muon in
$\beta$-MnO2 structure by examining the local field in the paramagnetic phase.
The local field provides useful information for identifying the muon stopping
site, because it is mainly determined by the spatial arrangement of the
nearest neighboring (nn) Mn $d$-electrons. As shown in Fig. 1, it is suggested
from our preliminary calculation using VASP that the Hartree potential for the
interstitial Mu+ exhibits minima around the center of the oxygen channels
along the $c$-axis direction. When an external magnetic field is applied
transverse to the initial muon polarization vector, the muon exhibits spin
precession at a frequency corresponding to the summation of the external and
internal fields. Thus, by evaluating the shift from the reference frequency of
the external field, one can extract the contribution of the internal field at
a particular stopping site from the TF-$\mu$SR spectrum.

Figure 2 shows the fast Fourier transform (FFT) of TF-$\mu$SR spectra under a
transverse field of 6 T obtained at various temperatures. One can see a single
Gaussian-like lineshape at around 813 MHz, where the slight broadening of
linewidth with decreasing temperature is probably due to the gradual
development of staggered magnetization related to the short range magnetic
order Sato ; Regls . To interpret the lineshape of these spectra, we
calculated frequency distribution in polycrystalline sample (powder pattern)
for several candidate muon sites Slichter .

Figure 2: Fast Fourier transform of the $\mu$SR time spectra measured under a
transverse field of 6 T.

The calculated powder patterns of putative muon sites are shown in Fig. 3 with
additional Gaussian broadening due to the local field inhomogeneity and the
truncation of Fourier transforms by the finite time window ($t\leq$ 10
$\mu$s). These muon sites have been selected by considering the electrostatic
potential maps (see Fig. 1) and the distance from nn O ions. The lineshape has
two edges determined by the anisotropy of the hyperfine coupling (A∥ and
A${}_{\mathrm{\bot}})$, which is related with the local $\chi$ induced by the
nn Mn spins. One can see that the muon site near a certain Mn ion, $e.g.$ 2$b$
site (0, 0, 1/2), shows a marked asymmetric feature with broad distribution.
Since the experimental spectra show no signs of these shoulder-like features,
implanted muons tend to be located away from Mn ions, at a high symmetry
position in the Mn sublattice. Comparison of the spectra with that simulated
for the O defect position (vacancy) also infers that muons have little chance
to find O vacancies within their lifetime.

Figure 3: Simulated powder patterns for the putative muon stopping sites in
$\beta$-MnO2 under a transverse field of 6 T. Each interstitial
crystallographic position is represented in a Wycoff letter: 4$c$ (1/2, 0, 0),
8$i$ (0.42667, 0.11185, 0), 16$k$ (0.48954,0.26299, 0.18754), 2$b$ (0, 0,
1/2). The black lines represent the calculated powder patterns. The red lines
represent the additional Gaussian broadening (see the main text) superimposed
on the powder patterns.

Considering previous neutron and NMR studies Filax ; JTES ; Pitel ; JMS ; NMR
, it is reasonable to assume that muons are localized near the center of the
oxygen channel in the rutile structure. Our result further suggests that muons
do not occupy low symmetry positions including oxygen vacancies. On the other
hand, there remains a possibility that muons are localized in Mn vacancies due
to the lack of frequency resolution (see the lineshape at the Mn defect
position in Fig. 3). Therefore, it requires additional information to resolve
this issue (see the next section).

###  III.2 Local fields in the magnetically ordered phase

Figure 4: (a) $\mu$SR time spectrum at 2.5 K under a zero external field. The
solid line represents the result of least-squares fitting (see the main text).
(b) Fast Fourier transform of the time spectrum.$f_{1}$ ($f_{2})$ indicates a
peak frequency. Figure 5: Temperature dependence of (a) the transverse
(longitudinal) asymmetry $A_{T}$ ($A_{L})$, (b) the transverse (longitudinal)
relaxation rate $\lambda_{T}$ ($\lambda_{L})$, (c) the precession frequencies
$\gamma_{\mu}B_{\mathrm{ave}}$ and $\gamma_{\mu}B_{\mathrm{dif}}$, and (d) the
initial phase $\varphi_{b}$ and $\varphi_{c}$. The open circle in (b)
indicates $\lambda_{L}$ obtained at J-PARC MUSE. The solid line in (c)
indicates the result of least-squares curve-fit by the power-law function (see
the main text).

$\beta$-MnO2 undergoes a helical magnetic transition at $T_{N}=$ 92 K with
incommensurate magnetic modulation [propagation vector $q_{m}=$ (0, 0,
$\sim$2/7)] without much temperature dependence SatoX ; Regls . Given that the
magnetic structure is already known, a spontaneous long-range magnetic order
is often helpful to identify the muon sites by comparing the observed field
distribution with that calculated for the candidate sites using the known
magnetic structure.

Figure 4(a) shows a ZF-$\mu$SR time spectrum at 2.5 K for the sample from the
same batch used in the TF measurements. One can see a clear oscillation in the
early time domain (<0.25 $\mu$s), implying the long-range magnetic order. The
FFT spectrum [Fig. 4 (b)] indicates that the oscillating signal is dominated
by two frequency components ($f_{1\,}\simeq$ 97 and $f_{2\,}\simeq$ 124 MHz).
While this feature is usually attributed to the presence of two distinct muon
sites, it turns out to be not necessarily the case in $\beta$-MnO2.

For an incommensurate magnetic structure under zero field, the local field
distribution $D(B_{\mathrm{loc}})$ is given by

$D\left(B_{\mathrm{loc}}\right)=\frac{2}{\pi}\frac{B_{\mathrm{loc}}}{\sqrt{B_{\mathrm{loc}}^{2}-B_{\mathrm{min}}^{2}}\sqrt{B_{\mathrm{max}}^{2}-B_{\mathrm{loc}}^{2}}}\,,$

where it has two characteristic cutoff field ($B_{\mathrm{min}}$ and
$B_{\mathrm{max}})$ Yao . Unfortunately, the muon polarization function in
time domain associated with $D(B_{\mathrm{loc}})$ cannot be obtained
analytically, we adopt an approximate distribution which is called shifted-
Overhauser distribution$D^{\prime}(B_{\mathrm{loc}})$,

$D^{\prime}\left(B_{\mathrm{loc}}\right)=\frac{1}{\pi}\frac{1}{\sqrt{B_{\mathrm{dif}}^{2}-\left(B_{\mathrm{loc}}-B_{\mathrm{ave}}\right)^{2}}}\,,$

where $B_{\mathrm{dif}}=$ ($B_{\mathrm{max}}-B_{\mathrm{min}})$/2 and
$B_{\mathrm{ave}}=$ ($B_{\mathrm{max}}+B_{\mathrm{min}})$/2, is used for the
incommensurate case Amato . According to this approximation, the oscillatory
part of the muon polarization function is written by the product of the zeroth
order Bessel function and the cosine function.

Based on the above approximation, we analyzed the spectrum by following curve-
fit function,

$A\left(t\right)=A_{T}J_{0}\left(\gamma_{\mu}B_{\mathrm{dif}}t+\phi_{\mathrm{b}}\right)\cos{\left(\gamma_{\mu}B_{\mathrm{ave}}t+\phi_{c}\right)\mathrm{exp}\left(-\lambda_{T}t\right)}+A_{L}\mathrm{exp}\left(-\lambda_{L}t\right)$

where the amplitude $A_{T}$ ($A_{L})$ is proportional to the fraction
perpendicular (parallel) to initial polarization direction, $J_{0}$ is a
Bessel function of the first kind, $\gamma_{\mu\,}(=$ 2$\pi$ $\times$ 135.53
MHz/T) is the muon gyromagnetic ratio, $\varphi_{b}$ ($\varphi_{c})$ is the
initial phase of each function, respectively, and $\lambda_{T}$
($\lambda_{L})$ is the transverse (longitudinal) relaxation rate. The solid
curve in Fig. 4(a) shows the result of least-squares curve-fit. The
temperature dependence of respective parameters below $T_{N}$ are shown in
Figs. 5.

The slight decrease of $A_{T}$ at lower temperatures is attributed to the
partial loss of initial muon polarization (which is discussed below). The
development of the magnetic order parameter with decreasing temperature
manifests in the evolution of the frequency $\gamma_{\mu}B_{\mathrm{ave}}$ as
shown in Fig. 5(c). The temperature dependence of
$\gamma_{\mu}B_{\mathrm{ave}}$ shows a critical behavior described by the
power law,
$\gamma_{\mu}B_{\mathrm{ave}}(T)\mathrm{\,\propto\,}(T_{N}-T)^{\beta}$ with
$T_{N}=$ 92.7(2) K and $\beta=$ 0.182(9). The value of critical exponent
$\beta$ is in good agreement with those obtained from the synchrotron x-ray
magnetic scattering SatoX  and the neutron scattering Regls . Although the
phase parameter $\varphi_{c}$ is nearly zero except around $T_{N}$,
$\varphi_{b}$ tends to be shift by $\sim$ $\pi$/4 compared with $\varphi_{c}$.
It may be related to the $\pi$/4 phase shift between the Bessel function and
the cosine function Yao .

The temperature dependence of the observed signal fraction,
[$A_{T}+A_{L}$]/$A_{0}$ (red circle) and the longitudinal fraction,
$A_{L}$/$A_{0}$ (blue circle) are shown in Fig. 6, where $A_{0}$ ($=$ 0.23) is
estimated from the full asymmetry above $T_{N}$. One can clearly see that the
signal fraction below $T_{N}$ is less than 90 % of the total volume fraction.
In contrast, the longitudinal fraction ($z$ component) is almost one third of
the total volume fraction, indicating that all the implanted muons are subject
to isotropic static magnetic order. These facts indicate that there is an
unresolved signal component (missing fraction) which is depolarized within the
dead time ($t$ <0.01 $\mu$s) of the experiment.

The missing fraction at 2.3 K is estimated to be about 15 % of the total
volume fraction. Such a fast depolarization is attributed to the strong
disorder of $D(B_{\mathrm{loc}})$. Considering that the possibility of muons
occupying O vacancy is unlikely as inferred from $D(B_{\mathrm{loc}})$ (at
least above $T_{N}$, see Fig.3), it is reasonable to attribute the origin of
the missing fraction to the muons occuping Mn vacancies that corresponds to
naturally included Ruetschi protons. This also leads to the hypothesis that
the residual 85 % of muons corresponds to interstitial hydrogens (Coleman
protons) which occupies the oxygen channels.

In close examination of the temperature dependence of $A_{L}$/$A_{0}$ in Fig.
6, one can find that it slightly declines from the one third line as
temperature decreases. This is attributed to the formation of paramagnetic Mu0
state which is subject to spin/charge exchange process with diamagnetic state
(Mu${}^{+}+$ e- $\to$ Mu${}^{0})$. The reduction of $A_{L}$ strongly suggests
that the electronic state bound to the interstitial muon is not polaronic
(i.e., the electron spin not being locked along the Mn spins), but rather
atomic to allow spin-triplet ($F=1$) or spin-singlet ($F=0$) state relative to
the muon spin direction. The precession frequency for the $F=0$ state usually
exceeds the limit determined by experimental time resolution Schenck .
Therefore, only residual polarization for the $F=1$ state can be observed,
resulting in a reduction of the longitudinal fraction by half. If Mu0 is
present, the longitudinal asymmetry $A_{L}$ is reproduced by following linear
combination

$A_{\mathrm{L}}=\frac{1}{3}\left(1-{\frac{1}{2}f}_{\mathrm{Mu}}\right)A_{\mathrm{0}},$

where $f_{\mathrm{Mu}}$ is the fractional yield of Mu0.

The temperature dependence of longitudinal signal fraction $f_{\mathrm{L}}$
(see Fig.6 inset) is adequately described by the thermal activation model
$f_{\mathrm{L}}$ =
$f_{\mathrm{d}}$exp$(-E_{\mathrm{a}}$/2$k_{\mathrm{B}}$$T$), where
$f_{\mathrm{d}}$ is the diamagnetic longitudinal fraction ($f_{\mathrm{L}}$ =
$f_{\mathrm{d}}$, when $f_{\mathrm{Mu}}$ = 0), $E_{\mathrm{a}}$ is the
characteristic energy, and $k_{\mathrm{B}}$ is Boltzmann’s constant. A fit of
$f_{\mathrm{L}}$ (blue line in the inset of Fig.6) below 50 K yields
$f_{\mathrm{d}}$ = 0.3420(8) and $E_{\mathrm{a}}$ = 0.0292(8) meV/K.
Considering the value of $f_{\mathrm{L}}$ = 0.316(1) at 2.3 K,
$f_{\mathrm{Mu}}$ is estimated to be 0.0760(4) (which may be a lower bound for
the Mu0 yield). This means that the electrons from the hydrogen impurity
levels are mostly promoted to the conduction band in most temperature ranges.
In other words, the donor level associated with hydrogen is estimated to be
extremely shallow, probably less than 0.1 meV. Since the diamagnetic muons in
the oxygen channels (corresponding to the Coleman proton) can be attributed to
Mu0 in its ionized state, they are presumed to serve as electron donors. For
these reasons, the interstitial hydrogen is expected to be a source of
$n$-type conductivity in $\beta$-MnO2. We note that the hydrogen concentration
(2.918 $\times$ 10${}^{19\,}$cm${}^{-3})$ obtained from our TDS measurement is
comparable with the carrier density ($n=$ 4.8 $\times$
10${}^{19\,}$cm${}^{-3})$ reported in the previous work Sato , which is
perfectly in line with this scenario.

Figure 6: Temperature dependence of the observed signal fraction (red circle)
and the longitudinal signal fraction $f_{\mathrm{L}}$ (blue circle). The inset
shows the enlargement below 50 K with the fitting curve using the thermal
activation model (see main text). The broken line represents the 1/3 of the
total fraction.

###  III.3 First-principles calculations

In order to gain more understanding of hydrogen states in $\beta$-MnO2, the
DFT calculations were performed using VASP for comparisons with the $\mu$SR
experiments. As the first step, we calculated the local field
$B_{\mathrm{loc}}$ at the most energetically stable position 4$c$ (1/2, 0, 0)
with the screw periodicity $q=$ 2$c$*/7 and Mn4+ moment size of 2.34$\mu_{B}$
Regls . Under zero-external field, $B_{\mathrm{loc}}$ is determined by a
vector sum of the magnetic dipolar field of the Mn ions and the Fermi contact
interaction with unpaired electron in the muon 1$s$ orbital Hayano ,

$B_{\mathrm{loc}}=\sum\limits_{i}{\hat{A}_{i}\mu_{i}}+\frac{8}{3}\pi
S_{\mathrm{e}}\delta\left(r\right)$

where $\mu_{i}$ is the magnetic moment of the $i$-th Mn ion located at
distance $r_{i}=$ ($x_{i}$, $y_{i}$, $z_{i})$ from the muon site,
$S_{\mathrm{e}}$ is the electron spin magnetic moment, and $\delta(r)$ is the
delta function, and

$\hat{{\bm{A}}_{i}}=A_{i}^{\alpha\beta}=\frac{1}{r_{i}^{3}}\left(\frac{3\alpha_{i}\beta_{i}}{r^{2}_{i}}-\delta_{\alpha\beta}\right),\>(\alpha,\beta=x,y,z)$

is the dipolar tensor. The Fermi contact term is basically negligible except
for the Mu0 state.

Figure 7(b) shows a simulated distribution of $B_{\mathrm{loc}}$ for the 4$c$
site, which is in remarkable agreement with the experimental result [Fig.
7(a)]. A comparative study with DFT calculation provides useful information
for identifying the muon’s location and state. For example, if all the nn Mn4+
($S=$ 3/2) ions bear the full moment (3$\mu_{B})$, the calculated distribution
of$B_{\mathrm{loc}}$ becomes broader, leading to the shift of central
frequency by $\sim$ 30 MHz (not shown), which is in disagreement with the
experimental result. This also supports the moment reduction of Mn ions
(2.34$\mu_{B})$ from the viewpoint of local magnetism probed by $\mu$SR.

We further extended the simulation of $B_{\mathrm{loc}}$ to investigate the
influence of structural relaxation by a muon/hydrogen insertion. The Brillouin
zone for a 2 $\times$ 2 $\times$ 3 supercell was sampled using a cutoff energy
of 500 eV and a 4 $\times$ 4 $\times$ 4 $k$-point mesh for calculations. The
structure optimization was performed using lattice parameters and all of the
atomic positions were allowed to relax until the forces acting on each ion
became less than 0.01 eV/Å. Note that the spin degree of freedom and the
electron correlation parameter $U$ are not incorporated, because the
difference between optimized atomic positions of Mn around the hydrogen is
small irrespective of the presence or absence of $U$(= 4 eV). Therefore, it is
not an issue in this study.

The result of calculation is shown in Fig. 8. The hydrogen occupies slightly
off-center of the oxygen channel where it forms O-H bond with a bond length of
0.99 Å and donating electron to the host lattice. The surrounding oxygen show
large displacement from their initial crystallographic positions due to
hydrogen bond. On the other hand, the surrounding cationic sub-lattice of Mn
ions shows a slight expansion due to the electrostatic repulsion between H+
and Mn4+ ion; resulting in 1.4$\sim$7.2% Mn-Mn bonds elongation compared to
those without H insertion. Taking into account the displacement of neighboring
six Mn ions from H, we have simulated $B_{\mathrm{loc}}$ using a
686-octahedron 7 $\times$ 7 $\times$ 7 supercell to reproduce the partially
distorted $\beta$-MnO2 lattice considering the period of the helical magnetic
structure.

When muon (or hydrogen) is bound to an oxygen, the donated electron can
localize on a neighboring transition metal $d$ orbital to form a charge-
neutral complex (i.e., polaron) Shimo2 ; Vilao ; Cr2O3 . In light of this
possibility, we simulated $B_{\mathrm{loc}}$ with two different charge
distributions: a polaronic electron distributed on the nn Mn ion only [Fig.
7(c)] or the nearby three Mn ions [Fig. 7(d)], that are not obtained from the
DFT calculation. Although subtle differences in shape are observed in
comparison with the non-distorted case [Fig. 7(b)], they still have a
characteristic two peak structure. Note that when the structure was optimized
by fixing hydrogen at the 4$c$ site, the Mn-H distance was uniformly extended
for the all Mn ions. Therefore, no significant change in the shape of the
local magnetic field was obtained.

Hydrogen generally takes an asymmetric O-H···O configuration with a short O-H
covalent bond and a long H···O hydrogen bond. However, it resides
equidistantly between two nearest oxygen atoms under high pressures. This
phenomenon is called H-bond symmetrization and often observed in high-pressure
phase of ice ice . The H-bond symmetrization is closely linked to zero-point
vibration, thus the isotope effect is expected to be more pronounced Hsym ;
Hira . For example, $\delta$-AlOOH, which has a distorted-rutile-type
framework, shows the H-bond symmetrization under high pressures with a
significant isotope effect; the deuteration induces a shift of the
symmetrization to a higher pressure AlOOH . It is commonly known that a
lighter isotope possesses a larger zero-point energy. The presumption that
muon locates in the middle of two nearest oxygen atoms to form a symmetric
O-$\mu$-O bond due to the isotope effect would be a possible scenario to
explain the muon localization at the 4$c$ site. We have evaluated the zero-
point motion energy of muon and hydrogen from the potential energy
distribution along O-4$c$-O direction by the DFT calculation using VASP (see
Supplemental Material Fig. S4). Given the energy difference between the ground
state ($E_{0}$) and the first excited state ($E_{1}$) of muon or hydrogen,
muon can tunnel between the two wells at a very fast timescale as far as the
muon spin precession frequency is much lower than the tunneling frequency.
This picture can be approximated as a kind of resonant state where muon
locates at the center of the tunnel. In such a case where the quantum nature
is pronounced, muon may not be able to simulate the behavior of dilute
hydrogen defects correctly.

Figure 7: (a) Local field distribution ($B_{loc})$ at 2.3 K under a zero
external field. (b) Simulated distribution of $B_{\mathrm{loc}}$ at the 4$c$
site. Simulated distribution of $B_{\mathrm{loc}}$ at the structurally-
optimized site with a polaronic electron (not obtained from the DFT
calculation) distributed on the nearest Mn ion only (c) or the nearby three Mn
ions (d). Figure 8: Deformation of the host lattice on hydrogen intercalation
inferred from the DFT (VASP) calculation. The broken line represents the H···O
hydrogen bond.

###  III.4 Electronic and magnetic states

Finally, we discuss the electronic and magnetic states of $\beta$-MnO2 based
on the results obtained so far. The formal magnetic moment of Mn4+ ions should
be 3$\mu_{B}$ as the localized spin system, because $\beta$-MnO2 shows an
insulating behavior in the low-temperature region ($T$ <50 K). However, the Mn
ordered moment is markedly reduced to 2.34$\mu_{B}$ Regls , and this fact has
been supported by the present $\mu$SR result as described above.

The discrepancy in the size of the Mn moment may be attributed to the
localized $e_{g}$ electrons originating from unintended impurity defects. If
the Hund coupling is fairly large, $e_{g}$ electrons can become localized at a
part of the Mn sites and perfectly polarized along the localized $t_{2g}$
moments at the expense of their itinerancy Sato . However, massive amounts of
impurities or oxygen defects would be required to induce such significant Mn
moment reduction (3$\mu_{B}\to$ 2.34$\mu_{B})$. Moreover, Hund coupling
usually enhances the magnetic moment, e.g., double-exchange interaction, thus
it would lead to a contradictory result. Although there is a possibility that
the delocalization effect of $t_{2g}$ electrons remains at low temperatures
Sato , we can rule it out because our $\mu$SR results clearly indicate that
the local field distribution of $\beta$-MnO2 is well reproduced by the
summation of the dipolar interaction from localized Mn moments.

Here, we propose an alternative scenario to explain the Mn moment reduction in
$\beta$-MnO2. Charge redistribution from a formal valence between a metal and
its ligand is a fundamental process in transition metal oxides Boc . For
example, there are slightly-trivalent Cu ions (Cu${}^{2+}\underline{L}$,
$\underline{L}$: a hole on the oxygen ligand sites) in hole-doped high-$T_{C}$
superconductors LH  or oxygen holes in Co3+/Co4+ mixed valence state in
LixCoO2 Miz . The degree of the charge transfer from O2$p$ orbital to
transition-metal 3$d$ orbital depends on the energy difference between both
orbitals, and also relates to their electronegativity. In the case of
Mn3+/Mn4+ mixed-valent oxides, significant charge transfer between Mn and O
ions has been observed in the x-ray Compton scattering Cmpt . It is also in
good agreement with the theoretical prediction that Mn–O bond is not purely
ionic but partially covalent DV . Such an intermediate character between ionic
and covalent Mn-O bond may play an important role in the catalytic activity
and the performance for electrochemical storage of manganese oxides.

Based on these facts, we assumed the charge-transfer process from the O2$p$
orbital to the Mn 3$d$-$t_{2g}$ orbital in $\beta$-MnO2. If the covalent
character of Mn-O bond suppresses the increase of Mn3+ ionic radius and
stabilizes a low-spin electronic configuration ($t_{2g}^{4})$ absolutely, the
Mn3+ moment would be significantly reduced (2$\mu_{B})$ compared to that
(4$\mu_{B})$ of a high-spin configuration ($t_{2g}^{3}e_{g}^{1})$. As a
reference, a low spin state of Mn3+ ion is often found in some tunnel-
structure manganese oxides Shannon ; LS . The experimental value of the
reduced Mn moment (2.34$\mu_{B})$ in $\beta$-MnO2 can be reproduced by the
charge transfer of about 0.1 electron from each of the surrounding six oxygen
ions (or even less than this, if the oxygen spin polarization is antiparallel
to Mn spins miya ; ishii ). Further studies of each wave-function character,
e.g., Compton scattering, will be required to support our argument.

##  IV Summary

We have studied the electronic and magnetic state of $\beta$-MnO2 in terms of
hydrogen impurities by $\mu$SR technique combined with the first-principles
calculations. We have found the following experimental results: (1) 85 % of
muons localize in the oxygen channels of the rutile structure as interstitial
protons, (2) residual 15% of muons are associated with lattice imperfection as
Mn vacancies, and (3) 7.6% of those muons (in the oxygen channels) form
charge-neutral Mu state at 2.3 K, which suggests that hydrogen interstitials
act as shallow donors whose ionization energy is less than 1 meV. Detailed
analyses combined with the first-principles calculations showed that implanted
muons may prefer the center of the oxygen channel to form a symmetric
O-$\mu$-O bond due to a large isotope effect.

## Acknowledgment

This work was supported by the MEXT Elements Strategy Initiative to Form Core
Research Centers, from the Ministry of Education, Culture, Sports, Science,
and Technology of Japan (MEXT) under Grant No. JPMXP0112101001. The $\mu$SR
experiments were performed under user programs (Proposal No. 2017A0113 and
2019MS02 at the Materials and Life Science Experimental Facility of the
J-PARC, and M1712 at TRIUMF). We would like to thank the staff of J-PARC and
TRIUMF for their technical support during the $\mu$SR experiment.

## References

  * (1) E. Hayashi, Y. Yamaguchi, K. Kamata, N. Tsunoda, Y. Kumagai, F. Oba and M. Hara, J. Am. Chem. Soc. 141, 890 (2019).
  * (2) D. Glover, B. Schumm, A. Kozawa, Handbook of Manganese Dioxides-Battery Grade, (The International Battery Material Association, Cleveland, 1989).
  * (3) A. F. Wells, Structural Inorganic Chemistry, (Oxford University Press, Oxford, 1975).
  * (4) P. Ruetschi, J. Electrochem. Soc. 131, 2737 (1984).
  * (5) P. Ruetschi and R. Giovanoli, J. Electrochem. Soc. 135, 2663 (1988).
  * (6) J. J. Coleman, Trans. Electrochem. Soc. 90, 545 (1946).
  * (7) F. Fillaux, C. H. Cachet, H. Ouboumour, J. Tomkinson, C. Lévy-Clément, and L. T. Yu, J. Electro-chem. Soc. 140, 585 (1993).
  * (8) M. Nagao, C. Pitteloud, T. Kamiyama, T. Otomo, K. Itoh, T. Fukunaga, and R. Kanno, J. Electro-chem. Soc. 152, E230 (2005).
  * (9) C. Pitteloud, M. Nagao, K. Itoh, and R. Kanno, J. Sol. Sta. Chem. 181, 467 (2008).
  * (10) W. I. Jung, M. Nagao, C. Pitteloud, K. Itoh, A. Yamada, and R. Kanno, J. Mater. Chem. 19, 800 (2009).
  * (11) Y. Paik, J. P. Osegovic, F. Wang, W. Bowden, C. P. Grey, J. Am. Chem. Soc. 123, 9367 (2001).
  * (12) C. G. Van de Walle, Phys. Status Solidi B 235, 89 (2003).
  * (13) C. G. Van de Walle, Phys. Rev. Lett. 85, 1012 (2000).
  * (14) J. M. Gil, H. V. Alberto, R. C. Vilão, J. Piroto Duarte, P. J. Mendes, L. P. Ferreira, N. Ayres de Campos, A. Weidinger, J. Krauser, Ch. Niedermayer, and S. F. J. Cox, Phys. Rev. Lett. 83, 5294 (1999).
  * (15) S. F. J. Cox, E. A. Davis, S. P. Cottrell, P. J. C. King, J. S. Lord, J. M. Gil, H. V. Alberto, Vilão, J. Piroto Duarte, N. Ayresde Campos, A. Weidinger, R. L. Lichti, and S. J. C. Irvine, Phys. Rev. Lett. 86, 2601 (2001).
  * (16) K. Shimomura, K. Nishiyama, and R. Kadono, Phys. Rev. Lett. 89, 255505 (2002).
  * (17) Z. Salman, T. Prokscha, A. Amato, E. Morenzoni, R. Scheuermann, K. Sedlak, and A. Suter, Phys. Rev. Lett. 113, 156801 (2014).
  * (18) K. Shimomura, R. Kadono, A. Koda, K. Nishiyama, and M. Mihara, Phys. Rev. B 92, 075203 (2015).
  * (19) R. C. Vilão, R. B. L. Vieira, H. V. Alberto, J. M. Gil, A. Weidinger, R. L. Lichti, B. B. Baker, P. W. Mengyan, and J. S. Lord, Phys. Rev. B 92, 081202(R) (2015).
  * (20) M. Hiraishi, K. M. Kojima, M. Miyazaki, I. Yamauchi, H. Okabe, A. Koda, R. Kadono, S. Matsuishi, and H. Hosono, Phys. Rev. B 93, 121201(R) (2016).
  * (21) H. Okabe, M. Hiraishi, S. Takeshita, A. Koda, K. M. Kojima, and R. Kadono, Phys. Rev. B 98, 075210 (2018).
  * (22) K. M. Kojima, M. Hiraishi, H. Okabe, A. Koda, R. Kadono, K. Ide, S. Matsuishi, H. Kumomi, T. Kamiya, and H. Hosono, Appl. Phys. Lett. 115, 203901 (2019).
  * (23) Y. Noda, S. Ohno, and S. Nakamura, Phys. Chem. Chem. Phys. 18, 13294 (2016).
  * (24) H. Sato, T. Enoki, M. Isobe and Y. Ueda: Phys. Rev. B 61, 3563 (2000).
  * (25) W. H. Baur, Acta Crystallogr. B 32, 2200 (1976).
  * (26) H. Sato, K. Wakiya, T. Enoki, T. Kiyama, Y. Wakabayashi, H. Nakao and Y. Murakami, J. Phys. Soc. Jpn. 70, 37 (2001).
  * (27) M. Regulski, R. Przeniosło, I. Sosnowska, and J.-U. Hoffmann, J. Phys. Soc. Jpn. 73, 3444 (2004).
  * (28) A. Suter and B. Wojek, “MUSRFIT: A free platform-independent framework for $\mu$SR data analysis,” Phys. Procedia 30, 69–73 (2012).
  * (29) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * (30) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11 169 (1996).
  * (31) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (32) K. M. Kojima, J. Yamanobe, H. Eisaki, S. Uchida, Y. Fudamoto,I. M. Gat, M. I. Larkin, A. Savici, Y. J. Uemura, P. P. Kyriakou, M. T. Rovers, and G. M. Luke, Phys. Rev. B 70, 094402 (2004).
  * (33) K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).
  * (34) C. P. Slichter, Principles of Magnetic Resonance 3rd edition, (Springer-Verlag, New York, 1990).
  * (35) A. Schenck, Muon Spin Rotation Spectroscopy, Principles and Applications in Solid State Physics (Adam Hilger Ltd., Bristol and Boston, 1985), p. 215.
  * (36) R. S. Hayano, Y. J. Uemura, J. Imazato, N. Nishida, K. Nagamine, T. Yamazaki, Y. Ishikawa, and H. Yasuoka, J. Phys. Soc. Jpn. 49, 1773 (1980).
  * (37) A. Yaouanc and P. Dalmas de Réotier, Muon Spin Rotation, Relaxation and Resonance: Applications to Condensed Matter (Oxford University Press, Oxford, 2011).
  * (38) A. Amato, P. Dalmas de Réotier, D. Andreica, A. Yaouanc, A. Suter, G. Lapertot, I. M. Pop, E. Morenzoni, P. Bonfà, F. Bernardini, and R. De Renzi, Phys. Rev. B 89, 184425 (2014).
  * (39) M. H. Dehn, J. K. Shenton, S. Holenstein, Q. N. Meier, D. J. Arseneau, D. L. Cortie, B. Hitti, A. C. Y. Fang, W. A. MacFarlane, R. M. L. McFadden, G. D. Morris, Z. Salman, H. Luetkens, N. A. Spaldin, M. Fechner, and R. F. Kiefl, Phys. Rev. X 10, 011036 (2020).
  * (40) P. Pruzan, J. C. Chervin, and B. Canny, J. Chem. Phys. 99, 9842 (1993).
  * (41) M. Benoit, D. Marx, and M. Parrinello, Nature 392, 258 (1998).
  * (42) M. Hiraishi, K. M. Kojima, H. Okabe, S. Takeshita, A. Koda, R. Kadono, R. Khasanov, S. Iimura, S. Matsuishi, and H. Hosono, Phys. Rev. B 101, 174414 (2020).
  * (43) A. Sano-Furukawa, T. Hattori, K. Komatsu, H. Kagi, T. Nagai, J. J. Molaison, A. M. dos Santos, and C. A. Tulk, Scientific Reports. 8, 15520 (2018).
  * (44) Y. Liu, T. Fujiwara, H. Yukawa, and M. Morinaga, Solid State Ionics 126, 209 (1999).
  * (45) A. E. Bocquet, T. Mizokawa, T. Saitoh, H. Namatame, and A. Fujimori, Phys. Rev. B 46, 3771 (1992).
  * (46) A. Fujimori, E. Takayama-Muromachi, Y. Uchida, and B. Okai, Phys. Rev. B 35, 8814(R) (1987).
  * (47) T. Mizokawa, Y. Wakisaka, T. Sudayama, C. Iwai, K. Miyoshi, J. Takeuchi, H. Wadati, D. G. Hawthorn, T. Z. Regier, and G. A. Sawatzky, Phys. Rev. Lett. 111, 056404 (2013).
  * (48) K. Suzuki, B. Barbiellini, Y. Orikasa, N. Go, H. Sakurai, S. Kaprzyk, M. Itou, K. Yamamoto, Y. Uchimoto, Yung Jui Wang, H. Hafiz, A. Bansil, and Y. Sakurai, Phys. Rev. Lett. 114, 087401(2015).
  * (49) R. D. Shannon and C. T. Prewitt, Acta Cryst. B 25, 925 (1969).
  * (50) A. D. Wadsley, Acta Cryst. 6, 433 (1953).
  * (51) M. Miyazaki, R. Kadono, M. Hiraishi, A. Koda, K. M. Kojima, K. Ohashi, T. Takayama, and H. Takagi, Phys. Rev. B 91, 155113 (2015).
  * (52) Y. Ishii, S. Horio, Y. Noda, M. Hiraishi, H. Okabe, M. Miyazaki, S. Takeshita, A. Koda, K. M. Kojima, R. Kadono, H. Sagayama, H. Nakao, Y. Murakami, and H. Kimura, Phys. Rev. B 101, 224436 (2020).