# Ambipolar Diffusion in the Lower Solar Atmosphere: MHD Simulations of a
Sunspot

Conor D. MacBride Astrophysics Research Centre, School of Mathematics and
Physics, Queen’s University Belfast, Belfast, BT7 1NN, UK  David B. Jess
Astrophysics Research Centre, School of Mathematics and Physics, Queen’s
University Belfast, Belfast, BT7 1NN, UK  Department of Physics and Astronomy,
California State University Northridge, Northridge, CA 91330, USA  Elena
Khomenko Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife,
Spain  Departamento de Astrofísica, Universidad de La Laguna, 38205 La Laguna,
Tenerife, Spain  Samuel D. T. Grant Astrophysics Research Centre, School of
Mathematics and Physics, Queen’s University Belfast, Belfast, BT7 1NN, UK
Conor D. MacBride<EMAIL_ADDRESS>
###### Abstract

Magnetohydrodynamic (MHD) simulations of the solar atmosphere are often
performed under the assumption that the plasma is fully ionized. However, in
the lower solar atmosphere a reduced temperature often results in only the
partial ionization of the plasma. The interaction between the decoupled
neutral and ionized components of such a partially ionized plasma produces
ambipolar diffusion. To investigate the role of ambipolar diffusion in
propagating wave characteristics in the photosphere and chromosphere, we
employ the Mancha3D numerical code to model magnetoacoustic waves propagating
through the atmosphere immediately above the umbra of a sunspot. We solve the
non-ideal MHD equations for data-driven perturbations to the magnetostatic
equilibrium and the effect of ambipolar diffusion is investigated by varying
the simulation to include additional terms in the MHD equations that account
for this process. Analyzing the energy spectral densities for simulations
with/without ambipolar diffusion, we find evidence to suggest that ambipolar
diffusion plays a pivotal role in wave characteristics in the weakly ionized
low density regions, hence maximizing the local ambipolar diffusion
coefficient. As a result, we propose that ambipolar diffusion is an important
mechanism that requires careful consideration into whether it should be
included in simulations, and whether it should be utilized in the analysis and
interpretation of particular observations of the lower solar atmosphere.

##  1 Introduction

It is becoming clear that the degree of plasma ionization in the solar
atmosphere must be accounted for in a wide variety of situations. Plasma with
a particularly large number of neutrals often has significantly modified
properties, resulting in changes to the physical processes that drive
important aspects of solar dynamics (Pandey & Wardle, 2008; Arber et al.,
2007; Ni et al., 2015; Khomenko et al., 2018, 2021; Martínez-Gómez et al.,
2021; Popescu Braileanu & Keppens, 2021), see Ballester et al. (2018) for a
review.

Within the solar atmosphere, plasma at photospheric and chromospheric heights
is partially ionized, with an ionization degree of approximately
$10^{-4}-10^{-1}$ (Ni et al., 2015; Vernazza et al., 1981). This is a result
of the photospheric and chromospheric plasma having lower temperatures and gas
pressure decreasing with height. Accounting for the physical processes
occurring in both the partially ionized and fully ionized regimes, as well as
across their interface, is important for investigating the effect of waves in
heating the corona (Pandey & Wardle, 2008).

The plasma is composed of three species, ions, electrons and neutrals,
interacting via elastic or non-elastic collisions, ion-neutral, electron-
neutral, and electron-ion. All the collision frequencies are a function of
temperature and number density of colliding particles, and their values
strongly depend on height. Partial ionization effects primarily arise due to
the substantial quantity of neutral species in the plasma, and their
insufficient collisional coupling to the charged species, leading to
additional physical processes that need to be taken into consideration. For
example, accounting for the presence of neutrals, and the degree of coupling
between the neutral and charged plasma species, significantly affects the
Biermann battery effect (Biermann, 1950) and the strength of the magnetic
field it generates within the plasma (Martínez-Gómez et al., 2021).

Understanding the effects of partial ionization in the solar atmosphere is
essential for the complete understanding of high-resolution observations. By
comparing the collisional frequencies between the plasma components to the
dynamical frequency of the plasma, the level of coupling between the plasma
components can be determined. Based on this level of coupling, particular
partial ionization effects become significant. As highlighted above, partially
ionized plasma behaves differently to fully ionized environments, hence
employing incorrect plasma conditions in modern numerical simulations may
produce unexpected results and render direct comparisons with observations
inappropriate. This can particularly be seen in the study of the dynamics of
prominences in the solar atmosphere due to Rayleigh-Taylor and Kelvin-
Helmholtz instabilities (Ballester et al., 2020). The abundance of neutrals in
the associated plasma has been shown to affect the growth rates and
instability thresholds, directly impacting the interpretation of any
observation (Soler et al., 2012; Díaz et al., 2012; Khomenko et al., 2014b;
Popescu Braileanu et al., 2021a, b).

In this study we focus on the process of ambipolar diffusion, and its
significance within the chromospheric layer directly above a sunspot.
Ambipolar diffusion occurs within a partially ionized plasma when the neutrals
are not fully coupled to the ionized component. The ionized component is
affected by the Lorentz force, resulting in a magnetic field which is frozen
into those plasma species. However, the decoupled neutral particles do not
experience the Lorentz force, and hence the moving ionized particles are in
constant collision with the neutrals. These collisions result in frictional
effects between the two components, providing a mechanism for magnetic and
mechanical energy to be dissipated, and hence creates a source of localized
atmospheric heating (Khomenko et al., 2018).

Numerical works by Shelyag et al. (2016); Popescu Braileanu & Keppens (2021)
investigated the effect of ambipolar diffusion on chromospheric waves. While
their simulations were constructed to represent the quiet Sun, Popescu
Braileanu & Keppens (2021) concluded that ambipolar diffusion has a
significant damping effect on waves, with areas of greater magnetic field
strength exhibiting more localized damping. Chromospheric waves have also been
studied in realistic solar magneto-convection simulations by Khomenko et al.
(2018). They find that, under the influence of ambipolar diffusion, waves in
the chromosphere have amplitude variations that depend on both frequency and
atmospheric height. Further to this, Khomenko et al. (2021) established that
ambipolar diffusion has the effect of dissipating vortex motions that occur
over small length scales at atmospheric heights above 200 km.

Previous ambipolar diffusion studies have mainly focused on the quiet Sun over
magnetic regions such as sunspot umbrae, where higher magnetic field strengths
and lower densities produce a stronger ambipolar diffusion effect on wave
energy propagation in the chromosphere (Khomenko et al., 2014a). The strong,
aligned magnetic fields above sunspots provide an ideal conduit for wave
propagation, allowing energy to be transferred through the solar atmosphere
(see, e.g., the review by Jess et al., 2015). The influence of ambipolar
diffusion in these regions has already been highlighted by Krishna Prasad et
al. (2017, henceforth referred to as KP17) in their observations of the mid-
photosphere to mid-chromosphere above a sunspot umbra using high-resolution
data obtained by the Rapid Oscillations in the Solar Atmosphere (ROSA; Jess et
al., 2010), the Hydrogen-alpha Rapid Dynamics camera (HARDcam; Jess et al.,
2012), and the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al.,
2014) instruments. KP17 found that intensity oscillations (which were employed
as a proxy for plasma density) in the computed power spectral densities
exhibited a power-law dependence at frequencies beyond the dominant $p$-mode
peak at $\approx 6.5$ mHz. This power-law dependence was observed across a
range of wavelength channels, with the power-law index varying across the
different imaging filters (and hence atmospheric heights), suggesting a
height-dependency in high-frequency wave attenuation due to the changing
physical mechanisms at different heights in the atmosphere. Furthermore, KP17
highlight that ambipolar diffusion may be a factor in more efficient
frequency-dependent damping at higher frequencies. Such studies highlight that
forming a detailed understanding of the effect of ambipolar diffusion in this
domain will allow for more precise interpretations of observational results,
and allow for more realistic models and simulations of sunspot atmospheres.

In this work, we carry out magnetohydrodynamic (MHD) simulations of
magnetoacoustic waves propagating through the atmosphere directly above a
sunspot, and analyze the effects that ambipolar diffusion has on the wave
properties by comparing simulations with and without ambipolar effects, using
identical starting conditions. In Section 2 we introduce the numerical code
employed and in Section 3 we detail the atmospheric model that was employed in
the simulation. In Section 4 we outline the results from our numerical
experiments and provide further insights and conclusions in Section 5.

##  2 Numerical solution

For this work we are using the Mancha3D code (Khomenko & Collados, 2006;
Felipe et al., 2010; Khomenko & Collados, 2012), which solves the three-
dimensional non-linear non-ideal MHD equations for perturbations to the
magnetostatic equilibrium. In particular, Mancha3D can be configured to either
include (or exclude) ambipolar diffusion terms in the embedded MHD equations.
When written in their conservative form, the equations Mancha3D is configured
to solve in our simulations are,

$\dfrac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\mathbf{v}\right)=0\ ,$
(1)

$\dfrac{\partial{\left(\rho\mathbf{v}\right)}}{\partial
t}+\nabla\cdot\left[\rho\mathbf{v}\mathbf{v}+\left(p+\dfrac{|\mathbf{B}|^{2}}{2\mu_{0}}\right)\mathbf{I}-\dfrac{\mathbf{B}\mathbf{B}}{\mu_{0}}\right]=\rho\mathbf{g}+\mathbf{S}(t)\
,$ (2)

$\dfrac{\partial e_{\textrm{int}}}{\partial
t}+\nabla\cdot\mathbf{v}e_{\textrm{int}}+(\gamma-1)e_{\textrm{int}}\nabla\mathbf{v}=\eta_{A}{|\mathbf{J}_{\perp}|}^{2}\
,$ (3)

$\dfrac{\partial\mathbf{B}}{\partial
t}=\nabla\times\left(\mathbf{v}\times\mathbf{B}-\eta_{A}\mathbf{J}_{\perp}\right)\
,$ (4)

where $t$ is time, $\rho$ is the density, $\mathbf{v}$ is the velocity, $p$ is
the gas pressure, $\mu_{0}$ is the vacuum permeability, $\mathbf{B}$ is the
magnetic field, $\mathbf{I}$ is the identity tensor, $\mathbf{g}$ is the
gravitational acceleration, $\mathbf{S}(t)$ is the time-dependent driving
force applied during the simulation, and $e_{\textrm{int}}$ is the internal
energy. In our experiments we used an ideal equation of state (EOS) and
therefore $e_{\textrm{int}}=p/(\gamma-1)$, with the heat capacity ratio
$\gamma=5/3$. The ideal EOS neglects the effects of ionization and molecular
dissociation. The main ionizing species in the chromosphere is hydrogen and
its efficient ionization starts after reaching temperatures of about 10 kK
(Saha, 1921). In our model we only reach temperatures of about 7 kK in the
uppermost part of the domain, therefore, the effects of ionization in the EOS
should not compromise our results, while allowing for the calculations to be
simplified. Later in the paper it will be shown that the effect of ambipolar
diffusion peaks at even lower heights, where the hydrogen stays mostly
neutral. Throughout our numerical calculations we use a constant mean
molecular weight of $\mu=1.28$, corresponding to a neutral solar gas mixture.

The current density, $\mathbf{J}$, is given by
$\mathbf{J}=\nabla\times\mathbf{B}/\mu_{0}$, and $\mathbf{J}_{\perp}$ is the
current density in the direction perpendicular to the magnetic field
$\mathbf{B}$ and is given by,

$\mathbf{J}_{\perp}=-\dfrac{(\mathbf{J}\times\mathbf{B})\times\mathbf{B}}{|\mathbf{B}|^{2}}\
.$ (5)

Each of these equations will also contain an additional artificial diffusion
term for numerical stability.

The ambipolar diffusion coefficient, $\eta_{A}$, is,

$\eta_{A}=\dfrac{\xi^{2}_{n}|\mathbf{B}|^{2}}{\alpha_{n}}\ ,$ (6)

where $\xi_{n}$ is the fraction of neutrals, and $\alpha_{n}$ is the neutral
collisional parameter. The units of $\eta_{A}$ are $\mu_{0}$ m2 s-1. For
simulations without ambipolar diffusion present, all terms with $\eta_{A}$ are
excluded from the MHD equations employed.

In the MHD equations, the $\eta_{A}$ terms are always associated with
$\mathbf{J}_{\perp}$ terms. As a result, ambipolar diffusion only occurs when
there is a significant electric current in the direction perpendicular to the
magnetic field. This anisotropy is due to the net diffusion velocity of all
plasma components, while only the ions are are affected by the Lorentz force
(Cowling, 1957; Braginskii, 1965). The ambipolar diffusion coefficient is
related to the Cowling resistivity, $\eta_{C}=\eta+\eta_{A}$, where $\eta$ is
the Ohmic diffusion coefficient. It is often treated as equivalent when $\eta$
can be ignored, which is the case of the simulations reported here (Leake et
al., 2005; Leake & Arber, 2006; Arber et al., 2007).

The fraction of neutrals, $\xi_{n}$, is given by
$\xi_{n}=\rho_{n}/(\rho_{c}+\rho_{n})$, where $\rho_{n}$ and $\rho_{c}$ are
the mass densities of neutral and charged components, respectively. The
collisional parameter, $\alpha_{n}$, is given by,

$\alpha_{n}=\rho_{i}\nu_{in}+\rho_{e}\nu_{en}\ ,$ (7)

where $\rho_{i}$ and $\rho_{e}$ are the mass densities of ions and electrons,
and $\nu_{in}$ and $\nu_{en}$ are the collisional frequencies between
ions/neutrals and electrons/neutrals, respectively. The electron, ion, and
neutral mass densities are computed by applying the Saha equation, taking
account of all the relevant chemical species of the solar atmosphere, namely
the first 92 elements from the periodic table with the abundances given in
Anders & Grevesse (1989).

The collisional frequencies are (Braginskii, 1965),

$\nu_{in}~{}=~{}n_{n}\sigma_{in}\sqrt{\dfrac{8k_{B}T}{\pi\mu_{in}}}~{}\propto~{}n_{n}\sqrt{T}\
,$ (8)
$\nu_{en}~{}=~{}n_{n}\sigma_{en}\sqrt{\dfrac{8k_{B}T}{\pi\mu_{en}}}~{}\propto~{}n_{n}\sqrt{T}\
,$ (9)

where the reduced ion/neutral and electron/neutral masses are
$\mu_{in}=m_{i}m_{n}/(m_{i}+m_{n})$ and $\mu_{en}=m_{e}m_{n}/(m_{e}+m_{n})$,
respectively, and the ion/neutral and electron/neutral collision cross
sections are ${\sigma_{in}=5\times 10^{-19}}$ m2 and ${\sigma_{en}=1\times
10^{-19}}$ m2, respectively (Huba, 2013).

##  3 Model and simulation configuration

The sunspot model is produced following the procedures detailed in Khomenko &
Collados (2008). In addition, the model is further modified using the
techniques described by Khomenko et al. (2009) to make sure that it remains
convection-stable. To ensure this, at large radial distances from the sunspot
axis, the model reduces to a quiet Sun model, modified to prevent convective
instability. The standard model ‘S’ from Christensen-Dalsgaard et al. (1996)
is used as the quiet Sun configuration. However, it is adapted to prevent
convective instabilities from developing below the photosphere during the
simulation run. Following the methods documented by Parchevsky & Kosovichev
(2007), the model is adjusted such that the square of the Brunt-Väisälä
frequencies are positive across all atmospheric heights, hence satisfying the
criterion for stability.

To perform 2.5-dimensional (2.5D) wave simulations, our three-dimensional
sunspot model first undergoes a coordinate transformation from cylindrical
coordinates $(r,\phi,z)$ to Cartesian coordinates $(x,y,z)$. Then, a vertical
slice is taken through the axis of the sunspot, along the $y$-axis, such that
the $B_{y}$ component of the magnetic field is equal to zero. Some of the key
characteristics of the model are shown in Figure 1.

Figure 1: Plot showing the key characteristics of the model. Non-
dimensionalized parameters are plotted throughout height at the sunspot axis.
The parameters are, density ($\rho_{0}$), gas pressure ($p_{g,0}$),
temperature ($T_{0}$), and magnetic pressure ($p_{m,0}$). These values have
all been normalized to their photospheric value at the temperature minimum on
the sunspot axis. Their photospheric values are $\rho_{\text{ph}}=4.735\times
10^{-6}$ kg m-3, $p_{g,\text{ph}}=112.4$ N m-2, $T_{\text{ph}}=2857$ K, and
$p_{m,\text{ph}}=2758$ N m-2.

The resulting two-dimensional simulation domain is $14.2$ Mm in the horizontal
direction and $2.8$ Mm in the vertical direction, of which $1$ Mm is above the
photosphere. The sunspot’s axis is located at the center of the simulation
domain. The separation between numerical grid points is $50$ km in the
horizontal direction and $20$ km in the vertical direction. For a 2.5D
simulation, we place a single grid point along the $y$-axis and set the vector
quantities in the MHD equations governing Mancha3D to use all three Cartesian
directions, though in practice, since $B_{y}$ is zero in the simulation plane,
only perturbations in the $x$\- and $z$-components of the velocity and
magnetic field vectors will develop. However, the code also computes
perpendicular currents, which in this case will be along the $y$-direction,
necessitating a 2.5D setup.

The Courant-Friedrichs-Lewy (CFL) numerical stability condition is used by
Mancha3D to determine the time step, which results in the time step in the
simulation with ambipolar diffusion typically varying between $0.001-0.008$
seconds. The dynamic time step used by the simulation is the minimum time from
the advection time step, the diffusion time step, and a fixed maximum time
step provided during initial configuration of the simulation (i.e., $0.008$
seconds). The advection time step will decrease as the maximum absolute value
of the Alfvén speed, the sound speed, and the plasma flow velocity increase.
The diffusion time step will generally decrease as the ambipolar ($\eta_{A}$),
and artificial ($\nu$) diffusion coefficients increase.

Figure 2: Diagram showing the grid layout of the simulation and the location
of key numerical features. The locations along the width of the sunspot model
are given in both pixel coordinates (lower) and physical coordinates (upper).
Similarly, the locations along the height of the sunspot model are given for
pixel coordinates (left) and physical coordinates (right). The vertical dashed
line is the axis of the sunspot, while the horizontal dashed line highlights
the photospheric height, which is taken to be at the temperature minimum along
the sunspot axis. The simulation is driven from the point where these lines
intersect, and this point is shown with a cross. The green shaded region is
the simulation domain which will be analyzed, while the region outside is only
used for numerical purposes. The blue diagonal hatching shows the upper and
lower PMLs, while the orange diagonal hatching shows the left and right
regions simulating periodic boundary conditions. The red and purple sine curve
demonstrates how the model is reflected horizontally to simulate the periodic
boundaries, with each dotted line showing how the region with the same color
of solid line is reflected horizontally.

Periodic boundary conditions are used for the left and right (furthest
distances away from the central axis of the sunspot umbra) boundaries. This is
implemented by reflecting the model horizontally over the right border and
then shifting the model by half the original width to the right such that the
sunspot remains in the center of the $x$-axis. This produces a model with 568
grid points along the horizontal direction, of which 284 grid points are
inside the central simulation domain. The entire numerical domain is large
enough such that waves exiting from the left and right boundaries do not re-
enter the simulation domain or affect its physics. For the upper and lower
boundaries, perfectly matching layers (PMLs; Berenger, 1994) are used to
absorb waves such that open boundaries are simulated. The upper PML is
comprised of 15 additional rows of grid points above the simulation domain
(300 km thick), while the lower PML is comprised of 10 additional rows below
(200 km thick). Inside the upper PML, gravity has been set to zero to assist
with numerical stability. By removing gravity, waves can propagate upwards
more freely, which improves the absorption of perturbations by the PML. To
remain consistent with a zero gravity atmosphere, density ($\rho_{0}$),
temperature ($T_{0}$), and gas pressure ($p_{g,0}$) must all be constants
within the zero gravity region. As both the MHD equations and the model have
been suitably adjusted for a zero gravity atmosphere, the simulation remains
in magnetohydrostatic (MHS) equilibrium, and there are no reflections at the
PML boundary. This strategy has been previously applied by Khomenko & Cally
(2012). The layout of the simulation is shown in detail in Figure 2.

Perturbations are introduced at the photosphere to drive waves through the
atmospheric layers directly above the sunspot. To best represent the effect of
waves in the sunspot atmosphere, we base our wave driver on observationally
acquired Doppler velocity information from the Helioseismic and Magnetic
Imager (HMI; Schou et al., 2012) onboard the Solar Dynamics Observatory (SDO;
Pesnell et al., 2012) spacecraft. The observational data was acquired from
14:36:45 – 17:36:45 UT on 2014 August 30, with the observations reprojected to
account for differential solar rotation using algorithms available in SunPy
(The SunPy Community et al., 2020; Mumford et al., 2022). These data observe
the same sunspot (part of active region NOAA 12149) that was observationally
studied by KP17.

Vector magnetic field information from the HMI/SDO instrument was calculated
using the Very Fast Inversion of the Stokes Vector (VFISV; Borrero et al.,
2011) algorithm, enabling the inclination angles of the sunspot magnetic field
to be extracted. The umbral pixel with the most vertical magnetic field
inclination was identified and a cross-cut through this location was used to
extract the one-dimensional Doppler velocity driver. To ensure the resulting
sunspot wave driver is devoid of long-term trends (e.g., to account for solar
rotation affecting the component of velocity directed along the spacecraft’s
line of sight), a mean quiet Sun velocity time series was extracted, fitted
with a low-order polynomial, and the resulting trend was subtracted from the
one-dimensional sunspot wave driver. A low-frequency sine function was then
subtracted from the driver such that the perturbations oscillate around zero,
ensuring the time series is stationary. The driver was also multiplied by a
constant coefficient of $10^{-3}$ to reduce the maximum perturbation value,
improving numerical stability, resulting in a driver with an RMS force density
of $0.0481$ N m-3. Finally, the driver is spline interpolated to decrease the
cadence of the observed velocities from 45 seconds down to 1 second.

This one-dimensional driver is applied as a perturbation with units of force
density in the MHD equations, as the $z$-component of $\mathbf{S}(t)$, at the
base of the photosphere on the axis of the sunspot, with a Gaussian reduction
to a zero perturbation around this point using a sigma of 3 grid points in
both directions. As the driver has a cadence of 1 second and Mancha3D utilizes
a much shorter time step, within Mancha3D, linear interpolation is applied
between successive driver values. In order to prevent numerical issues
associated with applying a strong driving force to an atmosphere initially at
rest, at the start of the simulation the driving force is multiplied by a
Gaussian, $1-\exp{\left(-t^{2}/(2\sigma^{2})\right)}$, where $\sigma=50$
seconds. This has the effect of gradually increasing the force from zero at
the start of the simulation. The final value of $S_{z}(t)$ is shown in Figure
3.

Figure 3: The observationally obtained driving force,
$\mathbf{S}(t)=S_{z}(t)\hat{\mathbf{z}}$, applied to the sunspot atmosphere
during the simulation around the sunspot’s central axis at photospheric
heights. The HMI/SDO Doppler velocity observations this driver is based upon
have been processed as outlined in Section 3. The horizontal dotted line
represents a zero driving force.

##  4 Results and Discussion

Figure 4: Plots of parameter perturbations during the simulation throughout
time ($x$-axis) and atmospheric height ($z$-axis) along the axis of the
sunspot. The $z$-component of velocity perturbations, $v_{z}$, are shown in
(a) and (b), and the temperature perturbations, $T$, are shown in (c) and (d)
for simulations with (panels a & c) and without (panels b & d) ambipolar
diffusion. Horizontal dotted lines show the location of the driver.

###  4.1 Simulated parameters

Two simulations are evaluated using our initial sunspot model across 3 hours
of solar evolution time, with one simulation modeling ambipolar diffusion,
whilst the other simulation has the ambipolar diffusion terms excluded.
Snapshots of the plasma parameters are saved every 10 seconds of simulation
time, resulting in a total of 1080 snapshots per simulation.

As a tool to monitor the simulation, we can inspect how the plasma parameters
vary throughout time at significant locations in the numerical domain, such as
along the axis of the sunspot. Figure 4 shows the value of the $z$-component
of velocity, $v_{z}$, and the temperature, $T$, for perturbations throughout
time and atmospheric height along the sunspot axis, hence creating a time-
distance diagram for each plasma parameter. In Figure 4, time-distance
diagrams with (panels a & c) and without (panels b & d) ambipolar diffusion
are shown for comparison. These plots also identify the photospheric layer
where the wave-driving perturbations are positioned.

###  4.2 Energy spectral density

In order to investigate the frequency dependence of ambipolar diffusion as a
function of geometric height, we calculate the autospectral density of the
simulated velocity-time series at each atmospheric height output by the
Mancha3D code. As can been seen in Figure 4, the output velocity data,
$v_{z}$, is non-stationary, i.e., its statistical properties vary throughout
time (Bendat & Piersol, 2010). Therefore, we must calculate energy spectra
rather than power spectra.

The energy autospectral density function, $\mathscr{G}_{xx}(f)$, is related to
the power autospectral density function, $G_{xx}(f)$, as,

$\mathscr{G}_{xx}(f)=LG_{xx}(f)\ ,$ (10)

where $L$ is the length of data available and $G_{xx}(f)$ is the one-sided
autospectral density function given by,

$G_{xx}(f)=2\lim_{L\rightarrow\infty}\dfrac{1}{L}E\Big{[}|X(f)|^{2}\Big{]}\ ,$
(11)

where $X(f)$ is the finite Fourier transform of the velocity-time series with
a sampling interval of $L$ seconds and $E[\dots]$ represents the expected
value (Bendat & Piersol, 2010). As we are working with a discrete time series,
we apply the Fast Fourier Transform (FFT) algorithm and also use an estimate
of $G_{xx}(f)$,

$\tilde{G}_{xx}(F)=\dfrac{2}{L}|X(f)|^{2}\ .$ (12)

This results in the one-sided energy autospectral density, which we call the
energy spectrum, namely (Stull, 1988),

$\mathscr{G}_{xx}(f)=2|X(f)|^{2}\ .$ (13)

We apply this equation independently to the time series extracted at each
pixel along the sunspot’s height axis ($z$-axis in the numerical domain). This
allows for the frequency response of waves to be compared across atmospheric
heights. When applying the FFT algorithm to the $z$-component of velocity,
$v_{z}(t)$, we exclude the initial 1000 seconds of the simulation in order to
give the driven oscillations time to strengthen. The maximum frequency
analyzed is limited by the 45 second cadence of the underlying HMI/SDO
observations, providing us with a Nyquist frequency of $\approx 11.1$ mHz. To
aid with the comparison to other studies, we divide all energy spectra by the
corresponding frequency sampling, $\Delta
f=1/L=1/(10\,800~{}\textrm{s}-1000~{}\textrm{s})\approx 0.102~{}\textrm{mHz}$,
to provide an energy spectral density (i.e., in units of km2 s-2 mHz-1) that
is suitable for scientific examination.

Furthermore, by calculating the energy spectral densities for both simulations
independently, this allows us to infer how wave behavior changes with
atmospheric height, and hence compare how the presence of ambipolar diffusion
affects these characteristics. Figures 5 & 6 show examples of energy spectral
densities plotted for a number of key atmospheric heights across both
simulations. At the location of the driver the energy spectral densities
between both simulations have similar values, while at the top of the
simulation domain the energy spectral densities exhibit more distinct
differences between the two simulations.

Figure 5: Plots of the energy spectral density at different heights within the
solar atmosphere along the $z$-axis of the sunspot for the simulation with
ambipolar diffusion terms included in the MHD equations. The top of the
simulation domain ($\approx 1$ Mm) is shown in (a), while the location of the
driver ($\approx 0$ Mm) is shown in (b). The vertical dashed line is at
$\approx 6.1$ mHz, after which the spectral points are marked with crosses
instead of dots. A power-law fit is applied to the points with frequencies
$>6.1$ mHz, where the fit is shown with a solid red line. Figure 6: Identical
plots to those shown in Figure 5, only now with no ambipolar diffusion
included in the relevant MHD equations.

###  4.3 Power-law indices

As highlighted in the work of KP17, observations of a sunspot umbra acquired
through different filters (corresponding to different heights in the solar
atmosphere) produced varying power-law indices for frequencies above the peak
value. In the study documented by KP17, the frequency that displayed dominant
spectral energy was $\approx 6.5$ mHz (corresponding to a period of 152
seconds). This frequency is consistent with the typical upper-
photospheric/chromospheric $p$-mode spectrum, where only frequencies greater
than the acoustic cutoff frequency (i.e., $>5$ mHz) are permitted to propagate
(see, e.g., Kanoh et al., 2016; Murawski et al., 2016; Felipe et al., 2018).
For frequencies higher than this dominant value, KP17 found power-law slopes
of the energy spectral densities to vary depending on the layer of the
atmosphere being sampled. In our present work, Figures 5 & 6 reveal a spectral
peak at $\approx 6.1$ mHz (periodicity of 164 seconds), which is consistent
with the observational work of KP17. At higher frequencies (i.e., $>6.1$ mHz),
we find the energy spectral densities decrease progressively up to the Nyquist
frequency of $\approx 11.1$ mHz.

By measuring the slopes of the energy spectral densities beyond the $\approx
6.1$ mHz enhancement, we are able to calculate a power-law index for this
frequency region at each atmospheric height, for both simulations, in order to
compare the effects of ambipolar diffusion on the resulting power-law slopes.
The energy spectral densities with a frequency beyond the enhancement are
selected, with a logarithmic function applied to the frequencies and energy
spectral densities such that the linear relationship can be determined. An
unweighted least squares first order polynomial fit is applied to the log-log
frequencies and energy spectral densities, with the power-law slope taken to
be the slope of the fitted first order polynomial. Examples of the fitted
power-law slopes are plotted for each simulation, at specific atmospheric
heights, in Figures 5 & 6, with the corresponding power-law indices displayed
as a function of atmospheric height in Figure 7. The work of KP17 focused on
observations spanning the mid-photosphere (e.g., data acquired in the Mg i b2
spectral line) through to the mid-chromosphere (e.g., H$\alpha$ and Ca ii K
observations). As a result, we focus our analysis between the atmospheric
heights of $500-1000$ km in order to examine the layers associated with the
photospheric/chromospheric interface.

In Figure 7 it can be seen that the power-law index for the simulation without
ambipolar diffusion increases monotonically with atmospheric height throughout
the $500-1000$ km range. This is not the case when ambipolar diffusion is
accounted for in the simulation, resulting in the power-law index decreasing
with atmospheric height between the heights of $680-840$ km in the atmosphere,
before increasing again above this height, as seen in Figure 7.

###  4.4 Ambipolar diffusion coefficient

To understand the differences between the power-law slopes arising in the
simulations with/without ambipolar diffusion, we can analyze and decompose the
ambipolar diffusion coefficient, $\eta_{A}$, in the MHD equations. As
discussed in Section 2, the coefficient is given by
$\eta_{A}=\xi^{2}_{n}|\mathbf{B}|^{2}/\alpha_{n}$. The magnetic field,
$|\mathbf{B}|$, and the temperature, $T$, do not vary significantly between
the range of atmospheric heights spanning $500-1000$ km, so their effects on
the power-law index differences within this range can be neglected.

In particular, the magnetic field squared, $|\mathbf{B}|^{2}$, within the
simulation with ambipolar diffusion varies approximately linearly from
$0.0062$ T2 at 500 km to $0.0056$ T2 at 1000 km, and has a standard deviation
of $0.00018$ T2. As a relatively constant value is maintained, the mean
magnetic field, $\overline{|\mathbf{B}|^{2}}=0.0059$ T2, is used in the
analysis. This approximation is applied by multiplying the ambipolar diffusion
coefficient, $\eta_{A}$, by $\overline{|\mathbf{B}|^{2}}/|\mathbf{B}|^{2}$,
which can be seen in Figure 8 to be relatively close to unity.

In Section 2, we introduced the neutral collisional parameter,
$\alpha_{n}=\rho_{i}\nu_{in}+\rho_{e}\nu_{en}$. We also showed that the
collisional frequencies are proportional to $n_{n}\sqrt{T}$. Therefore, the
neutral collisional parameter can be re-written as $\alpha_{n}\propto
n_{e}n_{n}\sqrt{T}$. Furthermore, the ambipolar diffusion term can be examined
analytically, where it can be shown that,

$\displaystyle\eta_{A}$ $\displaystyle\propto\xi_{n}^{2}\dfrac{1}{\alpha_{n}}$
(14)
$\displaystyle\propto\dfrac{n_{n}^{2}}{(n_{e}+n_{n})^{2}}\dfrac{1}{n_{e}n_{n}\sqrt{T}}$
(15)
$\displaystyle\propto\dfrac{n_{n}}{n_{e}(n_{e}+n_{n})^{2}}\dfrac{1}{\sqrt{T}}\
.$ (16)

Figure 7: Plot of the power-law indices of the energy spectral densities
within the frequency interval of $6.1-11.1$ mHz, calculated at different
heights within the solar atmosphere along the $z$-axis of the sunspot. The
black line is for the simulation without ambipolar diffusion (AD), while the
red line is for the simulation that contains AD terms in the associated MHD
equations.

In the case where there is a high degree of plasma ionization (i.e., $n_{e}\gg
n_{n}$), the ambipolar diffusion coefficient tends towards zero as expected.
However, for a partially ionized plasma with $n_{n}\approx n_{e}$ (i.e.,
$\xi_{n}\approx\frac{1}{2}$) or $n_{n}\gg n_{e}$ (i.e., $\xi_{n}\approx 1$),
the equation follows that,

$\displaystyle\eta_{A}$ $\displaystyle\propto\dfrac{1}{\alpha_{n}}$ (17)
$\displaystyle\propto\dfrac{1}{n_{e}n_{n}\sqrt{T}}\ .$ (18)

For the partially ionized plasma in our sunspot model, this results in a
proportionality between the ambipolar diffusion coefficient, $\eta_{A}$, and
collisional parameter, $\alpha_{n}$. Therefore, we would expect ambipolar
diffusion to have the most significant effect on partially ionized plasmas
when the collisional parameter is particularly small. In the case where the
temperature does not vary significantly, this corresponds to plasma where the
product of the number densities is small, which is the case in the ambipolar
diffusion simulation presented in this study. Relative to the variation in the
product of the number densities, temperature is approximately constant over
the height range of $500-1000$ km and the relationship simplifies to
$\alpha_{n}\propto n_{e}n_{n}$. In particular, the square root of the
temperature varies from a minimum of 67 K1/2 to a maximum of 86 K1/2, with a
mean of 78 K1/2 and a standard deviation of 3 K1/2, which is an insignificant
variation compared to variations of the number densities. The product of the
number densities, $n_{e}n_{n}$, varies from a minimum at $10^{33}$ to a
maximum at $10^{37}$, with a mean value of $2\times 10^{35}$ and a standard
deviation of $4\times 10^{35}$. The mean temperature is therefore used in the
analysis, with the ambipolar diffusion coefficient, $\eta_{A}$, adjusted by
multiplying it by $\sqrt{T}/\overline{\sqrt{T}}$. Figure 8 shows that the
effect of this adjustment on $\eta_{A}$ is small. This final adjustment leads
to,

$\eta_{A}\propto\dfrac{1}{n_{e}n_{n}}\ .$ (19)

The combined effect of all these approximations is shown in Figure 9 where the
ambipolar diffusion coefficient, $\eta_{A}$, is plotted alongside its
approximation,

$\eta_{A}^{\prime}=\eta_{A}\dfrac{\overline{|\mathbf{B}|^{2}}}{|\mathbf{B}|^{2}}\dfrac{1}{\xi_{n}^{2}}\dfrac{\sqrt{T}}{\overline{\sqrt{T}}}\
.$ (20)

Figure 8: Plots of the individual approximations applied to the ambipolar
diffusion coefficient, $\eta_{A}$, to produce $\eta_{A}^{\prime}$ as shown in
Equation 20. The approximation $1/\xi_{n}^{2}$ is plotted using a dashed black
line, $\overline{|\mathbf{B}|^{2}}/|\mathbf{B}|^{2}$ is plotted using a solid
blue line, and $\sqrt{T}/\overline{\sqrt{T}}$ is plotted with a dash-dotted
orange line. The solid black horizontal line highlights the value of these
approximations in relation to unity. All of the terms are extracted from the
Mancha3D simulation with ambipolar diffusion along the $z$-axis of the
sunspot. The mean of the terms are then taken throughout time to get their
average value at each atmospheric height. The terms with an overline use the
mean value calculated over all dimensions. Figure 9: Plots of the ambipolar
diffusion coefficient, $\eta_{A}$, along with its approximation,
$\eta_{A}^{\prime}$. The coefficient $\eta_{A}$ is plotted using a solid black
line against the left axis, while $\eta_{A}^{\prime}$ is plotted using a
dashed red line against the left axis. The difference between both lines,
$\eta_{A}^{\prime}-\eta_{A}$, is plotted using a solid blue line against the
right axis. All of the terms are extracted from the Mancha3D simulation with
ambipolar diffusion along the $z$-axis of the sunspot. The mean of the terms
are then taken throughout time to get their average value at each atmospheric
height.

Comparing these analytical results to outputs provided by the Mancha3D
simulation, we find close agreement with our theoretical considerations above.
Figures 8 & 9 highlight that it is the collisional parameter, $\alpha_{n}$,
rather than the neutral fraction, $\xi_{n}$, that is significant in
determining the value of the ambipolar diffusion coefficient, $\eta_{A}$, as
analytically shown in Equation 17. Furthermore, comparing Figure 9 to Figure
7, we are able to see that the reciprocal of the product of the number
densities, i.e., $\eta_{A}\propto 1/(n_{e}n_{n})$, becomes significantly large
only within the range of atmospheric heights where the greatest difference in
power-law index is observed between the simulations with and without ambipolar
diffusion. Therefore, this implies that the ambipolar diffusion coefficient,
$\eta_{A}$, which is proportional to the reciprocal of the collisional
parameter, $\alpha_{n}$, is responsible for the steepening of the power-law
index between the atmospheric heights of $680-840$ km.

The steeper downward power-law slopes suggest that ambipolar diffusion leads
to increased damping of higher frequency waves. Khomenko et al. (2018) show
this effect occurring in simulations of magneto-convection within the lower
solar atmosphere. The authors find that when ambipolar diffusion is modeled,
the power at higher frequencies is diminished to a greater extent than their
lower frequency counterparts. The observational work of KP17 also suggests
that ambipolar diffusion may enhance frequency-dependent damping, which is now
readily visualized in Figure 7.

In Figure 10, we plot the reciprocals of the terms $n_{e}$ and $n_{n}$
separately. On average, the reciprocal of the electron number density
($1/n_{e}$) has a bell-shaped curve with a peak at 700 km, the height where
there are fewest free electrons. On the other hand, the reciprocal of the
neutral number density ($1/n_{n}$) is small initially, but increases gradually
as atmospheric height increases, meaning that, on average, there are fewer
neutrals higher in the atmosphere. It is these two terms multiplied together
that results in the neutral collisional term, $1/\alpha_{n}$, exhibiting the
bell-shaped curve that is shown in Figure 9, which, on average, has a peak
value at 760 km in the atmosphere. As a result, we are able to conclude that
the maximal interplay between the reciprocals of the $n_{e}$ and $n_{n}$ terms
(maximizing the atmospheric distributions when there are fewest free electrons
and the smallest number of neutrals — i.e., corresponding to the least dense
plasma) provides the strongest effect on the wave power-law indices. These
effects are clearly visible in Figures 7, 9 & 10, where strong deviations in
the corresponding plots are evident at an atmospheric height of $\sim 760$ km.

Figure 10: Plots of the reciprocals of the number densities for electrons
($1/n_{e}$) and neutrals ($1/n_{n}$). The electron number density reciprocal
is represented by a solid black line against the left axis, whilst the the
neutral number density reciprocals plotted using a dashed red line against the
right axis. The reciprocals are calculated at every point along the sunspot
$z$-axis and the time-average value is calculated.

##  5 Conclusions

Magnetohydrodynamic wave simulations were evaluated in 2.5D for a cross-
section of the chromosphere above a sunspot umbra, with the energy spectral
densities calculated at each atmospheric height. Two simulations were produced
independently, one with ambipolar diffusion terms included in the embedded MHD
equations and one without. We compare the energy spectral densities between
the two simulations, focusing our analysis within the frequency range spanning
$6.1-11.1$ mHz, which can be fitted with a power-law relationship following
the observational work of KP17. We focus our analysis on the atmospheric
heights spanning $500-1000$ km, where previous wave studies of the mid-
photosphere to mid-chromosphere have been performed on the same sunspot
structure (KP17).

From the simulations, we provided evidence that increases in the ambipolar
diffusion coefficient, $\eta_{A}$, was a result of the collisional parameter,
$\alpha_{n}$, becoming smaller. We showed that in our sunspot model the
collisional parameter is proportional to $n_{e}n_{n}$, where $n_{e}$ and
$n_{n}$ are the number densities of the electrons and neutrals, respectively.
Using this relationship, we are now able to state that the process of
ambipolar diffusion in a sunspot atmosphere is significant, and needs to be
accounted for, when the value of $1/(n_{e}n_{n})$ becomes relatively large. We
find the value of the ambipolar diffusion coefficient ($\eta_{A}\propto
1/n_{e}n_{n}$) to be closely linked with variations seen between the two
simulations presented in this study. At atmospheric heights where there was an
insignificant difference between the energy spectral density power-law slopes
of both simulations, the ambipolar diffusion coefficient remained relatively
low. This resulted in the terms modeling ambipolar diffusion to have an
insignificant effect on the plasma parameters. However, at atmospheric heights
where significant variations were present in the power-law slopes of the
energy spectral densities, the ambipolar diffusion coefficient was found to be
particularly large. Under this regime (approximately $680-920$ km above the
solar surface), ambipolar diffusion caused the power-law slopes in the range
of $6.1-11.1$ mHz to become steeper, hinting at increased wave energy damping
in this environment. Interestingly, we do not find that the ionization
fraction, $\xi_{e}$, has a significant role in determining where ambipolar
diffusion is important in the sunspot atmosphere, even though it is part of
the ambipolar diffusion coefficient, $\eta_{A}$.

The importance of the product of the electron number density, $n_{e}$, and the
neutral number density, $n_{n}$, highlights that ambipolar diffusion is more
significant in less dense plasma. While in general we can expect ambipolar
diffusion to be more important higher in the atmosphere where the plasma is
less dense, we also need to consider the relative ionization of each of the
plasma components.

This work provides evidence of ambipolar diffusion having a significant effect
on the plasma properties within an umbral sunspot atmosphere and provides a
number of insights for future sunspot studies. The inclusion of ambipolar
diffusion terms in MHD simulations of sunspot atmospheres should be fully
considered, particularly when the ambipolar diffusion coefficient ($\eta_{A}$)
approximation ($1/n_{e}n_{n}$) is large. We also suggest that ambipolar
diffusion may produce a clear effect on the power-law indices of energy
spectral densities that is detectable in observations of the lower solar
atmosphere of sunspots.

The next stage of this work will be to acquire spectral observations (to
obtain Doppler velocity measurements instead of the intensity signatures
employed by KP17, ) of a sunspot across a range of wavelengths and analyze how
the energy spectral density power-law indices vary over atmospheric height. We
would expect to observe a steepening of the power-law slopes at similar
atmospheric heights to our simulation when ambipolar diffusion was included
($680-840$ km as shown in Figure 7). This would provide direct observational
evidence of ambipolar diffusion occurring within the lower chromosphere above
a sunspot.

Next generation spectral imaging instruments, such as the Visible Tunable
Filter (VTF; Schmidt et al., 2016) and the Diffraction Limited Near Infrared
Spectropolarimeter (DL-NIRSP), which will be available on the Daniel K. Inouye
Solar Telescope (DKIST; Tritschler et al., 2016)), will be essential for this
work. High-precision spectral observations over a fine spatial and temporal
grid are vital to calculate robust power-law slopes of energy spectral
densities to compare to the slopes from our simulations. In addition to the
observational instruments, advanced software tools will be required to extract
precise Doppler velocity measurements from the spectral images (MacBride et
al., 2021), which can be combined with cutting-edge inversion routines (e.g.,
Beck et al., 2019; de la Cruz Rodríguez et al., 2019; Ruiz Cobo et al., 2022)
to provide better height resolution of where the Doppler signatures are
formed, hence improving the reliability of comparisons between observations
and theory.

C.D.M. would like to thank the Northern Ireland Department for the Economy for
the award of a PhD studentship. D.B.J. is grateful to Invest NI and Randox
Laboratories Ltd. for the award of a Research & Development Grant (059RDEN-1).
D.B.J. also wishes to thank the UK Space Agency for a National Space
Technology Programme (NSTP) Technology for Space Science award (SSc 009).
D.B.J and S.D.T.G thanks the UK Science and Technology Facilities Council
(STFC) for the consolidated grant ST/T00021X/1. E.K. thanks the support by the
European Research Council through the Consolidator Grant
ERC-2017-CoG-771310-PI2FA and by the Spanish Ministry of Economy and the
Industry and Competitiveness through the grant PGC2018-095832-B-I00. Data
supporting this study are available upon request from the corresponding
author.

## References

  * Anders & Grevesse (1989) Anders, E., & Grevesse, N. 1989, Geochimica et Cosmochimica Acta, 53, 197, doi: 10.1016/0016-7037(89)90286-X
  * Arber et al. (2007) Arber, T. D., Haynes, M., & Leake, J. E. 2007, The Astrophysical Journal, 666, 541, doi: 10.1086/520046
  * Ballester et al. (2020) Ballester, J. L., Soler, R., Terradas, J., & Carbonell, M. 2020, Astronomy & Astrophysics, 641, A48, doi: 10.1051/0004-6361/202038220
  * Ballester et al. (2018) Ballester, J. L., Alexeev, I., Collados, M., et al. 2018, Space Science Reviews, 214, 58, doi: 10.1007/s11214-018-0485-6
  * Beck et al. (2019) Beck, C., Gosain, S., & Kiessner, C. 2019, The Astrophysical Journal, 878, 60, doi: 10.3847/1538-4357/ab1d4c
  * Bendat & Piersol (2010) Bendat, J. S., & Piersol, A. G. 2010, Random Data: Analysis and Measurement Procedures, 4th edn., Wiley Series in Probability and Statistics (Wiley)
  * Berenger (1994) Berenger, J.-P. 1994, Journal of Computational Physics, 114, 185, doi: 10.1006/jcph.1994.1159
  * Biermann (1950) Biermann, L. 1950, Zeitschrift Naturforschung Teil A, 5, 65
  * Borrero et al. (2011) Borrero, J. M., Tomczyk, S., Kubo, M., et al. 2011, Solar Physics, 273, 267, doi: 10.1007/s11207-010-9515-6
  * Braginskii (1965) Braginskii, S. I. 1965, Reviews of Plasma Physics, 1, 205
  * Christensen-Dalsgaard et al. (1996) Christensen-Dalsgaard, J., Däppen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286, doi: 10.1126/science.272.5266.1286
  * Cowling (1957) Cowling, T. G. 1957, Interscience Tracts on Physics and Astronomy, Vol. 4, Magnetohydrodynamics, 3rd edn. (Interscience Publishers)
  * de la Cruz Rodríguez et al. (2019) de la Cruz Rodríguez, J., Leenaarts, J., Danilovic, S., & Uitenbroek, H. 2019, Astronomy & Astrophysics, 623, A74, doi: 10.1051/0004-6361/201834464
  * De Pontieu et al. (2014) De Pontieu, B., Title, A. M., Lemen, J. R., et al. 2014, Solar Physics, 289, 2733, doi: 10.1007/s11207-014-0485-y
  * Díaz et al. (2012) Díaz, A. J., Soler, R., & Ballester, J. L. 2012, The Astrophysical Journal, 754, 41, doi: 10.1088/0004-637X/754/1/41
  * Felipe et al. (2010) Felipe, T., Khomenko, E., & Collados, M. 2010, The Astrophysical Journal, 719, 357, doi: 10.1088/0004-637X/719/1/357
  * Felipe et al. (2018) Felipe, T., Kuckein, C., & Thaler, I. 2018, Astronomy & Astrophysics, 617, A39, doi: 10.1051/0004-6361/201833155
  * Huba (2013) Huba, J. D. 2013, NRL Plasma Formulary (Naval Research Laboratory)
  * Jess et al. (2012) Jess, D. B., De Moortel, I., Mathioudakis, M., et al. 2012, The Astrophysical Journal, 757, 160, doi: 10.1088/0004-637X/757/2/160
  * Jess et al. (2010) Jess, D. B., Mathioudakis, M., Christian, D. J., et al. 2010, Solar Physics, 261, 363, doi: 10.1007/s11207-009-9500-0
  * Jess et al. (2015) Jess, D. B., Morton, R. J., Verth, G., et al. 2015, Space Science Reviews, 190, 103, doi: 10.1007/s11214-015-0141-3
  * Kanoh et al. (2016) Kanoh, R., Shimizu, T., & Imada, S. 2016, The Astrophysical Journal, 831, 24, doi: 10.3847/0004-637X/831/1/24
  * Khomenko & Cally (2012) Khomenko, E., & Cally, P. S. 2012, The Astrophysical Journal, 746, 68, doi: 10.1088/0004-637X/746/1/68
  * Khomenko & Collados (2006) Khomenko, E., & Collados, M. 2006, The Astrophysical Journal, 653, 739, doi: 10.1086/507760
  * Khomenko & Collados (2008) —. 2008, The Astrophysical Journal, 689, 1379, doi: 10.1086/592681
  * Khomenko & Collados (2012) —. 2012, The Astrophysical Journal, 747, 87, doi: 10.1088/0004-637X/747/2/87
  * Khomenko et al. (2014a) Khomenko, E., Collados, M., Díaz, A., & Vitas, N. 2014a, Physics of Plasmas, 21, 092901, doi: 10.1063/1.4894106
  * Khomenko et al. (2021) Khomenko, E., Collados, M., Vitas, N., & González-Morales, P. A. 2021, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379, 20200176, doi: 10.1098/rsta.2020.0176
  * Khomenko et al. (2014b) Khomenko, E., Díaz, A., de Vicente, A., Collados, M., & Luna, M. 2014b, Astronomy & Astrophysics, 565, A45, doi: 10.1051/0004-6361/201322918
  * Khomenko et al. (2009) Khomenko, E., Kosovichev, A., Collados, M., Parchevsky, K., & Olshevsky, V. 2009, The Astrophysical Journal, 694, 411, doi: 10.1088/0004-637X/694/1/411
  * Khomenko et al. (2018) Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2018, Astronomy & Astrophysics, 618, A87, doi: 10.1051/0004-6361/201833048
  * Krishna Prasad et al. (2017) Krishna Prasad, S., Jess, D. B., Van Doorsselaere, T., et al. 2017, The Astrophysical Journal, 847, 5, doi: 10.3847/1538-4357/aa86b5
  * Leake & Arber (2006) Leake, J. E., & Arber, T. D. 2006, Astronomy & Astrophysics, 450, 805, doi: 10.1051/0004-6361:20054099
  * Leake et al. (2005) Leake, J. E., Arber, T. D., & Khodachenko, M. L. 2005, Astronomy & Astrophysics, 442, 1091, doi: 10.1051/0004-6361:20053427
  * MacBride et al. (2021) MacBride, C. D., Jess, D. B., Grant, S. D. T., et al. 2021, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379, 20200171, doi: 10.1098/rsta.2020.0171
  * Martínez-Gómez et al. (2021) Martínez-Gómez, D., Popescu Braileanu, B., Khomenko, E., & Hunana, P. 2021, Astronomy & Astrophysics, 650, A123, doi: 10.1051/0004-6361/202039113
  * Mumford et al. (2022) Mumford, S. J., Freij, N., Stansby, D., et al. 2022, SunPy, Zenodo, doi: 10.5281/ZENODO.591887
  * Murawski et al. (2016) Murawski, K., Musielak, Z. E., Konkol, P., & Wiśniewska, A. 2016, The Astrophysical Journal, 827, 37, doi: 10.3847/0004-637X/827/1/37
  * Ni et al. (2015) Ni, L., Kliem, B., Lin, J., & Wu, N. 2015, The Astrophysical Journal, 799, 79, doi: 10.1088/0004-637X/799/1/79
  * Pandey & Wardle (2008) Pandey, B. P., & Wardle, M. 2008, Monthly Notices of the Royal Astronomical Society, 385, 2269, doi: 10.1111/j.1365-2966.2008.12998.x
  * Parchevsky & Kosovichev (2007) Parchevsky, K. V., & Kosovichev, A. G. 2007, The Astrophysical Journal, 666, 547, doi: 10.1086/520108
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Solar Physics, 275, 3, doi: 10.1007/s11207-011-9841-3
  * Popescu Braileanu & Keppens (2021) Popescu Braileanu, B., & Keppens, R. 2021, Astronomy & Astrophysics, 653, A131, doi: 10.1051/0004-6361/202140872
  * Popescu Braileanu et al. (2021a) Popescu Braileanu, B., Lukin, V. S., Khomenko, E., & de Vicente, Á. 2021a, Astronomy & Astrophysics, 646, A93, doi: 10.1051/0004-6361/202039053
  * Popescu Braileanu et al. (2021b) —. 2021b, Astronomy & Astrophysics, 650, A181, doi: 10.1051/0004-6361/202140425
  * Ruiz Cobo et al. (2022) Ruiz Cobo, B., Quintero Noda, C., Gafeira, R., et al. 2022, Astronomy & Astrophysics, 660, A37, doi: 10.1051/0004-6361/202140877
  * Saha (1921) Saha, M. N. 1921, Proceedings of the Royal Society of London Series A, 99, 135, doi: 10.1098/rspa.1921.0029
  * Schmidt et al. (2016) Schmidt, W., Schubert, M., Ellwarth, M., et al. 2016, in Ground-Based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, Vol. 9908 (Edinburgh, United Kingdom: SPIE), 99084N, doi: 10.1117/12.2232518
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Solar Physics, 275, 229, doi: 10.1007/s11207-011-9842-2
  * Shelyag et al. (2016) Shelyag, S., Khomenko, E., de Vicente, A., & Przybylski, D. 2016, The Astrophysical Journal, 819, L11, doi: 10.3847/2041-8205/819/1/L11
  * Soler et al. (2012) Soler, R., Díaz, A. J., Ballester, J. L., & Goossens, M. 2012, The Astrophysical Journal, 749, 163, doi: 10.1088/0004-637X/749/2/163
  * Stull (1988) Stull, R. B. 1988, An Introduction to Boundary Layer Meteorology, Atmospheric Sciences Library (Kluwer Academic Publishers)
  * The SunPy Community et al. (2020) The SunPy Community, Barnes, W. T., Bobra, M. G., et al. 2020, The Astrophysical Journal, 890, 68, doi: 10.3847/1538-4357/ab4f7a
  * Tritschler et al. (2016) Tritschler, A., Rimmele, T. R., Berukoff, S., et al. 2016, Astronomische Nachrichten, 337, 1064, doi: 10.1002/asna.201612434
  * Vernazza et al. (1981) Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, The Astrophysical Journal Supplement Series, 45, 635, doi: 10.1086/190731


# On activation in solid ionic electrolytes

K. Noori  NUS Institute for Functional Intelligent Materials, 4 Science Drive
2, Singapore 117544, Singapore  B. A. Olsen  Department of Physics, Lewis &
Clark College, Portland, Oregon, 97219, USA  A. Rodin  Yale-NUS College, 16
College Avenue West, 138527, Singapore  Centre for Advanced 2D Materials,
National University of Singapore, 117546, Singapore  Department of Materials
Science and Engineering, National University of Singapore, 117575, Singapore

###### Abstract

Ionic conductivity in solid electrolytes is commonly expected to exhibit
Arrhenius dependence on temperature, determined by a well-defined activation
energy. Consequently, a standard approach involves calculating this energy
using quasi-static methods and using the Arrhenius form to extrapolate the
numerical results from one temperature range to another. Despite the ubiquity
of this Arrhenius-based modeling, disagreements frequently arise between
theory and experiment, and even between different theoretical studies. By
considering a tractable minimal model, we elucidate the reason behind the
breakdown of the Arrhenius conductivity form. This breakdown is driven by non-
trivial phase-space boundaries between conducting and non-conducting regimes,
and depends on the kinetic properties of the system.

##  I Introduction

Ionic motion through solid electrolytes [1, 2] plays a central role in the
operation of advanced batteries [3, 4, 5]. One of the most important
quantities associated with the performance of these electrolytes is their
conductivity, generally accepted to take the form of
$\sigma=\sigma_{0}e^{-E_{a}/k_{B}T}/T$, [1, 3, 4, 5] where $T$ is the
temperature, $k_{B}$ is the Boltzmann constant, and $E_{a}$ is referred to as
the activation energy. In this formulation, ionic conductivity in solids
increases with temperature; the exponential term determines the number of
charge carriers, while $T^{-1}$ arises from the temperature dependence of
mobility. This dependence of $\sigma$ on $T$ is confirmed experimentally, with
the most explicit demonstration being a linear relation between $\log(\sigma
T)$ and $1/T$ [6, 7, 8, 9, 10], and others finding a linear dependence of
$\log(\sigma)$ on $1/T$ [11, 12, 13, 14, 15].

When current-carrying ions move through the solid electrolyte, they navigate a
nontrivial potential landscape produced by the stationary ions of the solid.
Traversing this landscape requires overcoming local energy maxima, which is
achievable only for sufficiently energetic mobile ions. Assuming that the
ionic energies are Boltzmann-distributed and the energy barrier is given by
$E_{a}$, the number of ions in the system that can overcome the barrier at any
given time is $\propto e^{-E_{a}/k_{B}T}$, giving rise to the Arrhenius
dependence of $\sigma$. This statement can be rephrased if the ionic motion is
assumed to be ergodic: the probability that a given ion manages to overcome
the energy barrier is $\propto e^{-E_{a}/k_{B}T}$.

This fact has implications for the interpretation of numerical studies on
conductivity. For example, molecular dynamics (MD) simulations—especially in
their ab initio formulation (AIMD)—are often used to estimate $\sigma$. Owing
to the high computational cost of AIMD, simulations are typically restricted
to a small set of temperatures; conductivity values at other temperatures are
then extrapolated using the Arrhenius relation. Another common approach is to
use nudged elastic band (NEB) simulations to compute the migration energy
barrier, $E_{\mathrm{NEB}}$, along the presumed ionic migration pathway. This
value is assumed to be equivalent to $E_{a}$ and thus serves as a proxy for
$\sigma$. However, there is experimental evidence of a breakdown in the linear
relation between $\log(\sigma)$ and $1/T$ in common solid electrolyte
materials [16]. Moreover, recent AIMD studies with denser temperature sampling
over extended simulation times using machine-learning assisted potentials have
also shown a departure from Arrhenius scaling [17].

To understand this deviation, note that the Arrhenius model assumes that the
energy barrier is independent of temperature. While this assumption is
generally true for electronic semiconductors, where the gap does not vary with
temperature, thermal vibrations of the lattice in ionic conductors lead to
significant fluctuations of the potential energy surface. These fluctuations
depend on temperature and preclude the existence of a well-defined and
temperature-independent activation energy.

In this work, we develop a minimal model to demonstrate how even the simplest
ionic systems deviate from Arrhenius scaling. In Sec. II, we propose a
theoretical explanation of why the Arrhenius behavior will generally not arise
in ionic conductors over extended temperature ranges. To substantiate this
model, in Sec. III, we consider a numerically-tractable system and demonstrate
the breakdown of the Arrhenius scaling in the presence of thermal
fluctuations. We conclude with a discussion of our results and their
implications for state-of-the-art numerical techniques in Sec. IV.

We perform all our calculations using JULIA [18] and make our code available
at https://github.com/rodin-physics/oscillator-activation. Our plots are
visualized using Makie.jl, [19] employing a color scheme suitable for color-
blind readers, developed in Ref. [20].

##  II Theory

Consider a large system with a single mobile ion. At a particular time, let
the velocity and position of the ion be given by
$\left(\mathbf{V},\mathbf{R}\right)$, six phase space coordinates in total.
Additionally, let the velocities and positions of all the other ions in the
system be described by a $6N$-dimensional vector $\bm{\xi}$, where $N$ is the
number of ions in the solid. For any $\mathbf{R}$, the mobile ion lies within
a polyhedron whose vertices correspond to the solid’s atoms (see Fig. 1 for a
2D analog). As the system evolves in time, the ion will either cross the face
of the polyhedron toward which it is headed or be reflected. Thus, in the
phase space of $(\mathbf{V},\mathbf{R},\bm{\xi})$, there are two groups of
points: the ones which result in the ion traversing the polyhedron face and
the ones that result in its reflection. The escape probability is then given
by the number of $(\mathbf{V},\mathbf{R},\bm{\xi})$ escape configurations,
weighted by their probability $\propto e^{-\beta
E(\mathbf{V},\mathbf{R},\bm{\xi})}$, where $E(\mathbf{V},\mathbf{R},\bm{\xi})$
is the energy of the configuration and $\beta=1/k_{B}T$. We underscore that
the energy of the entire system follows the Boltzmann distribution, as
expected for a macroscopic system at some temperature $T$. Defining the set of
the escape configurations by $\mathcal{A}$ (for activation), we get the
following expression for the escape probability

$\mathcal{P}=\left[\sum_{\mathbf{V},\mathbf{R},\bm{\xi}}e^{-\beta
E(\mathbf{V},\mathbf{R},\bm{\xi})}\right]^{-1}\sum_{(\mathbf{V},\mathbf{R},\bm{\xi})\in\mathcal{A}}e^{-\beta
E(\mathbf{V},\mathbf{R},\bm{\xi})}\,.$ (1)

From Eq. (1), it is not immediately obvious that $\mathcal{P}\propto e^{-\beta
E_{a}}$, as would be the case if the conductivity had the Arrhenius form.
However, for $N\rightarrow\infty$, it is difficult to provide a definitive
proof one way or another. Therefore, to make progress, we make a major
simplification.

Figure 1: _Mobile particle moving through a lattice._ Lattice ions
sufficiently close to the mobile particle, included in the coordinate subset
$\bm{\zeta}$, directly affect the motion of the particle. All other (distant)
ions’ coordinates are included in $\bm{\psi}$. The interaction between the
distant ions and the mobile particle is sufficiently weak to be neglected when
computing the particle’s trajectory.

We assert that the escape success is determined by a subset of $\bm{\xi}$
coordinates, corresponding to the solid’s ions sufficiently close to the
mobile ion’s position. This assertion allows us to redefine
$\bm{\xi}\rightarrow(\bm{\zeta},\bm{\psi})$, where $\bm{\zeta}$ are the
coordinates that determine the escape and $\bm{\psi}$ are the remaining ones,
as shown in Fig. 1. Moreover, we assume that the newly-defined $\bm{\zeta}$
couples weakly to $\bm{\psi}$, allowing us to neglect the influence of
$\bm{\psi}$ on the time evolution of $\bm{\zeta}$ for sufficiently short
times. Consequently, $\bm{\psi}$ only plays a role in the probability of a
configuration $(\mathbf{V},\mathbf{R},\bm{\zeta})$, but not determining which
configurations result in the ion’s escape. Hence, redefining $\mathcal{A}$, we
have

$\displaystyle\mathcal{P}$
$\displaystyle=\left[\sum_{\mathbf{V},\mathbf{R},\bm{\zeta},\bm{\psi}}e^{-\beta
E(\mathbf{V},\mathbf{R},\bm{\zeta},\bm{\psi})}\right]^{-1}\sum_{(\mathbf{V},\mathbf{R},\bm{\zeta})\in\mathcal{A},\bm{\psi}}e^{-\beta
E(\mathbf{V},\mathbf{R},\bm{\zeta},\bm{\psi})}$
$\displaystyle=\sum_{(\mathbf{V},\mathbf{R},\bm{\zeta})\in\mathcal{A}}P(\mathbf{V},\mathbf{R},\bm{\zeta})\,,$
(2)

where we formally traced out $\bm{\psi}$ and defined the probability of the
$(\mathbf{V},\mathbf{R},\bm{\zeta})$ configuration
$P(\mathbf{V},\mathbf{R},\bm{\zeta})$.

Importantly, $P(\mathbf{V},\mathbf{R},\bm{\zeta})$ cannot necessarily be
assumed to have a Boltzmann dependence on the energy of the reduced system
containing the coordinates $(\mathbf{V},\mathbf{R},\bm{\zeta})$. Moreover,
even if it does, $\mathcal{A}$ can take a complicated shape in the phase space
so that the sum of $P(\mathbf{V},\mathbf{R},\bm{\zeta})$ over $\mathcal{A}$
should not be expected to give $e^{-\beta E_{a}}$, in general. As an
illustrative example, the following Section explores the simplest mobile ion-
lattice configuration given in Fig. 2.

Figure 2: _Single-mass lattice model._ A particle of mass $M$ is constrained
to motion along $x$ and interacts via a coupling $U$ with a gate harmonic
oscillator of mass $m$, which is constrained to motion along $y$. The axis of
the mobile particle is separated from the equilibrium position of the
oscillator by an impact parameter $s$. The gate oscillator has spring constant
$k$, with mode occupation consistent with a temperature $T$.

##  III Numerical simulations

Figure 3: _Phase space._ In a system with $M=1,\lambda=2,\Phi_{0}=10$,
initialized with various $\mathcal{E}_{\mathrm{particle}}$,
$\mathcal{E}_{\mathrm{oscillator}}$, and $\phi$, the mobile particle passes
the oscillator or is reflected. The plots show the domains where the particle
is reflected (orange) and transmitted (green; this is the region
$\mathcal{A}$). Each panel shows a portion of the
$\mathcal{E}_{\mathrm{particle}}$-$\mathcal{E}_{\mathrm{oscillator}}$ space
for a fixed $\phi$. The simulations are performed for $s=0$ (head-on; a–h) and
$s=2$ (offset; i–p).

###  III.1 Equations of motion

Figure 2 shows a system where a particle of mass $M$ moves in one dimension
perpendicular to a harmonic oscillator with mass $m$ and force constant $k$,
interacting with it via a repulsive coupling $U$. If the mobile particle is
launched towards the oscillator from some initial position $x$ with some
initial speed $\dot{x}$, it will either pass it or be reflected. Using the
notation from Sec. II, $\bm{\zeta}$ is a two-component vector containing the
displacement and speed of the oscillator, $y$ and $\dot{y}$ respectively.
Whether or not the particle passes the oscillator depends on the initial
conditions of $x$, $\dot{x}$, $y$, and $\dot{y}$, partitioning the four-
dimensional space $\left\\{x,\dot{x},y,\dot{y}\right\\}$ into two regions, as
discussed in Sec. II. To explore the geometry of this phase space and study
the role of temperature in the transmission, we start by writing down the
Lagrangian for the two-mass system:

$L=\frac{1}{2}M\dot{x}^{2}+\frac{1}{2}m\dot{y}^{2}-\frac{1}{2}ky^{2}-U(x,y)\,,$
(3)

where the oscillator and the mobile particle are constrained to move in
single, perpendicular directions.

From the oscillator frequency $\Omega=\sqrt{k/m}$, we define a time scale
$t_{\mathrm{osc}}=2\pi/\Omega$, a length scale
$l_{\mathrm{osc}}=\sqrt{\hbar/m\Omega}$ (corresponding to the quantum
oscillator length), and an energy scale
$E_{\mathrm{osc}}=\hbar\Omega=kl_{\mathrm{osc}}^{2}$. Redefining $x\rightarrow
xl_{\mathrm{osc}}$, $y\rightarrow yl_{\mathrm{osc}}$, $M\rightarrow mM$, and
$U\rightarrow UE_{\mathrm{osc}}$ then dividing Eq. (3) by $E_{\mathrm{osc}}$,
we get a dimensionless Lagrangian

$L=\frac{1}{8\pi^{2}}M\dot{x}^{2}+\frac{1}{8\pi^{2}}\dot{y}^{2}-\frac{1}{2}y^{2}-U(x,y)\,,$
(4)

where the time derivatives are taken with respect to
$\tau=t/t_{\mathrm{osc}}$. The corresponding equations of motion are

$\displaystyle\frac{M\ddot{x}}{4\pi^{2}}$
$\displaystyle=-\frac{d}{dx}U(x,y)\,,$
$\displaystyle\frac{\ddot{y}}{4\pi^{2}}$
$\displaystyle=-y-\frac{d}{dy}U(x,y)\,.$ (5)

###  III.2 Transmission phase space

For convenience, we express $y$ and $\dot{y}$ using the amplitude and the
phase of the oscillator $y_{0}$ and $\phi$, respectively. Although the phase
space for the two-mass system is four-dimensional, if the mobile particle
starts far enough from the oscillator so that the the interaction vanishes,
modifying the initial separation is equivalent to changing $\phi$, reducing
the effective dimensionality of the phase space.

To illustrate the shape of $\mathcal{A}$ containing all the transmitting
configurations, we choose $U(z)=U_{0}e^{-z^{2}/2\lambda^{2}}$ with
$z=\sqrt{x^{2}+(y-s)^{2}}$, where $s$ is the impact parameter, $\lambda=2$,
and $U_{0}=10$. The motivation behind this choice of $U(z)$ is its simplicity:
one can easily change the strength and extent of the interaction by adjusting
$U_{0}$ and $\lambda$, as was done in other works demonstrating dissipation in
similar systems [21, 22, 23]. We set $M=1$ and use two different impact
parameters ($s=0$ and $s=2$). Next, we take a range of particle energies
$1/10\leq\mathcal{E}_{\mathrm{particle}}\leq 12$, oscillator energies
$0\leq\mathcal{E}_{\mathrm{oscillator}}\leq 20$, and oscillator phases
$\phi\in[0,\pi/4,\dots 7\pi/4]$. Here
$\mathcal{E}_{\mathrm{particle}}=E_{\mathrm{particle}}/E_{\mathrm{osc}}$ and
$\mathcal{E}_{\mathrm{oscillator}}=E_{\mathrm{oscillator}}/E_{\mathrm{osc}}$.
For each triplet
$(\mathcal{E}_{\mathrm{particle}},\mathcal{E}_{\mathrm{oscillator}},\phi)$, we
get the particle’s initial speed, as well as the oscillator’s initial speed
and position. Then, we launch the particle towards the oscillator starting at
$x=-20$, evolving the system in time using the fifth order Runge-Kutta
procedure with the time step $\delta\tau=5\times 10^{-4}$. At each time step,
the velocity and the position of the particle are checked and the time
evolution stops when either the velocity of the particle becomes negative
(reflection) or the particle passes the oscillator (transmission).

We plot the results in Fig. 3 as a collection of domain maps in the
$\mathcal{E}_{\mathrm{particle}}$–$\mathcal{E}_{\mathrm{oscillator}}$ space
for each of the $\phi$’s. We see that the shape of $\mathcal{A}$ is highly
nontrivial, even in this simple model. To calculate the transmission
probability, one needs to integrate $P(\dot{x},y_{0},\phi)$ over
$\mathcal{A}$:

$\mathcal{P}=\sum_{(\dot{x},x,\dot{y},y)\in\mathcal{A}}P(\dot{x},x,\dot{y},y)\equiv\sum_{(\dot{x},y_{0},\phi)\in\mathcal{A}}P(\dot{x},y_{0},\phi)\,.$
(6)

Given the shape of the phase space, even if we assume that the total system
energy is Boltzmann-distributed with $P(\dot{x},y_{0},\phi)\propto
e^{-\mathcal{E}/\omega_{T}}$ (with $\omega_{T}=k_{B}T/E_{\mathrm{osc}}$) so
that constant-$\mathcal{E}$ surfaces form quarter-circles on the
$\mathcal{E}_{\mathrm{particle}}$–$\mathcal{E}_{\mathrm{oscillator}}$ plane,
it is unlikely that the integral will exhibit an Arrhenius behavior. To be
certain, we “integrate” over the phase space by conducting a series of four
numerical experiments, where we send the mobile particle towards the
oscillator and calculate its transmission probability, making different
assumptions about the system.

###  III.3 Stationary oscillator

For the first two experiments, the oscillator starts from rest with zero
displacement so that $y_{0}=0$, making $\phi$ irrelevant. By eliminating two
out of three variables in Eq. (6), we have
$\mathcal{P}=\sum_{\dot{x}\in\mathcal{A}}P(\dot{x})$. In this simple case, for
a repulsive interaction, the particle must have a certain minimum speed so
that $\mathcal{P}=\sum_{\dot{x}\geq\dot{x}_{\mathrm{min}}}P(\dot{x})$,
corresponding to some minimum kinetic energy. If the particle energy follows
the Boltzmann distribution $\propto
e^{-\mathcal{E}_{\mathrm{particle}}/\omega_{T}}$, the probability of
transmission will take the Arrhenius form, given by
$e^{-\mathcal{E}_{\mathrm{min}}/\omega_{T}}$. As we show below, the value of
$\mathcal{E}_{\mathrm{min}}$ is not always easily apparent.

If the oscillator is frozen in place and not allowed to move in response to
its interaction with the particle, the minimum kinetic energy that the
particle must start with to overcome the potential barrier is
$\mathcal{E}_{\mathrm{min}}=U(s)$. To establish the baseline, our first
experiment aims to confirm this known result by considering a range of
$\omega_{T}$’s and sampling 10,000 energies from the Boltzmann distribution at
each $\omega_{T}$. For each of the energies, the particle is launched towards
the fixed oscillator from $x=-20$ and the system is evolved following the same
Runge-Kutta procedure as in Sec. III.2. By counting how many of the 10,000
particles make it past the fixed oscillator at each $\omega_{T}$, we obtain
the transmission probability as a function of temperature, which is plotted in
Fig. 4(a) for the same two impact parameters as in Sec. III.2. The experiment
was repeated three times for different particle masses (1/5, 1, and 5) to
confirm that the result is mass-independent, as expected. By using the same
random seed to generate the initial conditions for the three masses, we
obtained identical results, as can be seen from the perfect overlapping of the
scatter plot symbols. The numerical results show an excellent agreement with
the analytic Arrhenius form, validating our computation approach.

Figure 4: _Stationary oscillator._ A particle starting from $x=-20$, is
launched towards the oscillator, which is initially at rest, either head-on
($s=0$) or offset to the side ($s=2$). The speed of the particle is obtained
by generating an energy $\mathcal{E}$ from $e^{-\mathcal{E}/\omega_{T}}$ and
equating it to the kinetic energy $M\dot{x}^{2}/8\pi^{2}$ to calculate
$\dot{x}$. Three masses are used ($1/5$, $1$, and $5$, referred to as “light”,
“medium”, and “heavy”, respectively). The system is evolved in time until the
particle either passes the oscillator or is reflected. For each mass at each
temperature, 10,000 realizations are performed and the fraction of particles
passing the oscillator is recorded. In panel (a), the oscillator is held
fixed, resulting in a conservative potential for the particle. In panel (b),
the system evolves in accordance with Eq. (5). Solid lines correspond to
$e^{-\Phi(s)/\omega_{T}}$; the dashed line is given by
$e^{-\mathcal{E}_{*}/\omega_{T}}$, where $\mathcal{E}_{*}$ is the minimized
potential energy when the particle is at $x=0$.

For the second experiment, we allow the oscillator to be deflected in response
to the particle. By performing the same steps as for Fig. 4(a), we obtain the
transmission probability and plot the results in Fig. 4(b). Because the
oscillator can only move in the direction perpendicular to the particle’s
trajectory, for $s=0$, there is no force component applied on the oscillator
due to the vanishing derivative, so it remains stationary. Comparing the
results for $s=0$ in panels (a) and (b) of Fig. 4, we confirm that they are
identical due to the same random seed. For a finite impact parameter, on the
other hand, the particle can push the mass laterally, reducing the maximum
value of the interaction between the two objects and lowering the potential
barrier.

In the spirit of NEB, one way to estimate the reduced barrier is to relax the
system by bringing the particle to the point of the closest approach to the
oscillator, minimizing the potential energy by varying $y$, and using this
value for the barrier. Explicitly, we minimize the potential energy by solving
$U^{\prime}(y_{*}-s)+y_{*}=0$. In case there are multiple stable solutions,
corresponding to the oscillator mass being on either side of the particle
trajectory, we take the one that gives the smaller energy. Using this
minimized value as the activation energy
$\mathcal{E}_{\mathrm{min}}=\mathcal{E}_{*}=\Phi(y_{*}-s)+y_{*}^{2}/2$ in the
Arrhenius relation, we plot the “relaxed” theoretical curve in Fig. 4(b) as a
dashed line, along with the “frozen” result for an unrelaxed barrier
$e^{-U(2)/\omega_{T}}$.

Comparing our numerical results for $s=2$ to the Arrhenius plots with the two
activation energies, we observe that the theoretical curves bound the
numerical results. Additionally, performing a power-law fit of
$-\ln\mathcal{P}$ vs. $1/\omega_{T}$ for each mass’s numerical data indicates
that the dependence has an Arrhenius form. This form indicates that there
indeed exists an $\mathcal{E}_{\mathrm{min}}$, but it is mass-dependent and
lies between the relaxed and unrelaxed values. The mass dependence can be
easily understood intuitively: for a given kinetic energy, heavier particles
move slower compared to the light ones. The slow speed allows the oscillator
to move out of the particle’s way and assume a configuration that is similar
to the relaxed one. Lighter particles, on the other hand, move faster and the
oscillator does not have time to respond so that the particles need to
overcome the unrelaxed (frozen oscillator) potential barrier.

There are two key points to take away from these experiments. First, if only
one of the system’s degrees of freedom is probabilistic and follows the
Boltzmann distribution, the transmission probability will take the Arrhenius
form. Second, one cannot generally determine the activation energy accurately
while ignoring kinetic effects, as can be seen from Fig. 4(b). Consequently,
the activation energy obtained from quasi-static methods, such as NEB, is
likely to disagree with experiment-derived values. Next, we show that the
existence of multiple probabilistic degrees of freedom, as would be the case
in real systems, causes the transmission probability to deviate from the
Arrhenius form.

###  III.4 Moving oscillator

Figure 5: _Moving oscillator._ Simulations following the same protocol as
outlined in the caption of Fig. 4, but with the oscillator not starting from
rest. Instead, it is initialized with the initial position
$\sqrt{2n+1}\sin\phi$ and velocity $2\pi\sqrt{2n+1}\cos\phi$. In panel (a),
the probability distribution of the initial configurations is
$P(\dot{x},y_{0},\phi)=P_{\mathrm{particle}}(\dot{x})P_{\mathrm{oscillator}}(y_{0})/2\pi$.
In panel (b), the total energy of the system is distributed as
$e^{-\mathcal{E}/\omega_{T}}$ with the fraction of energy distributed between
the particle and the oscillator being random and uniform. In both cases,
$\phi$ is uniformly distributed between $0$ and $2\pi$.

For the next two experiments, the oscillator is initialized with a finite
speed and displacement for probabilistic $y_{0}$ and $\phi$. In the first
experiment, the particle and the oscillator are treated as two independent
objects following independent Boltzmann distributions. This scenario would
arise if, prior to the particle being launched towards the oscillator, both
objects were kept in contact with separate heat baths. Following the same
procedure as in Sec. III.3, 10,000 realizations are performed at each
$\omega_{T}$ with the particle energy sampled from the Boltzmann distribution
$e^{-\mathcal{E}_{\mathrm{particle}}/\omega_{T}}$, while the oscillator energy
is generated from $e^{-n/\omega_{T}}$ for integer $n$’s. The phase $\phi$ is
picked from a uniform distribution $[0,2\pi]$. The amplitude of the oscillator
then becomes $\zeta=\sqrt{2n+1}$ so that the initial position is
$\zeta\sin\phi$ and the initial speed is $2\pi\zeta\cos\phi$.

The simulations are performed for the same masses and impact parameters used
in Fig. 4. The results for the transmission probability are given in Fig.
5(a), along with the theoretical Arrhenius lines using relaxed and frozen
activation energies. Since the oscillator does not start from rest with
$y_{0}=0$ and can be deflected by the particle even at $s=0$, the Arrhenius
line for the relaxed configuration is also included.

The second experiment treats the particle-oscillator system as a single whole
and assumes that the composite system’s total energy is Boltzmann distributed.
To create such a configuration, one could imagine that the particle’s motion
is bounded by hard walls positioned far from the oscillator. When the particle
bounces off these walls, it exchanges energy with an external bath and later
interacts with the oscillator so, after a long time, the two-body system
acquires the temperature of the bath. Here, the initial conditions for the
10,000 runs at each $\omega_{T}$ are generated by first picking the total
energy $\mathcal{E}$ from the exponential distribution
$e^{-\mathcal{E}/\omega_{T}}$. After generating $\mathcal{E}$, we randomly
pick $n$ between $0$ and $\lfloor\mathcal{E}\rfloor$ and set the energy of the
particle to $\mathcal{E}-n$. Finally, $\phi$ is chosen from a uniform
distribution. The results of this experiment are given in Fig. 5(b).

The two experiments represent maximally different scenarios: either the two
degrees of freedom are completely independent or their composite system
follows a thermal distribution. In a real system, the situation is between the
two extremes. The statistics of the $\bm{\zeta}$ coordinates are not
completely independent from the mobile particle. At the same time, the mobile
ion and its neighborhood do not comprise an isolated system with a Boltzmann
energy distribution. Figure 5 shows that the Arrhenius formula does not
capture the transmission probability for either case for both impact
parameters, regardless of which activation energy is used. We can see that,
for either distribution, while the results for one impact parameter match a
theory prediction reasonably well, the other impact parameter results clearly
deviate from the Arrhenius form. Notably, the impact parameter that is matched
better by a theoretical result is different between the two panels, as is the
method of obtaining the theoretical prediction (relaxed vs. frozen). A closer
look reveals that even the better-matched data exhibit a slight curvature,
which would lead to a temperature-dependent activation energy.

The outcome of the experiments supports the discussion in Sec. II that the
transmission probability will, in general, not follow the Arrhenius form.
There are several implications stemming from this result. First, extrapolating
the results obtained using AIMD to different temperatures is precarious since
the dependence of the activation energy on temperature is generally not known.
Moreover, even if one assumes that the activation energy is constant over a
limited temperature range, it is not evident which temperature the NEB results
should correspond to.

##  IV Conclusions

In this work, we have explored the idea of activation in solid state ionic
electrolytes. By constructing a minimal model and visualizing the
corresponding three-dimensional phase space, we have shown that the portion of
this space that results in activation has a non-trivial shape. Consequently,
even if the probability of each point in the phase space follows the Boltzmann
distribution, it is unlikely that summing over all the points in the
activation region after weighting them with the appropriate Boltzmann factor
will give rise to an Arrhenius dependence associated with activation. To
support this argument, we conducted a series of numerical experiments using
our minimal model to demonstrate that one gets the Arrhenius curve only if
there is a single Boltzmann-distributed degree of freedom and, at best, one
can expect a slowly-varying temperature-dependent “activation energy.” The
form of this variable activation energy is non-trivial even for our simple
configuration, making it difficult to reliably extrapolate the numerical
results obtained at one temperature to different temperature ranges. Finally,
we have shown that one cannot accurately obtain the activation energy by
neglecting the kinetic component of the ionic motion, as is commonly done in
the quasi-static NEB calculations.

The strengths of existing numerical methods make them extremely useful for
rapidly screening large numbers of viable candidates for further experimental
scrutiny. While this work does not undermine the general usefulness of these
methods, based on our findings, we suggest that caution must be taken when
interpreting data from numerical simulations using the Arrhenius form of the
$\sigma$–$T$ relation. Specifically, it is not guaranteed a priori that such a
relation exists for any given system, especially over a wide temperature
range. Moreover, a single, temperature-independent activation energy for ionic
migration is not likely to exist in general. This work points the way to areas
of possible refinement for the theoretical simulation and study of ionic
migration in solid electrolytes.

###### Acknowledgements.

A. R. and K. N. acknowledge the National Research Foundation, Prime Minister
Office, Singapore, under its Medium Sized Centre Programme. A. R. acknowledges
support by Yale-NUS College (through Start-up Grant). B. A. O. acknowledges
support from the M. J. Murdock Charitable Trust.

## References

  * Mahan and Roth [1976] G. D. Mahan and W. L. Roth,  _Superionic Conductors_ , 1st ed. (Plenum Press, 1976).
  * Mehrer [2007] H. Mehrer,  _Diffusion in Solids_ , 1st ed. (Springer, 2007).
  * Bachman  _et al._ [2016] J. C. Bachman, S. Muy, A. Grimaud, H. H. Chang, N. Pour, S. F. Lux, O. Paschos, F. Maglia, S. Lupart, P. Lamp, L. Giordano, and Y. Shao-Horn, Inorganic solid-state electrolytes for lithium batteries: Mechanisms and properties governing ion conduction, Chemical Reviews 116, 140 (2016).
  * Manthiram  _et al._ [2017] A. Manthiram, X. Yu, and S. Wang, Lithium battery chemistries enabled by solid-state electrolytes, Nat. Rev. Mats. 2, 1 (2017).
  * Famprikis  _et al._ [2019] T. Famprikis, P. Canepa, J. A. Dawson, M. S. Islam, and C. Masquelier, Fundamentals of inorganic solid-state electrolytes for batteries, Nature Materials 18, 1278 (2019).
  * Rogez  _et al._ [2000] J. Rogez, P. Knauth, A. Gamier, H. Ghobarkar, and O. Schaf, Determination of the crystallization enthalpies of lithium ion conducting alumino-silicate glasses, J. Non-Cryst. 262, 177 (2000).
  * Wilkening  _et al._ [2008] M. Wilkening, V. Epp, A. Feldhoff, and P. Heitjans, Tuning the li diffusivity of poor ionic conductors by mechanical treatment: High li conductivity of strongly defective litao3 nanoparticles, J. Phys. Chem. C 112, 9291 (2008).
  * Chu  _et al._ [2016] I. H. Chu, H. Nguyen, S. Hy, Y. C. Lin, Z. Wang, Z. Xu, Z. Deng, Y. S. Meng, and S. P. Ong, Insights into the performance limits of the li7p3s11 superionic conductor: A combined first-principles and experimental study, ACS Appl. Mater. Interfaces 10.1021/acsami.6b00833 (2016).
  * Deng  _et al._ [2017] Y. Deng, C. Eames, B. Fleutot, R. David, J. N. Chotard, E. Suard, C. Masquelier, and M. S. Islam, Enhancing the lithium ion conductivity in lithium superionic conductor (lisicon) solid electrolytes through a mixed polyanion effect, ACS Applied Materials and Interfaces 9, 7050 (2017).
  * Krauskopf  _et al._ [2017] T. Krauskopf, C. Pompe, M. A. Kraft, and W. G. Zeier, Influence of lattice dynamics on na+ transport in the solid electrolyte na3ps4-xsex, Chem. Mater 29, 8859 (2017).
  * Deng  _et al._ [2015] Y. Deng, C. Eames, J. N. Chotard, F. Laleìre, V. Seznec, S. Emge, O. Pecher, C. P. Grey, C. Masquelier, and M. S. Islam, Structural and mechanistic insights into fast lithium-ion conduction in li4sio4-li3po4 solid electrolytes, Journal of the American Chemical Society 137, 9136 (2015).
  * Kato  _et al._ [2016] Y. Kato, S. Hori, T. Saito, K. Suzuki, M. Hirayama, A. Mitsui, M. Yonemura, H. Iba, and R. Kanno, High-power all-solid-state batteries using sulfide superionic conductors, Nat. Energy 1, 1 (2016).
  * Fall  _et al._ [2019] B. Fall, P. Prakash, M. R. Gau, S. L. Wunder, A. Venkatnathan, and M. J. Zdilla, Experimental and theoretical investigation of the ion conduction mechanism of tris(adiponitrile)perchloratosodium, a self-binding, melt-castable crystalline sodium electrolyte, Chem. Mater. 10.1021/acs.chemmater.9b02853 (2019).
  * Yashima  _et al._ [2021] M. Yashima, T. Tsujiguchi, Y. Sakuda, Y. Yasui, Y. Zhou, K. Fujii, S. Torii, T. Kamiyama, and S. J. Skinner, High oxide-ion conductivity through the interstitial oxygen site in ba7nb4moo20-based hexagonal perovskite related oxides, Nat. Commum. 12, 1 (2021).
  * Hu  _et al._ [2023] L. Hu, J. Zhu, C. Duan, J. Zhu, J. Wang, K. Wang, Z. Gu, Z. Xi, J. Hao, Y. Chen, J. Ma, J. X. Liu, and C. Ma, Revealing the pnma crystal structure and ion-transport mechanism of the li3ycl6 solid electrolyte, Cell Reports Physical Science 4, 101428 (2023).
  * Byeon and Kim [2021] Y.-W. Byeon and H. Kim, Review on interface and interphase issues in sulfide solid-state electrolytes for all-solid-state li-metal batteries, Electrochem 2, 452 (2021).
  * Qi  _et al._ [2021] J. Qi, S. Banerjee, Y. Zuo, C. Chen, Z. Zhu, M. Holekevi Chandrappa, X. Li, and S. Ong, Bridging the gap between simulated and experimental ionic conductivities in lithium superionic conductors, Materials Today Physics 21, 100463 (2021).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing *, Society for Industrial and Applied Mathematics 59, 10.1137/141000671 (2017).
  * Danisch and Krumbiegel [2021] S. Danisch and J. Krumbiegel, Makie.jl: Flexible high-performance data visualization for julia, Journal of Open Source Software 6, 3349 (2021).
  * Wong [2011] B. Wong, Color blindness, Nature Methods 8, 441 (2011).
  * Rodin  _et al._ [2022] A. Rodin, B. A. Olsen, M. Choi, and A. Tan, Microscopic theory of thermalization in one dimension with nonlinear bath coupling, Physical Review Research 4, 033057 (2022).
  * Mahalingam  _et al._ [2023a] H. Mahalingam, Z. W. Yap, B. A. Olsen, and A. Rodin, Minimal model of drag in one-dimensional crystals, Phys. Rev. Res. 5, 013053 (2023a).
  * Mahalingam  _et al._ [2023b] H. Mahalingam, B. A. Olsen, and A. Rodin, Dissipation and diffusion in one-dimensional solids, Phys. Rev. Res. 5, 033044 (2023b).