# Local-time formula for dissipation in solid ionic electrolytes

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  B. A. Olsen  Department of Physics, Lewis &
Clark College, Portland, Oregon, 97219, USA  A. Ustyuzhanin  Constructor
University, Bremen, Campus Ring 1, 28759, Germany  Institute for Functional
Intelligent Materials, National University of Singapore, 4 Science Drive 2,
Singapore 117544, Singapore  A. Maevskiy  Institute for Functional Intelligent
Materials, National University of Singapore, 4 Science Drive 2, Singapore
117544, Singapore  K. Noori  Institute for Functional Intelligent Materials,
National University of Singapore, 4 Science Drive 2, Singapore 117544,
Singapore

###### Abstract

When ions move through solids, they interact with the solid’s constituent
atoms and cause them to vibrate around their equilibrium points. This
vibration, in turn, modifies the potential landscape through which the mobile
ions travel. Because the present-time potential depends on past interactions,
the coupling is inherently non-local in time, making its numerical and
analytical treatment challenging. For sufficiently slow-moving ions, we
linearize the phonon spectrum to show that these non-local effects can be
ignored, giving rise to a drag-like force. Unlike the more familiar drag
coefficient in liquids, the drag takes on a matrix form due to the crystalline
structure of the framework. We numerically simulate trajectories and
dissipation rates using both the time-local and non-local formulas to validate
our simplification. The time-local formula dramatically reduces the
computational cost of calculating the motion of mobile particle through a
crystalline framework and clearly connects the properties of the material to
the drag experienced by the particle.

##  I Introduction

Understanding ionic motion through solids [1, 2] is essential to developing
high-performance solid-state electrolytes necessary for the operation of
solid-state batteries [3, 4, 5]. A substantial fraction of theoretical work
focuses on elucidating the mechanisms behind ionic current flow [6, 7, 8, 9,
10, 11, 12, 13] and how the properties of the solid framework impact the ionic
transport, in particular the ability of the current-carrying ions to overcome
energy barriers within the host lattice [14, 15, 16, 17]. A complementary and
equally important question is how these ions lose energy to the solid
framework during their motion. Relating this energy loss and, by extension,
the resistance to the material properties is crucial for designing
application-grade electrolytes. The fundamental difficulty of predicting
material resistivity is that the dynamics of the mobile ion and the framework
occur on similar time scales. As a result, the mobile particle both displaces
and experiences recoil from the framework, which is continuously responding to
the mobile particle’s earlier motion, leading to equations of motion that are
not local in time. In this work, we develop an approximation scheme to
circumvent this difficulty, so the dissipation is determined from the
framework’s unperturbed potential landscape and its vibrational modes.

Earlier work [18] introduced an analytic formula for dissipation in three-
dimensional (3D) crystals using a local-time approximation. The results showed
that the drag experienced by the mobile particles depends on the curvature of
the potential generated by the lattice. Additionally, the dissipation rate is
inversely proportional to the cube of speed of sound in the material and also
to the density of the crystal. In other words, the dissipation is lowest in
stiff, dense materials. Due to the complexity of the problem, these results
were not validated using numerical demonstrations.

Following the analytical work of Ref. [18], dissipation has been explored in
the simplest crystal functioning as a solid electrolyte: an infinite one-
dimensional (1D) chain of identical masses interacting with a single mobile
particle moving along the chain [19, 20]. These studies showed that
dissipation is nonlinear, decreases at high particle speeds, and can also be
nonmonotonic in particle speed depending on the system parameters. In
particular, if the phonon spectrum is gapped due to the harmonic confinement
of the chain masses, the dissipation is exponentially suppressed at low speeds
in the absence of thermal vibration.

There are three key takeaways from this work we would like to highlight.
First, by validating the results of Ref. [18] numerically, we confirm the
dependence of dissipation on material properties. Second, by obtaining a time-
local equation of motion, we are able to dramatically accelerate simulations
of the ionic motion. This acceleration is useful because ionic motion is a
probabilistic process and the quality of numerical experiments improves with
simulation time. Finally, by starting with a simpler 1D system and gradually
building up to a more complex 3D configuration, we elucidate the features of
the problem that make the time-local treatment possible, paving the way for
our future work, where we will include thermal motion of the framework.

In this work, we make use of the techniques developed for 1D systems as the
dimensionality constraints are gradually relaxed in order to validate the
local-time approximation formalism. As the first step, we discuss a general
model of single-particle motion through a crystal in Sec. II. With the
formalism established, readers interested primarily in the final result can
skip directly to Eq. (26) and consult earlier portions of Sec. IV for
clarifications as needed. For readers interested in the technical aspects of
the formalism, we suggest proceeding in order. Starting with a 1D chain whose
masses can move in three dimensions and can support longitudinal and
transverse phonon modes, we use numerical and analytical tools to demonstrate
how the system properties determine the energy dissipation of the mobile ion,
with the results given in Sec. III. By exploring the high- and low-speed
limits, we establish the range of validity for the local-time approximation.
Next, by focusing on a prototypical 3D system in Sec. IV, we show that, unlike
the 1D case, where the dissipation diverges at low speeds without harmonic
confinement, the drag experienced by slow-moving ions is proportional to their
speed. Finally, we present and validate a time-local formula that connects the
framework potential landscape to the dissipation. We discuss some of the key
choices and assumptions used in this paper in Sec. V. We propose some avenues
for extension to thermal systems and conclude in Sec. VI.

All of our computations were performed using julia [21]. The plots were made
with the Makie.jl package [22] using a color scheme designed for colorblind
readers [23]. The scripts used for computing and plotting can be found at
https://github.com/rodin-physics/cubic-lattice-ion-motion.

##  II Model

Similar to earlier work, [18, 24, 19, 20] we start with a Lagrangian
describing the motion of a mobile particle of mass $M$ through a framework of
identical masses $m$:

$L=\frac{M}{2}\dot{\mathbf{R}}^{T}\dot{\mathbf{R}}+\frac{m}{2}\dot{\mathbf{r}}^{T}\dot{\mathbf{r}}-\mathbf{r}^{T}\frac{\tensor{V}}{2}\mathbf{r}-U(\mathbf{r},\mathbf{R})\,.$
(1)

Here, $\mathbf{r}=\bigoplus_{j=1}\mathbf{r}_{j}$ and $\mathbf{R}$ are vectors
whose elements are the displacements of all the framework masses from their
equilibrium locations and the position of the mobile particle, respectively.
The first two terms in Eq. (1) give the kinetic energy of the system, and
$\mathbf{r}^{T}\tensor{V}\mathbf{r}/2$ is the internal potential energy for
the framework using the harmonic approximation. Finally,
$U(\mathbf{r},\mathbf{R})$ describes the interaction between the framework and
mobile particle.

Equation (1) makes two important simplifications. First, it focuses on a
single mobile particle. In real ionic conductors, there are multiple current-
carrying ions interacting with each other directly (at small separations) and
indirectly (by deforming the framework). At the same time, the number of
framework ions exceeds the number of mobile ions, resulting in a significant
separation between the current-carrying ions. Hence, Eq. (1) is applicable
when the interaction between mobile ions can be neglected. The second
simplification is the use of the same value $m$ for all framework masses.
Naturally, most materials are composed of multiple atomic types with different
masses. We will demonstrate in Sec. IV, however, that this simplification,
while making the study much more transparent, does not pose a limitation to
our final result, which depends on the density of the system and not the
individual masses.

For convenience, we express the coupling matrix in terms of a characteristic
spring constant $K$: $\tensor{V}=K\tensor{\Xi}$, to be determined at a later
stage. We can then introduce a characteristic frequency $\Omega=\sqrt{K/m}$,
express masses in terms of $m$: $M=\mu\times m$, and express the evolution
time $t$ of the system in terms of the number of periods: $t=\tau\times
2\pi/\Omega$. Finally, we can express all lengths in terms of the quantum
oscillator length associated with $\Omega$, $l=\sqrt{\hbar/m\Omega}$ so that
$\mathbf{R}=\bm{\sigma}l$ and $\mathbf{r}=\bm{\rho}l$. This rescaling turns
the Lagrangian into

$\mathcal{L}=\frac{\mu}{8\pi^{2}}\dot{\sigma}^{2}+\frac{1}{8\pi^{2}}\dot{\rho}^{2}-\bm{\rho}^{T}\frac{\tensor{\Xi}}{2}\bm{\rho}-\Phi(\bm{\rho},\bm{\sigma})\,,$
(2)

where we divided both sides of Eq. (1) by the characteristic energy scale
$\hbar\Omega$, defining $\Phi=U/\hbar\Omega$. The corresponding equations of
motion for the mobile particle and the framework are

$\displaystyle\frac{\mu}{4\pi^{2}}\ddot{\bm{\sigma}}$
$\displaystyle=-\nabla_{\bm{\sigma}}\Phi\left(\bm{\rho},\bm{\sigma}\right)\,,$
(3) $\displaystyle\frac{1}{4\pi^{2}}\ddot{\bm{\rho}}$
$\displaystyle=-\tensor{\Xi}\bm{\rho}-\nabla_{\bm{\rho}}\Phi\left(\bm{\rho},\bm{\sigma}\right)\,,$
(4)

which can be solved in a variety of ways. In this work, we numerically solve
the equations of motion using a Runge-Kutta method. To explore the dynamics of
the framework analytically, however, we also find a formal solution to Eq.
(4). Dropping the last term yields an eigenvalue equation

$\frac{1}{4\pi^{2}}\ddot{\bm{\rho}}=-\omega_{j}^{2}\bm{\rho}=-\tensor{\Xi}\bm{\rho}\,,$
(5)

with normalized eigenvectors $\bm{\varepsilon}_{j}$ and corresponding
eigenvalues $\omega_{j}^{2}$. The general result for $\bm{\rho}$ becomes
$\tensor{\varepsilon}\bm{\zeta}(\tau)$, where $\bm{\zeta}(\tau)$ is a column
vector of normal coordinates giving the amplitude of each mode and
$\tensor{\varepsilon}=[\bm{\varepsilon}_{1},\bm{\varepsilon}_{2},\dots]$ is a
row of column vectors $\bm{\varepsilon}_{j}$. Due to the lattice periodicity,
the eigenvector $\varepsilon_{j}$ is given by
$[e^{i\mathbf{L}_{1}\cdot\mathbf{q}_{j}},e^{i\mathbf{L}_{2}\cdot\mathbf{q}_{j}},\dots]\otimes\bm{\eta}_{j}/\sqrt{N}$.
Here, $\mathbf{L}_{k}$ are the coordinates of the $N\rightarrow\infty$ unit
cells, $\mathbf{q}_{j}$ is the momentum associated with the $j$th mode, and
$\bm{\eta}_{j}$ is a $3n$-dimensional vector with $n$ equal to the number of
masses per unit cell.

Using $\bm{\rho}(\tau)=\tensor{\varepsilon}\bm{\zeta}(\tau)$, Eq. (4) becomes

$\frac{1}{4\pi^{2}}\ddot{\bm{\zeta}}=-\tensor{\omega}^{2}\bm{\zeta}-\tensor{\varepsilon}^{\dagger}\nabla_{\bm{\rho}}\Phi\left(\bm{\rho},\bm{\sigma}\right)\,,$
(6)

where
$\tensor{\omega}^{2}=\tensor{\varepsilon}^{\dagger}\tensor{\Xi}\tensor{\varepsilon}$
(using $\tensor{\varepsilon}^{-1}=\tensor{\varepsilon}^{\dagger}$) is a
diagonal matrix of the squared eigenfrequencies. For each $\zeta_{j}$, the
equation of motion takes the form of a forced harmonic oscillator
$\ddot{\zeta}_{j}=-4\pi^{2}\omega_{j}^{2}\zeta_{j}-4\pi^{2}f_{j}$, which can
be solved using the Green’s function

$G_{j}(\tau,\tau^{\prime})=\frac{\sin\left[2\pi\omega_{j}(\tau-\tau^{\prime})\right]}{2\pi\omega_{j}}\Theta(\tau-\tau^{\prime})\,,$
(7)

to obtain

$\displaystyle\zeta_{j}(\tau)$ $\displaystyle=\zeta_{j}^{H}(\tau)$
$\displaystyle-2\pi\int^{\tau}d\tau^{\prime}\frac{\sin\left[2\pi\omega_{j}(\tau-\tau^{\prime})\right]}{\omega_{j}}\bm{\varepsilon}^{\dagger}_{j}\nabla_{\bm{\rho}}\Phi[\bm{\rho}(\tau^{\prime}),\bm{\sigma}(\tau^{\prime})]\,,$
(8)

where $\zeta_{j}^{H}(\tau)$ is the homogeneous solution. Finally, from Eq.
(8), we get

$\displaystyle\bm{\rho}(\tau)$
$\displaystyle=\tensor{\varepsilon}\bm{\zeta}^{H}(\tau)$
$\displaystyle-2\pi\int^{\tau}d\tau^{\prime}\tensor{\Gamma}(\tau-\tau^{\prime})\nabla_{\bm{\rho}}\Phi[\bm{\rho}(\tau^{\prime}),\bm{\sigma}(\tau^{\prime})]\,,$
(9) $\displaystyle\tensor{\Gamma}(\tau)$
$\displaystyle=\sum_{j}\bm{\varepsilon}_{j}\bm{\varepsilon}_{j}^{\dagger}\frac{\sin\left(2\pi\omega_{j}\tau\right)}{\omega_{j}}\,.$
(10)

Equation (9) gives $\bm{\rho}(\tau)$ as a combination of a homogeneous
trajectory, determined by the initial conditions, and a “memory term”
originating from the framework’s interaction with the mobile particle.

As the particle moves through the framework, it loses energy by doing work on
the framework masses. In the absence of thermal motion, the work done up to
time $\tau$ is given by [20]

$\displaystyle\Delta(\tau)$
$\displaystyle=\int^{\tau}d\tau^{\prime}\left\\{-\nabla_{\bm{\rho}}\Phi\left[\bm{\rho}(\tau^{\prime}),\bm{\sigma}(\tau^{\prime})\right]\right\\}\cdot\dot{\bm{\rho}}(\tau^{\prime})$
$\displaystyle=2\pi^{2}\sum_{j}\left|\int^{\tau}d\tau^{\prime}\nabla_{\bm{\rho}}\Phi\left[\bm{\rho}(\tau^{\prime}),\bm{\sigma}(\tau^{\prime})\right]\cdot\bm{\varepsilon}_{j}e^{2\pi
i\omega_{j}\tau^{\prime}}\right|^{2}\,.$ (11)

In the first line, the term inside the curly braces is the force on the
framework masses and $\dot{\bm{\rho}}(\tau^{\prime})$ their velocity so that
the product of these two quantities gives the power exerted by the particle on
the framework. The second line is obtained by substituting the time derivative
of the second term in Eq. (9) for the velocity vector. In the following two
sections, we use Eq. (11) as the starting point to study the analytical form
of the dissipation.

##  III 1D System

Figure 1: _Schematic of 1D model._ A single mobile particle of mass $M$,
confined to motion along the $z$-axis, interacts via a potential $U$ with a
single mass which is part of an infinite periodic 1D chain, each of mass $m$.
Each chain mass experiences harmonic confinement along all three directions,
with spring constants $\kappa_{x},\kappa_{y},\kappa_{z}$. The chain masses are
also coupled to their nearest neighbors with effective spring constants
$k_{x},k_{y},k_{z}$. In our simulations, the modes are degenerate, so
$\kappa_{\parallel}=\kappa_{\perp}$ and $k_{\parallel}=k_{\perp}$.

To develop intuition for this form of dissipation, we start with a relatively
simple system composed of equally-spaced identical masses arranged in a line,
as shown in Fig. 1. At equilibrium, each mass is confined by a harmonic
potential in all three orthogonal directions. Although the chain’s masses can
move in three dimensions, its vibrational modes are characterized by one-
dimensional momentum. This configuration is distinct from the scenario
considered in Refs. [19, 20], where the masses’ motion was restricted to the
direction parallel to the chain. In this section, unless otherwise stated, “1D
chain” refers to a system as shown in Fig. 1.

Assuming that the 1D chain of masses runs along the $z$ axis, the force
constant in the longitudinal direction is set to $\kappa_{\parallel}$ and for
each of the transverse directions to $\kappa_{\perp}$. Additionally, each mass
is connected to its immediate neighbors, setting the force constant for the
displacement along (perpendicular to) the chain to $k_{\parallel}$
($k_{\perp}$), giving rise to a single longitudinal branch and two branches of
transverse modes. The transverse modes of the chain are distinct from the
flexural modes observed in systems like graphene, which have a quadratic
dispersion at low momenta. Here, the dispersion for all three branches,
labeled by the index $u$, is given by

$\omega_{u,j}^{2}=\kappa_{u}+4k_{u}\sin^{2}\left(\frac{\pi j}{N}\right)\,.$
(12)

Following earlier work [24, 19, 20], the interaction between the mobile
particle and the framework is given by the sum of pairwise interaction terms,
each dependent on the separation between the mobile particle and individual
chain masses. The mass-particle coupling is set to
$\Phi(\bm{\sigma}-\bm{\rho}_{j})=\Phi_{0}e^{-|\bm{\sigma}-\bm{\rho}_{j}|^{2}/2\lambda^{2}}$
for the purposes of demonstration as this potential form is easily tractable
analytically and allows for an easy adjustment of the interaction extent and
magnitude.

The pairwise interaction and the Gaussian form are major simplifications of
the much more complicated $\Phi(\bm{\rho},\bm{\sigma})$ in real materials,
where _ab initio_ calculations are typically required to obtain this energy.
We discuss the implications of using the simple Gaussian form for more
realistic systems in Sec. V.

Figure 2: _Dissipation in 1D._ After passing a single framework mass, the
mobile particle, with initial speed $\sigma_{0}$, loses energy due to its
interaction with that mass. The analytic data were obtained from Eq. (14), the
numeric results were computed by integrating Eq. (13) using quadratures, and
the simulation was performed as described in the main text using Eqs. (3)–(4).
For all panels, $k_{\parallel}=k_{\perp}=99/4$, the chain length was set to
100 with periodic boundary conditions, $\delta\tau=10^{-4}$, and
$\Phi_{0}=0.01$. $\bm{\sigma}_{0}=s\hat{y}$ with $s=0$ for the head-on cases
(left column) and $s=1$ for the offset cases (right column). In the top row,
$\kappa_{\parallel}=\kappa_{\perp}=1$; in the bottom row,
$\kappa_{\parallel}=\kappa_{\perp}=0$. For $\kappa\neq 0$, we see dissipation
at low speeds vanishes. Conversely, for the unconfined chain,
$\Delta\propto\dot{\sigma}^{-1}$ at low speeds, leading to flat
$\dot{\sigma}\Delta$ curves. At high speeds, the dissipation decreases for all
configurations.

To explore the dissipative behavior of this system, we performed several sets
of numerical simulations with different system parameters. For simplicity, we
kept the confining potential and the mass coupling isotropic so that
$\kappa=\kappa_{\perp}=\kappa_{\parallel}$ and $k=k_{\perp}=k_{\parallel}$.
With the framework initially at rest, we introduced a particle with mass
$\mu=1$ and speed $\dot{\sigma}$ moving in the positive $z$ direction a
distance $7\lambda$ to the left of the interacting mass with one of two impact
parameters $\bm{\sigma}_{0}=s\hat{y}$: $s=0$ or $1$. We then evolved the
system following Eqs. (3) and (4) using the fifth order Runge-Kutta method
with time step $\delta\tau=10^{-4}$, terminating when the particle was more
than $7\lambda$ to the right of the interacting mass’s initial position. We
computed the loss $\Delta$ by taking the difference between the initial and
final kinetic energies of the particle, given by
$\mu\dot{\sigma}^{2}/8\pi^{2}$. By repeating the simulation for a range of
initial speeds, we found the speed-dependent loss for a particular set of
framework and interaction parameters, with results shown in Fig. 2.

We can gain a better understanding of the nontrivial loss profiles in Fig. 2
by taking advantage of the Gaussian interaction to analyze Eq. (11).
Unfortunately, in its present form, Eq. (11) is not tractable because we need
$\bm{\sigma}(\tau)$ and $\bm{\rho}(\tau)$ which require solving the equations
of motion in Eqs. (3) and (4). However, if we knew the trajectories, we could
directly compute the dissipation. Making Eq. (11) useful requires several
simplifying assumptions. First, we take $\bm{\rho}$ to be sufficiently small
to replace $\bm{\rho}\rightarrow\bm{0}$ inside the interaction term of Eq.
(11). Second, we assume the velocity of the particle is constant during the
time window of interest, leading to
$\bm{\sigma}(\tau)=\bm{\sigma}_{0}+\tau\dot{\bm{\sigma}}_{0}=s\hat{y}+\dot{\sigma}_{0}\tau\hat{z}$.
This approach was previously used in Refs. [19, 20] in the context of “true”
one-dimensional systems where it was shown that this approximation provides
reasonable results even when the velocity varies substantially. Extending the
time integration limits to $(-\infty,\infty)$ to obtain the energy lost by the
particle during a single “collision” with a chain mass and taking the time
integral yields

$\displaystyle\Delta$
$\displaystyle=\frac{2\pi^{2}}{N}\sum_{j}\left|\frac{1}{\dot{\sigma}}\bm{\eta}_{j}\cdot\left[\nabla_{\perp}+i\hat{z}\frac{2\pi\omega_{j}}{\dot{\sigma}}\right]\tilde{\Phi}_{2\pi\omega_{j}/\dot{\sigma}}\left(\bm{\sigma}_{0}\right)\right|^{2}$
$\displaystyle=\frac{2\pi^{2}}{N}\sum_{j}\left|\nabla_{\perp}\tilde{\Phi}_{2\pi\omega_{j}/\dot{\sigma}}\left(\bm{\sigma}_{0}\right)\cdot\bm{\eta}_{j}\frac{1}{\dot{\sigma}}\right|^{2}$
$\displaystyle\quad+\frac{2\pi^{2}}{N}\sum_{j}\left|\tilde{\Phi}_{2\pi\omega_{j}/\dot{\sigma}}\left(\bm{\sigma}_{0}\right)\frac{2\pi\omega_{j}}{\dot{\sigma}^{2}}\hat{z}\cdot\bm{\eta}_{j}\right|^{2}\,.$
(13)

where $\tilde{\Phi}_{p}(x,y)=\int dz\,e^{ipz}\Phi(x,y,z)$ is the
$z$-coordinate Fourier transform of $\Phi(x,y,z)$. For the Gaussian
interaction considered here,
$\tilde{\Phi}_{2\pi\omega_{j}/\dot{\sigma}}(\bm{\sigma}_{0})=\Phi_{0}\sqrt{2\pi}\lambda\exp\left(-\frac{\sigma_{0}^{2}+4\pi^{2}\omega_{j}^{2}\lambda^{4}/\dot{\sigma}^{2}}{2\lambda^{2}}\right)$,
setting $\nabla_{\perp}\rightarrow-\bm{\sigma}_{0}/\lambda^{2}$. The two terms
of Eq. (13) correspond to two dissipation channels: the first for transverse
modes and the second for longitudinal. The $\bm{\sigma}_{0}\cdot\bm{\eta}_{j}$
term in the transverse channel arising after the replacement of the gradient
operator, combined with the cylindrical symmetry of the system, leads to only
the transverse branch with polarization along $\bm{\sigma}_{0}$ contributing
to dissipation.

Equation (13) can be evaluated numerically by computing $\bm{\eta}_{j}$ and
$\omega_{j}$ from the trivial diagonal dynamical matrix with Eq. (12) for each
direction on the diagonal. We plot the numerically computed values in Fig. 2
and see an excellent agreement with the simulation results. The mild deviation
seen at small speeds is a consequence of the kinetic energy ($\approx 0.05$
for the slowest speed considered, $\dot{\sigma}=2$) becoming comparable to
$\Phi_{0}=0.01$, violating the constant-speed assumption we employed above.

The simplicity of the system allows us to take the integral in Eq. (13)
analytically. Because the system contains a single mass per unit cell, the
polarization vectors $\bm{\eta}_{j}$ contain two 0’s and a single 1
corresponding to the polarization direction, leading to

$\displaystyle\Delta$
$\displaystyle=\Lambda^{2}\Phi_{0}^{2}e^{-\frac{\sigma_{0}^{2}}{\lambda^{2}}}\frac{\sigma_{0}^{2}}{\lambda^{4}}\int_{0}^{\pi}d\theta\,e^{-\Lambda^{2}\omega_{\perp}^{2}(\theta)}$
$\displaystyle+\frac{\Lambda^{4}}{\lambda^{2}}\Phi_{0}^{2}e^{-\frac{\sigma_{0}^{2}}{\lambda^{2}}}\int_{0}^{\pi}d\theta\,\omega^{2}_{\parallel}(\theta)e^{-\Lambda^{2}\omega_{\parallel}^{2}(\theta)}$
$\displaystyle=\pi\Lambda^{2}\Phi_{0}^{2}e^{-\frac{\sigma_{0}^{2}}{\lambda^{2}}}\frac{\sigma_{0}^{2}}{\lambda^{4}}e^{-\Lambda^{2}(\kappa_{\perp}+2k_{\perp})}I_{0}(2\Lambda^{2}k_{\perp})$
$\displaystyle+\pi\frac{\Lambda^{4}}{\lambda^{2}}\Phi_{0}^{2}e^{-\frac{\sigma_{0}^{2}}{\lambda^{2}}}e^{-\Lambda^{2}(\kappa_{\parallel}+2k_{\parallel})}$
$\displaystyle\times\left[(\kappa_{\parallel}+2k_{\parallel})I_{0}\left(2\Lambda^{2}k_{\parallel}\right)-2k_{\parallel}I_{1}\left(2\Lambda^{2}k_{\parallel}\right)\right]\,,$
(14)

where $\Lambda=2\pi\lambda/\dot{\sigma}$ and $I_{n}$ are modified Bessel
functions of the first kind. Plotting Eq. (14) in Fig. 2, we see an
essentially perfect agreement with the numerical integrals, indicating good
convergence.

At high speeds, $\Lambda\rightarrow 0$, resulting in

$\Delta_{\mathrm{fast}}=\pi\frac{\Lambda^{2}}{\lambda^{2}}\Phi_{0}^{2}e^{-\frac{\sigma_{0}^{2}}{\lambda^{2}}}\left[\frac{\sigma_{0}^{2}}{\lambda^{2}}+\Lambda^{2}(\kappa_{\parallel}+2k_{\parallel})\right]\,,$
(15)

showing that the dissipation due to the transverse channel scales as
$1/\dot{\sigma}^{2}$, while that of the longitudinal channel scales as
$1/\dot{\sigma}^{4}$. Naturally, the transverse channel is only possible for
nonzero $\sigma_{0}$ and becomes dominant when
$\frac{\sigma_{0}^{2}}{\lambda^{2}}>\Lambda^{2}(\kappa_{\parallel}+2k_{\parallel})$,
corresponding to
$\dot{\sigma}>2\pi\lambda^{2}\sqrt{\kappa_{\parallel}+2k_{\parallel}}/\sigma_{0}$.
The suppression of $\Delta$ with speed is clearly illustrated by Fig. 2.

Conversely, at low speeds, $\Lambda\rightarrow\infty$, and

$\displaystyle\Delta_{\mathrm{slow}}$ $\displaystyle=\pi\Phi_{0}^{2}\Lambda
e^{-\frac{\sigma_{0}^{2}}{\lambda^{2}}}\frac{1}{2\sqrt{\pi}\lambda^{2}}$
$\displaystyle\times\Bigg{[}\frac{\sigma_{0}^{2}}{\lambda^{2}}\frac{e^{-\Lambda^{2}\kappa_{\perp}}}{\sqrt{k_{\perp}}}+\left(\Lambda^{2}\kappa_{\parallel}+\frac{1}{2}\right)\frac{e^{-\Lambda^{2}\kappa_{\parallel}}}{\sqrt{k_{\parallel}}}\Bigg{]}\,.$
(16)

Equation (16) demonstrates a qualitatively different dissipation behavior
depending on whether the masses are confined ($\kappa>0$) or not. For finite
$\kappa$, the dissipation becomes exponentially suppressed due to the
$\exp\left[-\left(2\pi\lambda\sqrt{\kappa}/\dot{\sigma}\right)^{2}\right]$
term. Physically, the time that it takes the particle to pass the chain mass
$\sim\lambda/\dot{\sigma}$ is much greater than the period of the slowest
vibrational modes $\sim 1/\sqrt{\kappa}$. Therefore, there is no coupling to
the chain modes and, consequently, no dissipation.

Without confinement, $\Delta_{\mathrm{slow}}\propto\dot{\sigma}^{-1}$.
Physically, at $\kappa=0$, there are zero-frequency modes, corresponding to
translational motion, and the slow particle can couple to them. Therefore, the
particle traveling at a constant speed can be thought as “dragging” the entire
infinitely long chain along with itself. This divergence is the consequence of
the framework’s low dimensionality: unlike 3D systems, where the density of
states of the zero-frequency modes vanishes, it remains finite in 1D.

For $\kappa=0$, Eq. (14) indicates that the low-speed limit is applicable when
$2\Lambda^{2}k\gg 1$, corresponding to
$2\pi^{2}\lambda^{2}\omega_{\mathrm{max}}^{2}/\dot{\sigma}^{2}\gg 1$, where
$\omega_{\mathrm{max}}=2\sqrt{k}$ is the maximum frequency of the band. The
$z$-direction Fourier transform of the interaction term also shows that the
particle’s coupling to a mode $j$ is exponentially suppressed by
$\exp\left(-2\pi^{2}\lambda^{2}\omega_{j}^{2}/\dot{\sigma}^{2}\right)$. Thus,
at slow speeds, only low-frequency modes contribute to dissipation so that the
actual dispersion of the high-frequency modes is irrelevant, allowing us to
linearize the phonon dispersion. Returning to Eq. (11) with the approximation
$\bm{\rho}\rightarrow\bm{0}$ and integrating over the linearized spectrum,
then taking the $(-\infty,\infty)$ time integral yields

$\displaystyle\Delta$
$\displaystyle=\frac{2\pi^{2}}{N}\sum_{j}\sum_{u=\perp,\parallel}\int_{-\infty}^{\infty}d\tau^{\prime}d\tau
e^{2\pi i(\tau^{\prime}-\tau)2\sqrt{k_{u}}\frac{\pi j}{N}}$
$\displaystyle\times\partial_{u}\Phi\left[\bm{\sigma}(\tau^{\prime})\right]\partial_{u}\Phi\left[\bm{\sigma}(\tau)\right]$
$\displaystyle=4\pi^{3}\sum_{u=\perp,\parallel}\int_{-\infty}^{\infty}d\tau^{\prime}d\tau\partial_{u}\Phi\left[\bm{\sigma}(\tau^{\prime})\right]\partial_{u}\Phi\left[\bm{\sigma}(\tau)\right]$
$\displaystyle\times\delta\left[4\pi^{2}(\tau^{\prime}-\tau)\sqrt{k_{u}}\right]$
$\displaystyle=\sum_{u=\perp,\parallel}\frac{\pi}{\sqrt{k_{u}}}\int_{-\infty}^{\infty}d\tau\left[\partial_{u}\Phi\left[\bm{\sigma}(\tau)\right]\right]^{2}$
$\displaystyle=\pi^{3/2}\Phi_{0}^{2}\frac{e^{-\frac{\sigma_{0}^{2}}{\lambda^{2}}}}{\lambda\dot{\sigma}}\left(\frac{1}{2\sqrt{k_{\parallel}}}+\frac{\sigma_{0}^{2}}{\lambda^{2}}\frac{1}{\sqrt{k_{\perp}}}\right)\,,$
(17)

which is identical to Eq. (16) with $\kappa=0$. The penultimate line in Eq.
(17) shows that, for slow particles, the loss expression becomes local in
time. Based on this result, we will linearize the phonon dispersion for a more
complicated 3D system in the following section. This approach will allow us to
bypass the more complicated approach in Eq. (13) to directly obtain a local-
time result for the dissipation.

##  IV 3D System

Figure 3: _Schematic of 3D system._ We consider four different scenarios for
interactions between the mobile ion and 3D the framework, which is a simple
cubic lattice in our simulations. In (a), the ion interacts with a single
framework mass, and has impact parameter $s=0$. In (b), the ion interacts with
a single framework mass, but has impact parameter $\mathbf{s}=\hat{y}$. In
(c), the ion interacts with four of the framework masses on a single face of
the unit cell perpendicular to the ion’s velocity, and passes through the
center of that face. In (d), the ion interacts with all eight framework masses
in a single unit cell, and passes through the center of two opposite faces.

Before delving into the details of energy loss by a mobile particle in 3D, we
explore how the system dimensionality impacts dissipation using a set of
simulations similar to the ones performed in Sec. III. To keep our discussion
as transparent as possible, we choose a cubic lattice for the 3D framework. In
the harmonic approximation, the framework masses can be viewed as connected by
springs with natural length equal to the framework separation at equilibrium.
Thus, the pairwise potential energy arising from the displacement of two
masses is

$\displaystyle P_{mn}$
$\displaystyle=\frac{k_{mn}}{2}\left(D_{0}-|\mathbf{D}_{0}-\bm{\rho}_{m}+\bm{\rho}_{n}|\right)^{2}$
$\displaystyle\approx\frac{k_{mn}}{2}\left[\hat{\mathbf{D}}_{0}\cdot\left(\bm{\rho}_{m}-\bm{\rho}_{n}\right)\right]^{2}\,,$
(18)

where $\mathbf{D}_{0}$ is the vector connecting the equilibrium positions of
framework masses $m$ and $n$. The approximation in the second line holds when
$D_{0}\gg|\bm{\rho}_{m}-\bm{\rho}_{n}|$. Assuming that each mass couples only
to its first and second nearest neighbors with force constants $k_{1}$ and
$k_{2}$ respectively, the dynamical matrix becomes [25]

$\displaystyle\mathcal{D}(\mathbf{q})$
$\displaystyle=2k_{1}\begin{pmatrix}1-\cos q_{x}&0&0\\\ 0&1-\cos q_{y}&0\\\
0&0&1-\cos q_{z}\end{pmatrix}$ $\displaystyle+2k_{2}\begin{pmatrix}2-\cos
q_{x}\cos q_{y}-\cos q_{x}\cos q_{z}&\sin q_{x}\sin q_{y}&\sin q_{x}\sin
q_{z}\\\ \sin q_{y}\sin q_{x}&2-\cos q_{y}\cos q_{x}-\cos q_{y}\cos q_{z}&\sin
q_{y}\sin q_{z}\\\ \sin q_{z}\sin q_{x}&\sin q_{z}\sin q_{y}&2-\cos q_{z}\cos
q_{x}-\cos q_{z}\cos q_{y}\end{pmatrix}\,,$ (19)

where $q_{x,y,z}=2\pi\left[1,2,\dots,N_{x,y,z}\right]/N_{x,y,z}$.

Equation (19) demonstrates why the next-nearest neighbor coupling is required:
keeping only $k_{1}$ decouples the Cartesian coordinates and gives rise to
modes with polarization only along the crystal’s axes, producing degenerate
zero-frequency modes along $\Gamma X$, $\Gamma Y$, and $\Gamma Z$ directions.
Introducing $k_{2}$ eliminates these nodal lines. To keep the energy scales
similar to the 1D case, we set $k_{1}=15$ and $k_{2}=5$. The framework masses
are not confined, resulting in a phonon spectrum that is not gapped.

Identically to the 1D case, the mobile particle interacts with a single
framework mass via a Gaussian term and is constrained to move in one
dimension, parallel to one of the framework’s axes, as shown in Fig. 3(a) and
(b). The simulation protocol is the same as in Sec. III, with the particle
starting $7\lambda$ to the left of the interacting mass and the Runge-Kutta
evolution terminating once the particle reaches the position $7\lambda$ to the
right of it. The results of these simulations are shown in Fig. 4.

Figure 4: _Dissipation in 3D with a single interaction._ When a mobile ion
with initial speed $\dot{\sigma}_{0}$ passes through a unit cell of a simple
cubic lattice, as shown in Fig. 3(a) and (b), it loses energy due to
interaction with a single mass. We computed numeric results by integrating Eq.
(13) using cubatures and obtained the the low-speed asymptotic behavior from
Eq. (23), showing that $\Delta\propto\dot{\sigma}$ at low speeds. The
simulation was performed using Eqs. (3)–(4) with the details given in the main
text. For both panels, the nearest (next-nearest) neighbor force constant was
$k_{1}=15$ ($k_{2}=5$), the system contains $50\times 50\times 50$ masses with
periodic boundary conditions, uses time step $\delta\tau=10^{-3}$, and has
Gaussian interactions with amplitude $\Phi_{0}=0.01$.
$\bm{\sigma}_{0}=s\hat{y}$ with $s=0$ for the head-on collisions and $s=1$ for
the offset configurations. The low-speed behavior is captured well by the
analytic formula for a range of simulation parameters.

Just as for the 1D case, the loss decays at high speeds. A crucial difference
occurs at low speeds: $\Delta/\dot{\sigma}_{0}$ approaches a constant,
indicating that the dissipation is proportional to the particle speed, instead
of diverging as its reciprocal (as is the case for 1D). This proportionality
to speed is reminiscent of drag and is a more intuitive outcome than the
divergent behavior for the 1D chain.

We confirm our simulation results with numeric calculations using Eq. (13),
where we computed the eigenstates using the dynamical matrix in Eq. (19) and
integrated over the Brillouin zone using cubatures. The results of these
numeric calculations are shown as open squares in Fig. 4, and are nearly
identical to the simulation results. They differ only at very small speeds,
where the constant-$\dot{\sigma}$ assumption we used to obtain Eq. (13) is
violated.

To explore the low-speed limit, we integrate Eq. (11) by parts over both time
variables assuming that $\bm{\rho}\rightarrow\bm{0}$ and that the mobile
particle interacts with a single framework mass:

$\displaystyle\Delta$
$\displaystyle=\frac{2\pi^{2}}{N}\sum_{j}\int_{-\infty}^{\infty}d\tau
d\tau^{\prime}\nabla\Phi\left[\bm{\sigma}(\tau^{\prime})\right]^{T}\bm{\eta}_{j}$
$\displaystyle\times\bm{\eta}_{j}^{T}\nabla\Phi\left[\bm{\sigma}(\tau)\right]e^{2\pi
i\omega_{j}(\tau^{\prime}-\tau)}$
$\displaystyle=\frac{1}{N}\sum_{j}\frac{1}{2\omega_{j}^{2}}\int_{-\infty}^{\infty}d\tau
d\tau^{\prime}\dot{\bm{\sigma}}(\tau^{\prime})\cdot\mathbf{H}\Phi\left[\bm{\sigma}(\tau^{\prime})\right]\bm{\eta}_{j}$
$\displaystyle\times\bm{\eta}_{j}^{T}\mathbf{H}\Phi\left[\bm{\sigma}(\tau)\right]\cdot\dot{\bm{\sigma}}(\tau)e^{2\pi
i\omega_{j}(\tau^{\prime}-\tau)}\,,$ (20)

where $\mathbf{H}$ denotes the Hessian operator.

To perform the summation over modes, we first write the sum as an integral:

$\frac{1}{2N}\sum_{j}\frac{\bm{\eta}_{j}\bm{\eta}_{j}^{T}}{\omega_{j}^{2}}e^{2\pi
i\omega_{j}\tau}=\frac{1}{16\pi^{3}}\int
d\mathbf{q}\frac{\bm{\eta}_{\mathbf{q}}\bm{\eta}_{\mathbf{q}}^{T}}{\omega_{\mathbf{q}}^{2}}e^{2\pi
i\omega_{\mathbf{q}}\tau}\,.$ (21)

Next, we write $q_{x}\rightarrow q\cos\phi\sin\theta$, $q_{y}\rightarrow
q\sin\phi\sin\theta$, and $q_{z}\rightarrow q\cos\theta$ and linearize the
spectrum by expanding the dynamical matrix $\mathcal{D}(\mathbf{q})$ to
leading order in $q$, which is $q^{2}$, so that
$\mathcal{D}(\mathbf{q})\approx q^{2}\mathcal{S}(\theta,\phi)$. In this limit,
the eigenvectors of $\mathcal{D}(\mathbf{q})$ and $\mathcal{S}(\theta,\phi)$
are identical and can be labeled by their branches and directions:
$\bm{\eta}_{u,\theta,\phi}$. The frequencies $\omega_{\mathbf{q}}$, on the
other hand, are $2\pi q\sqrt{k_{u,\theta,\phi}}$, where $k_{u,\theta,\phi}$ is
the eigenvalue corresponding to $\mathcal{S}(\theta,\phi)$’s eigenvector
$\bm{\eta}_{u,\theta,\phi}$, leading to

$\displaystyle\frac{1}{2N}\sum_{j}\frac{\bm{\eta}_{j}\bm{\eta}_{j}^{T}}{\omega_{j}^{2}}e^{2\pi
i\omega_{j}\tau}$ $\displaystyle=$
$\displaystyle\sum_{u}\int_{0}^{\infty}dq\int_{0}^{\pi}d\theta\sin\theta\oint
d\phi\frac{\bm{\eta}_{u,\theta,\phi}\bm{\eta}_{u,\theta,\phi}^{T}}{16\pi^{3}k_{u,\theta,\phi}}e^{2\pi
i\sqrt{k_{u,\theta,\phi}}q\tau}$ $\displaystyle=$
$\displaystyle\sum_{u}\int_{0}^{\pi}d\theta\sin\theta\oint
d\phi\frac{\bm{\eta}_{u,\theta,\phi}\bm{\eta}_{u,\theta,\phi}^{T}}{16\pi^{2}k_{u,\theta,\phi}}\delta\left(2\pi\sqrt{k_{u,\theta,\phi}}\tau\right)\,.$
(22)

Combining Eqs. (20) and (22), we can write the energy loss in terms of a
local-time recoil matrix $\mathcal{M}$:

$\displaystyle\Delta$
$\displaystyle=\int_{-\infty}^{\infty}\underbrace{d\tau\dot{\bm{\sigma}}(\tau)}_{d\bm{\sigma}}\cdot\,\mathbf{H}\Phi\left[\bm{\sigma}(\tau)\right]\mathcal{M}\mathbf{H}\Phi\left[\bm{\sigma}(\tau)\right]\cdot\dot{\bm{\sigma}}(\tau)\,,$
$\displaystyle\mathcal{M}$
$\displaystyle=\sum_{u}\int_{0}^{\pi}d\theta\sin\theta\oint
d\phi\frac{\bm{\eta}_{u,\theta,\phi}\bm{\eta}_{u,\theta,\phi}^{T}}{4\left(2\pi\sqrt{k_{u,\theta,\phi}}\right)^{3}}\,.$
(23)

By writing $d\tau\dot{\bm{\sigma}}\rightarrow d\bm{\sigma}$, we can identify
the remaining portion of the integrand as the dissipative force which is both
time-local and proportional to the particle velocity. For the cubic lattice,
numerical integration over $\theta$ and $\phi$ yields an $\mathcal{M}$ that is
diagonal with identical entries for each nonzero element, as expected from the
system symmetry.

Using this recoil matrix $\mathcal{M}$, we then compute the Hessian of the
potential and, assuming a constant speed so that
$\bm{\sigma}=\dot{\bm{\sigma}}_{0}\tau+\bm{\sigma}_{0}$, take the time
integral in Eq. (23) numerically to give the low-speed loss. The results are
plotted as straight lines in Fig. 4 for each combination of $\lambda$ and $s$,
demonstrating that, at low speeds, Eq. (23) agrees very well with the
numerical calculations from Eq. (13).

Figure 5: _Dissipation in 3D lattice with multiple interactions._ As a mobile
particle moves through the middle of a unit cell, it loses energy due to its
interaction with (a) four or (b) eight masses in a cubic lattice, as shown in
Fig. 3(c) and (d). The system parameters are the same as in Fig. 4, but with
$\Phi_{0}$ reduced to $0.0025$. The low-speed asymptotic behavior is given by
Eq. (23), showing that neglecting the phase difference in the phonon modes of
adjacent atoms is valid. $\Delta$’s nonmonotonicity with $\lambda$ arises from
a combined effect of two opposing phenomena: narrower potentials have a
greater second derivative (enhancing dissipation) but the particle spends less
time being acted on by the potential (diminishing dissipation).

Our analysis thus far has focused on a rather artificial scenario, where the
particle interacts with only a single framework mass. However, including
multiple ion couplings in Eq. (11) substantially complicates the expression by
introducing a phase prefactor to the vector $\bm{\eta}$, as discussed in Sec.
II. Fortunately, our single-interaction results give an important hint that
leads to a dramatic simplification. As we have shown, at low particle speeds,
the dissipation is dominated by low-frequency, long-wavelength framework
modes. For these modes, the phase difference between neighboring masses is
small, allowing us to neglect it and use Eq. (23) with an interaction term
$\Phi$ that includes multiple framework masses.

To validate this simplification, we performed two sets of simulations. For the
first set, we initialized the particle normally incident to one face of a
cubic unit cell, and set it to interact with the four framework masses
defining that face, as shown in Fig. 3(c). As before, the initial position of
the particle was $7\lambda$ away from the face of the unit cell and particle
was constrained to move in one dimension, passing through the middle of the
face. For the second simulation set, the particle, normally incident to one of
the unit cell faces, traversed the entire unit cell while interacting with all
eight masses at the corners of the unit cell, see the schematic in Fig. 3(d).
In both cases, the simulations terminated once the particle reached a distance
of $7\lambda$ from the face containing the interacting masses, at which point
its kinetic energy was calculated and subtracted from the initial value. The
results for the energy losses in these two scenarios are shown in Fig. 5. In
addition to the simulation results, we also computed $\Delta$ using Eq. (23)
for the appropriate form of $\Phi$, shown as horizontal lines in Fig. 5,
demonstrating that the asymptotic limit does an excellent job of capturing the
behavior of $\Delta$ at low speeds.

Through a series of increasingly realistic interaction scenarios, we have seen
that local-time dynamics of the framework in response to its interaction with
the particle describe dissipation well. Although these configurations remain
somewhat contrived since the particle moves along a high-symmetry direction,
one should keep in mind that their goal is to establish the validity of the
local-time approach, not to study realistic particle trajectories. The fact
that the local-time treatment works well regardless whether the particle
passes through the potential minimum in the middle of the unit cell or
directly through the framework mass demonstrates the robustness of the
approach. Additionally, without including thermal motion of the lattice, which
is the subject of our next work, any choice regarding the particle’s
trajectory is bound to be artificial.

Now that we have confirmed the validity of Eq. (23) at low speeds, we can
explore the implications of a regime where the force is time-local. To make
use of this simplified formula, we rewrite Eq. (3) for a small framework
deflection by expanding the interaction to leading order in $\bm{\rho}$ and by
using Eq. (9), leading to

$\displaystyle\frac{\mu}{4\pi^{2}}\ddot{\bm{\sigma}}(\tau)$
$\displaystyle=-\nabla_{\bm{\sigma}}\left[\Phi\left[\bm{0},\bm{\sigma}(\tau)\right]-\nabla_{\bm{\rho}}\Phi\left[\bm{0},\bm{\sigma}(\tau)\right]\cdot
2\pi\int^{\tau}d\tau^{\prime}\sum_{j}\bm{\varepsilon}_{j}\bm{\varepsilon}_{j}^{\dagger}\frac{\sin\left[2\pi\omega_{j}(\tau-\tau^{\prime})\right]}{\omega_{j}}\nabla_{\bm{\rho}}\Phi\left[\bm{0},\bm{\sigma}(\tau^{\prime})\right]\right]\,.$
(24)

Integrating the last term of Eq. (24) over $\tau^{\prime}$ by parts gives

$\displaystyle\int^{\tau}d\tau^{\prime}\frac{\sin\left[2\pi\omega_{j}(\tau-\tau^{\prime})\right]}{\omega_{j}}\nabla_{\bm{\rho}}\Phi[\bm{0},\bm{\sigma}(\tau^{\prime})]$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi\omega_{j}^{2}}\nabla_{\bm{\rho}}\Phi[\bm{0},\bm{\sigma}(\tau)]$
$\displaystyle-$
$\displaystyle\int^{\tau}d\tau^{\prime}\frac{\cos\left[2\pi\omega_{j}(\tau-\tau^{\prime})\right]}{2\pi\omega_{j}^{2}}\nabla_{\bm{\sigma}}\nabla_{\bm{\rho}}\Phi[\bm{0},\bm{\sigma}(\tau^{\prime})]\cdot\dot{\bm{\sigma}}(\tau^{\prime})\,.$
(25)

When reinserted back into Eq. (24), we see that the second term can be
identified as giving rise to the dissipative force due to its proportionality
to the velocity. This allows us to replace it by the term from the integrand
of Eq. (23) multiplying the velocity. The remaining portion can be combined
with the first term in Eq. (24), resulting in a softened potential so that

$\displaystyle\frac{\mu}{4\pi^{2}}\ddot{\bm{\sigma}}=$
$\displaystyle-\nabla_{\bm{\sigma}}\left[\Phi\left(\bm{\sigma}\right)-\nabla_{\bm{\rho}}\Phi\left(\bm{0},\bm{\sigma}\right)\cdot\sum_{j}\frac{\bm{\varepsilon}_{j}\bm{\varepsilon}_{j}^{\dagger}}{\omega_{j}^{2}}\nabla_{\bm{\rho}}\Phi(\bm{0},\bm{\sigma})\right]$
$\displaystyle-\mathbf{H}\Phi(\bm{\sigma})\mathcal{M}\mathbf{H}\Phi(\bm{\sigma})\cdot\dot{\bm{\sigma}}\,,$
(26)

which is the main result of this work. In this form, we can see that the
local-time interaction with the framework gives rise to two phenomena:
potential softening and dissipation. Because both of these terms involve sums
$\sim\omega_{j}^{-2}$, stiff systems (with large $\omega_{j}$) will be less
susceptible to deformation, resulting in a reduced softening and dissipation.
Qualitatively, this result matches our physical intuition.

Figure 6: _Local-time dissipation._ A particle constrained to move along one
of the edges of the unit cell with insufficient energy to escape the local
minimum for $\lambda=1/2$ and $\Phi_{0}=1$. The full solution from Eqs.
(3)–(4) is obtained using the same procedure as all other simulations. The
time-local solution is obtained by integrating Eq. (26) using the fifth order
Runge-Kutta method with $\delta\tau=10^{-3}$. The amplitude decay is similar
for the two solutions, and the phase difference is related to small
displacements of the framework that are neglected by the time-local approach.
The amplitude increase for the full solution at $\tau\approx 25$ is a
consequence of the finite size of the system. The calculation times on a
single core of a 2020 M1 MacBook Pro for the full solution and the time-local
formalism were about 4 hours and about 9 seconds, respectively.

We numerically explore this time-local behavior using the same cubic lattice
as in Figs. 4 and 5, but extended to include $100\times 100\times 100$ masses.
The mobile particle interacts with eight framework masses in a line with
equilibrium positions $\alpha(0,0,j)$ for integer $0\leq j\leq 7$ and the
lattice constant $\alpha=2$. We initialize a mobile particle halfway between
the 4th and the 5th masses with speed $\dot{\sigma}$ too small to overcome the
energy barrier. We then compute the trajectory of the particle using the full
simulation and Eq. (26), then plot the resulting position and velocity
evolution in Fig. 6. Although a more physical trajectory would pass through
the middle of the unit cell, similar to the setup in Fig. 5, it would have
resulted in a reduced energy barrier and, thus, a much smaller range of
possible initial speeds. Particles moving at these slow speeds would rapidly
dissipate their energy after only a few approaches toward the nearby masses.
In order to showcase how non-local time effects can be ignored even when the
particle repeatedly interacts with the same masses over an extended period of
time, we chose a less-physical setup.

In Fig. 6, we see that the time-local result, both for position and speed,
agrees well with the full solution. The main disparity is the phase difference
in the two sets of trajectories, which becomes more pronounced after the
particle performs several oscillations around its minimum and the framework
masses start to move. Despite this phase difference, the envelope of the decay
agrees very well between the two simulations. The mild disagreement at larger
values of $\tau$ can also be attributed, in part, to the finite size of the
system: when we reduced the system size to $80\times 80\times 80$, the
disagreement increased. Since these simulations are computationally costly, we
chose not to increase the system size beyond $100\times 100\times 100$.
Naturally, in a real system, thermal motion of the framework will prevent the
mobile ion from following such a simple path, so the large-$\tau$ effects will
be less relevant.

Up to now, we have been working using dimensionless equations. To connect our
results to real materials, we reintroduce units to Eq. (26):

$M\ddot{\mathbf{R}}=-\nabla_{\mathbf{R}}\left[U\left(\mathbf{R}\right)-\nabla_{\mathbf{r}}U\left(\bm{0},\mathbf{R}\right)\cdot\sum_{j}\frac{1}{m}\frac{\bm{\varepsilon}_{j}\bm{\varepsilon}_{j}^{\dagger}}{\Omega_{j}^{2}}\nabla_{\mathbf{r}}U(\bm{0},\mathbf{R})\right]-\mathbf{H}U(\mathbf{R})\left(\frac{2\pi}{m\Omega^{3}}\mathcal{M}\right)\mathbf{H}U(\mathbf{R})\cdot\dot{\mathbf{R}}\,,$
(27)

with

$\displaystyle\frac{2\pi}{m\Omega^{3}}\mathcal{M}$
$\displaystyle=\frac{a^{3}}{m}\frac{1}{(4\pi)^{2}}\sum_{u}\int
dS\frac{\bm{\eta}_{u,\theta,\phi}\bm{\eta}_{u,\theta,\phi}^{T}}{\left(a\Omega\sqrt{k_{u,\theta,\phi}}\right)^{3}}$
$\displaystyle=\frac{1}{\rho}\frac{1}{(4\pi)^{2}}\sum_{u}\int
dS\frac{\bm{\eta}_{u,\theta,\phi}\bm{\eta}_{u,\theta,\phi}^{T}}{v_{u,\theta,\phi}^{3}}\,,$
(28)
$\displaystyle\sum_{j}\frac{1}{m}\frac{\bm{\varepsilon}_{j}\bm{\varepsilon}_{j}^{\dagger}}{\Omega_{j}^{2}}$
$\displaystyle=V^{-1}\sum_{j}V\frac{\bm{\varepsilon}_{j}\bm{\varepsilon}_{j}^{\dagger}}{m\Omega_{j}^{2}}=V^{-1}\,.$
(29)

Here the integral is taken over the solid angle, $a$ is the lattice constant,
$v_{u,\theta,\phi}$ is the direction-dependent speed of the acoustic branch
$u$, and $\rho=m/a^{3}$ is the framework density.

To make the physical significance of the softening term clearer, we start with
the equation of motion for the framework
$m\ddot{\mathbf{r}}=-V\mathbf{r}-\nabla_{\mathbf{r}}U(\mathbf{r},\mathbf{R})$.
Setting the left-hand side to zero, we obtain the equation which determines
$\mathbf{r}$ for which the interaction with the particle is balanced by the
elastic forces. For small framework deformation, we get
$\mathbf{r}\approx-V^{-1}\nabla_{\mathbf{r}}U(\mathbf{0},\mathbf{R})$. Thus,
the softened potential in Eq. (27) is the leading order Taylor expansion of
the relaxed configuration.

##  V Discussion

In the course of our analysis, we have employed a number of simplifications to
make the problem tractable. Having validated these simplifications
numerically, let us now address their implications in the context of more
general systems.

Our use of the Gaussian interaction was driven primarily by convenience. At
first glance, this form might appear highly artificial because it allows, in
principle, for two objects to have zero separation–an unphysical situation for
real systems where two ions cannot occupy the same position. Nevertheless,
this choice is less contrived than it might appear as combinations of Gaussian
functions are commonly used to describe orbitals in molecular simulations
[26]. By using a single Gaussian function, we treat our ions as spherically
symmetric and do not consider “lobes” in their associated electronic clouds
that can arise due to bond formation. Including multiple Gaussian functions in
$\Phi$ makes the problem more complex analytically (e.g., by having Eq. (13)
contain a sum of Fourier transforms), but does not change the approach or the
main conclusions.

The cubic geometry was also chosen for its simplicity. Although real systems
generally have a more complex structure, the fact that our approach works well
for different dimensionalities supports its robustness. In fact, the only set
of simulations that explicitly reflect the cubic geometry are the ones given
in Fig. 5(b): Fig. 4 describes the particle’s interaction with a single
framework mass, while the square described by the four masses in Fig. 5(a)
could belong to a face of a monoclinic unit cell. Crucially, the local-time
formula makes no reference to the lattice geometry as the dissipation matrix
$\mathcal{M}$ depends on the speed of sound and material density, while the
interaction terms are pairwise and fairly short-range, preventing the mobile
particle from “seeing” the long-range order of the lattice.

Finally, although our study focuses on a single mobile particle, its results
are not limited to the case of a single current-carrying ion in an electrolyte
system. Instead, we can think of the model as describing a scenario where the
mobile ions are dilute enough so that their interaction can be neglected. When
performing nudged elastic bands calculations, supercells containing a few unit
cells in each direction are required to achieve convergence. Consequently,
even separations of a few unit cells between the mobile ions are sufficient to
view them as independent. The role of concerted motion or blocking merits a
careful investigation in follow-up work.

##  VI Summary

Starting with a simplified 1D setup, in this work we demonstrated that for
sufficiently slow-moving particles, the dissipation in solids is dominated by
low-frequency, long-wavelength modes. Using this result, we derived a time-
local formula for dissipation in 3D solids, where the dissipative term is
proportional to the speed of the moving particle, reminiscent of drag in
liquids. We showed that the drag depends on the curvature of the potential
landscape through which the mobile particle moves, as well as the stiffness of
the framework lattice. We validated this newly-derived formula by comparing
its results to trajectories computed using the full equations of motion.
Although our formula makes it possible to calculate the trajectory of an
individual mobile particle, its chief utility lies in its ability to elucidate
the key features of a good ionic electrolyte by doing away with the complexity
associated with time nonlocality. Specifically, because the potential
landscape and speed of sound can be computed using _ab initio_ methods, our
formula can be used as a screening tool for potential electrolyte materials.

A natural extension of this work would investigate the role of thermal
fluctuations in the time-local regime. In this work, we assumed that
interaction-induced displacements of the framework are small, allowing us to
keep only the first-order correction when calculating the interactions. In a
framework dominated by thermal vibrations, we suspect that neglecting higher-
order contributions to the interaction-driven displacements will be an even
better approximation. Thermal vibrations will effectively reduce the memory of
the system, further justifying the local-time approach, but for a more
realistic class of systems.

###### Acknowledgements.

A.R. acknowledges the National Research Foundation, Prime Minister Office,
Singapore, under its Medium Sized Centre Programme and the support by Yale-NUS
College (through Start-up Grant). B. A. O. acknowledges support from the M. J.
Murdock Charitable Trust. A.U. was supported by the Ministry of Education,
Singapore, under its funding for the Research Centre of Excellence Institute
for Functional Intelligence Materials and by the National Research Foundation,
Singapore under its AI Singapore Programme (AISG Award No: AISG3-RP-2022-028).

## 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).
  * 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).
  * Marcolongo and Marzari [2017] A. Marcolongo and N. Marzari, Ionic correlations and failure of nernst-einstein relation in solid-state electrolytes, Physical Review Materials 1, 1 (2017).
  * He  _et al._ [2017] X. He, Y. Zhu, and Y. Mo, Origin of fast ion diffusion in super-ionic conductors, Nature Communications 8, 1 (2017).
  * Stefano  _et al._ [2019] D. D. Stefano, A. Miglio, K. Robeyns, Y. Filinchuk, M. Lechartier, A. Senyshyn, H. Ishida, S. Spannenberger, D. Prutsch, S. Lunghammer, D. Rettenwander, M. Wilkening, B. Roling, Y. Kato, and G. Hautier, Superionic diffusion through frustrated energy landscape, Chem 5, 2450 (2019).
  * 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. 31, 8850 (2019).
  * Poletayev  _et al._ [2022] A. D. Poletayev, J. A. Dawson, M. S. Islam, and A. M. Lindenberg, Defect-driven anomalous transport in fast-ion conducting solid electrolytes, Nat. Mater 21, 1066 (2022).
  * Carvalho  _et al._ [2022] A. Carvalho, S. Negi, and A. H. C. Neto, Direct calculation of the ionic mobility in superionic conductors, Sci. Rep. 12, 19930 (2022).
  * 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).
  * Wang  _et al._ [2015] Y. Wang, W. D. Richards, S. P. Ong, L. J. Miara, J. C. Kim, Y. Mo, and G. Ceder, Design principles for solid-state lithium superionic conductors, Nature Materials 14, 1026 (2015).
  * Bo  _et al._ [2016] S. H. Bo, Y. Wang, J. C. Kim, W. D. Richards, and G. Ceder, Computational and experimental investigations of na-ion conduction in cubic na3pse4, Chemistry of Materials 28, 252 (2016).
  * 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).
  * Brenner  _et al._ [2020] T. M. Brenner, C. Gehrmann, R. Korobko, T. Livneh, D. A. Egger, and O. Yaffe, Anharmonic host-lattice dynamics enable fast ion conduction in superionic agi, Physical Review Materials 4, 115402 (2020).
  * Rodin  _et al._ [2022a] A. Rodin, K. Noori, A. Carvalho, and A. H. C. Neto, Microscopic theory of ionic motion in solids, Phys. Rev. B 105, 224310 (2022a).
  * 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).
  * 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, 65 (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._ [2022b] 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 (2022b).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Saunders College Publishing, 1976) p. 826.
  * Cramer [2004] C. J. Cramer,  _Essentials of Computational Chemistry, Theories and Models_, 2nd ed. (Wiley, 2004).


# Graphene valley filter using a line defect

D. Gunlycke  Naval Research Laboratory, Washington, D.C. 20375, USA  C. T.
White  Naval Research Laboratory, Washington, D.C. 20375, USA

###### Abstract

With its two degenerate valleys at the Fermi level, the band structure of
graphene provides the opportunity to develop unconventional electronic
applications. Herein, we show that electron and hole quasiparticles in
graphene can be filtered according to which valley they occupy without the
need to introduce confinement. The proposed valley filter is based on
scattering off a recently observed line defect in graphene. Quantum transport
calculations show that the line defect is semitransparent and that
quasiparticles arriving at the line defect with a high angle of incidence are
transmitted with a valley polarization near 100%.

###### pacs:

73.22.Pr, 73.61.Wp, 73.63.Bd, 85.75.-d

Owing to its exceptional electronBolo08_1  and thermalBala08_1  transport
properties, grapheneNovo04_1  is a promising material for use in advanced low-
energy electronic applications. As graphene is a semimetal with no band
gap,Wall47_1  it cannot in its intrinsic form be used as a replacement
material for silicon and other semiconductors in conventional electronics.
However, with its bands at the Fermi level defining two conical valleys,
graphene might instead offer novel electronic applications. Before such
applications can be realized, the control over transport in graphene needs to
be improved. Understanding the effects a recently observed line defectLahi10
and grain boundariesHuan11  have on the transport properties might be the key.
While there are several theoretical studies that have investigated the latter
grain boundaries,Malo10 ; Yazy10 ; Liu10 ; Yazy10a ; Huan11  the line defect
has so far remained relatively unexplored, even though the line defect is
always straight and its adjoining grains are aligned, arguably making the line
defect a more suitable structure for controlled transport in graphene.

This Letter investigates the transport properties of the atomically precise,
self-assembled graphene line defect shown in Fig. 1 and shows that this line
defect is semitransparent and can be used as a valley filter.

Figure 1: Extended line defect in graphene. (a) Arrangement of the carbon
atoms around the line defect highlighted in gray. The structure exhibits
translational symmetry along the defect with a primitive cell shown in beige.
The two sublattices in graphene indicated by blue and green atoms reverse upon
reflection at the line defect, as can be clearly seen in the enlarged overlay.
(b) Scanning tunneling microscope image of the line defect in graphene (Image
adapted from Ref. [Lahi10, ]).

Electron and hole quasiparticles can either transmit through the line defect
without changing direction or reflect following the law of specular
reflection. The transmission and reflection probabilities depend on the valley
degree of freedom, thus allowing the quasiparticles to be filtered according
to their valley degree of freedom with a polarization near 100% for
quasiparticles arriving at the line defect with a high angle of incidence.
This filter is much different from the valley filter proposed by Rycerz et
al.,Ryce07  which relies on the isolation of a few one-dimensional channels in
a narrowly confined region. The latter filter has not yet been demonstrated,
presumably due to challenges in fabricating the structure, which requires a
sub-10 nm constriction with saturated zigzag edges. In contrast, the filter
proposed herein relies on the two-dimensional geometry of graphene and its
required structure has already been observed,Lahi10  as demonstrated by the
micrograph in Fig. 1. The valley filter is expected to be a central component
in valleytronics,Ryce07  just as the spin filter is central in spintronics.
Electronics that make use of the two valleys in graphene is attractive because
the valleys are separated by a large wave vector, making valley information
robust against scattering from slowly varying potentials,Ando98  including
scattering caused by intravalley acoustic phonons that often limit coherent
low-bias devices to low-temperature operation. Therefore, the valley
information generated by the filter proposed herein could in principle be
preserved even in a diffusive charge transport regime.

Figure 2: Valley states scattering off the line defect. (a) Energies close to
the Fermi level in graphene have dark contours and are located in the corners
of the first Brillouin zone contained within the gold hexagon. The two
valleys, K and K′, are identified as the two disjoint low-energy regions in
the reciprocal primitive cell enclosed by the blue rectangle. (b) An incident
quasiparticle state is defined by the valley index $\tau$ and wave vector
$\vec{q}$, where the latter points in the direction $\hat{q}$ given by the
angle of incidence $\alpha$. (c) Owing to energy and momentum conservation
along the line defect, there are only two nonevanescent scattered states
allowed. (d) The sublattice symmetric $|+\rangle$ and antisymmetric
$|-\rangle$ components of the incident state $|\Phi_{\tau}\rangle$ are
transmitted and reflected, respectively. The thickness of each arrow indicates
the probability the quasiparticle will follow the respective path.

Consider a low-energy electron (hole) quasiparticle with energy $\varepsilon$
and valley index $\tau$ approaching the line defect from the left (right) at
the angle of incidence $\alpha$. Asymptotically far from the line defect, the
quasiparticle occupies a graphene state $|\Phi_{\tau}\rangle$, where $\tau=\pm
1$ is a valley index. See Fig. 2. Let the quasiparticle wave vector
$\vec{q}=(q_{x},q_{y})$ be the wave vector measured from the center of the
occupied valley, located at $\vec{K}_{\tau}=4\pi\tau\hat{y}/3a$, where $a$ is
the graphene lattice constant and $\hat{y}$ is the unit vector along the line
defect. To first order in $q\equiv|\vec{q}|$, the nearest-neighbor tight-
binding Hamiltonian in grapheneWall47_1  can be expressedSlon57  as
$H_{\tau}=\hbar v_{F}\left(q_{x}\sigma_{y}+\tau q_{y}\sigma_{x}\right)$, where
$v_{F}=\sqrt{3}|\gamma|a/2\hbar$ is the Fermi velocity with the nearest-
neighbor hopping parameter $\gamma\approx-2.6$ eV, and $\sigma_{x}$ and
$\sigma_{y}$ are Pauli matrices. The Hamiltonian has energy eigenvalues
$E=\eta\varepsilon$, where $\varepsilon=\hbar v_{F}q$ and $\eta=+1$ ($-1$) if
the quasiparticle is an electron (hole). From the quasiparticle energy
dispersion $\varepsilon$, it follows that the quasiparticle group velocity is
$v_{F}\hat{q}$, where $\hat{q}$ is the unit vector in the direction $\vec{q}$.
Because the quasiparticle travels in the direction $\vec{q}$, $q_{x}=\eta
q\cos\alpha$ and $q_{y}=\eta q\sin\alpha$. Using these relations, the
eigenstate of the graphene Hamiltonian for a given $\tau$ and $\alpha$ can be
expressed as

$|\Phi_{\tau}\rangle=\frac{1}{\sqrt{2}}\left(|A\rangle+ie^{-i\tau\alpha}|B\rangle\right),$
(1)

where $|A\rangle$ and $|B\rangle$ refer to the two sublattices in graphene.

The structure in Fig. 1 exhibits a useful symmetry line through the line
defect. In the limit $q\rightarrow 0$, the reflection operator commutes with
the graphene translation operator perpendicular to the line defect. Therefore,
symmetry-adapted states $|\pm\rangle$ can be constructed that are simultaneous
eigenstates of the graphene Hamiltonian and the reflection operator. As the
reflection operator maps $A$ sites onto $B$ sites, and vice versa, it can be
represented by the operator $\sigma_{x}$ acting on the two sublattices. From
the eigenstates of $\sigma_{x}$, one obtains

$|\pm\rangle=\frac{1}{\sqrt{2}}\left(|A\rangle\pm|B\rangle\right).$ (2)

The graphene state (1) expressed in the symmetry-adapted basis is

$|\Phi_{\tau}\rangle=\frac{1+ie^{-i\tau\alpha}}{2}|+\rangle+\frac{1-ie^{-i\tau\alpha}}{2}|-\rangle.$
(3)

The full Hamiltonian describing the system can be divided into three terms,
$\mathcal{H}_{\tau}=H_{\tau}+H_{D}+V$, which represent graphene, the isolated
line defect, and the interaction between graphene and the line defect,
respectively. As each term commutes with the reflection operator, the full
Hamiltonian must commute with the reflection operator, and thus, the
eigenstates of $\mathcal{H}_{\tau}$ in the symmetry-adapted basis are either
symmetric or antisymmetric about the line defect. Antisymmetric states have a
node at the line defect, and as a result, there are no matrix elements within
the nearest-neighbor model coupling the left and right sides. Therefore,
antisymmetric states cannot contribute to any transmission across the line
defect. As shown below, however, there are two symmetric states at the Fermi
level without a node on the line defect. As these states are extended
eigenstates of the full Hamiltonian, they carry quasiparticles across the line
defect without scattering. Thus, we can conclude that the transmission
probability of the quasiparticle approaching the line defect is

$T_{\tau}=\left|\langle+|\Phi_{\tau}\rangle\right|^{2}=\frac{1}{2}\left(1+\tau\sin\alpha\right).$
(4)

As the sum of the transmission probabilities in Eq. (4) over the two valleys
is $\sum_{\tau}T_{\tau}=1$, we can also conclude that the line defect is
semitransparent. The semitransparency follows from the relation
$\langle+|\Phi_{-\tau}\rangle=\langle-|\Phi_{\tau}\rangle^{*}$ and the
normalization of $|\Phi_{\tau}\rangle$. See Fig. 2. Fig. 3 shows that the
transmission probability of a quasiparticle varies significantly with its
angle of incidence $\alpha$.

Figure 3: The probability that an incident quasiparticle at the Fermi level
with valley index $\tau$ and angle of incidence $\alpha$ will transmit through
the line defect.

At a high angle of incidence, there is almost full transmission or reflection,
depending on the valley index $\tau$.

Owing to the semitransparency, given an unpolarized beam of incident
quasiparticles, the probability $P_{\tau}$ that a transmitted quasiparticle
has valley index $\tau$ is given by $T_{\tau}$. The polarization
$\mathcal{P}\equiv\langle\tau\rangle=P_{+1}-P_{-1}$ of the transmitted beam is

$\mathcal{P}=\sin\alpha.$ (5)

This expression shows that an unpolarized beam of quasiparticles approaching
the line defect at a high angle of incidence will lead to outgoing transmitted
and reflected beams that are almost completely polarized.

The results presented above are based on symmetry arguments that are valid
only in the limit as the quasiparticle energy $\varepsilon\rightarrow 0$. The
results, however, hold to an excellent approximation as long as
$\varepsilon\ll\hbar v_{F}/a\approx 2.3$ eV. To show this, we have performed
numerical transport calculations that are treated exactly within the nearest-
neighbor tight-binding model. The transmission probability is shown in Fig. 4
as a function of the components of the wave vector $\vec{q}$.

Figure 4: Transmission probability for a quasiparticle with finite energy. (a)
The transmission probability for a quasiparticle approaching the line defect
at an angle of incidence $\alpha$ is indicated by the brightness, ranging from
0 to 1 in increments of 0.05. For the angle $\alpha$ shown, the transmission
is near one if the quasiparticle has valley index $\tau=+1$, as illustrated on
the right, where the filled (open) circle represents an electron (hole)
quasiparticle. (b) The transmission probability for the corresponding
quasiparticle with valley index $\tau=-1$.

Tracing the transmission probability in the figure along any constant energy
contour, which is almost perfectly circular due to the approximately conic
dispersion, yields a dependence on the angle of incidence that is virtually
indistinguishable from that in Fig. 3, thus confirming that almost full
polarization can be achieved as long as $qa\ll 1$ or $\varepsilon\ll\hbar
v_{F}/a$. To further test the robustness of the valley filter, we have
performed calculations with interactions across the line defect and with
potentials on the line defect and their neighboring sites. None of these tests
led to results qualitatively different from those presented in Fig. 4.

Because the structure in Fig. 1(a) exhibits a symmetry plane through the
center of the primitive cell, there is time-reversal symmetry in the direction
along the line defect. This time-reversal symmetry implies that the
transmission probability of a quasiparticle with valley index $-\tau$ can be
obtained from that of a quasiparticle with valley index $\tau$ by letting
$q_{y}\rightarrow-q_{y}$. This relationship between $\tau=\pm 1$ can be seen
in Fig. 4. Note, however, that the transmission probabilities in Fig. 4 are
not symmetric about $q_{y}$, and thus one can conclude that the scattering of
a quasiparticle depends on the valley index, which is a necessary requirement
for a valley filter. As both graphene and the line defect have Hamiltonians
exhibiting electron–hole symmetry, one might expect that
$T_{\tau}(\vec{q})=T_{\tau}(-\vec{q})$. The scattering does, in general
however, not obey electron–hole symmetry. That the condition is satisfied by
Eq. (4) is in part a consequence of the neglect of evanescent waves. These
evanescent waves are accounted for in the numerically obtained transmission
probability in Fig. 4. Note, for instance, that
$T_{\tau}(0.005\pi/a,0.01\pi/a)\neq T_{\tau}(-0.005\pi/a,-0.01\pi/a)$. To
understand the lack of electron–hole symmetry in the combined graphene–line
defect system, we note that the structure in Fig. 1(a) is not bipartite; in
particular, the sites participating in the pentagons at the line defect cannot
be divided into two types where one type has only nearest neighbors of the
other type.

To gain further insight into the numerical calculations and how they lead to
Eq. (4), it is useful to perform the transport calculations analytically in
the limit of small $q$. As there is translational symmetry along the line
defect, one can focus on those atoms within the primitive cell shown in Fig.
1. The Hamiltonian of the isolated line defect is then

$H_{D}=\gamma\left(\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right).$ (6)

Next, we seek a retarded self energy $\Sigma$ that accounts for the coupling
of the two line defect atoms to the semi-infinite portion of graphene on each
side of the line defect. As shown in Fig. 1, there are two atoms neighboring
the line defect on each side. Expressed in the basis of the atoms parallel to
the line defect, the graphene state with valley index $\tau$ is given by

$|\tau\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\ e^{-2\pi
i\tau/3}\end{array}\right).$ (7)

As the basis contains two atoms belonging to the same sublattice, the graphene
state above can be folded onto another graphene state with the same wave
vector in the full system. This latter state, however, is evanescent near the
Fermi level and can be neglected. Requiring $\Sigma$ to be retarded fixes the
relative phase between the atoms on the line defect and their neighbors,
resulting in the relation $\Sigma\langle B|\Phi_{\tau}\rangle=\gamma\langle
A|\Phi_{\tau}\rangle|\tau\rangle\langle\tau|$, from which one obtains

$\Sigma=-\frac{i\gamma}{2}e^{i\tau\alpha}\left(\begin{array}[]{cc}1&e^{2\pi
i\tau/3}\\\ e^{-2\pi i\tau/3}&1\end{array}\right).$ (8)

Equipped with $H_{D}$ describing the interactions within the line defect and
$\Sigma$ describing the coupling to the semi-infinite graphene on each side,
one can calculate the retarded Green function on the line defect,
$G=\left(\eta\varepsilon I-H_{D}-2\Sigma\right)^{-1}$, where $I$ is the unit
matrix. To zeroth order in $q$,

$G=\frac{-\gamma^{-1}}{1+ie^{i\tau\alpha}}\left(\begin{array}[]{cc}ie^{i\tau\alpha}&1-ie^{i\tau(\alpha+2\pi/3)}\\\
1-ie^{i\tau(\alpha-2\pi/3)}&ie^{i\tau\alpha}\end{array}\right).$ (9)

The probability that the quasiparticle will transmit through the line defect
is given by $T_{\tau}=\langle\tau|\Gamma G\Gamma G^{\dagger}|\tau\rangle$,
where $\Gamma\equiv i\left(\Sigma-\Sigma^{\dagger}\right)$. Inserting Eqs.
(7–9) into this equation, one recovers Eq. (4) exactly.

In the initial analysis leading to Eq. (4), an assertion was made that there
are two symmetric states at the Fermi level without a node on the line defect.
That claim can now be verified using Eqs. (1, 6, 8). According to Eq. (1), a
symmetric state must satisfy $\tau\alpha=\pi/2$. When this condition is
satisfied, one finds that the determinant
$\operatorname{det}\left(H_{D}+2\Sigma\right)=0$, which implies that there are
exactly two symmetric states at the Fermi level, one for each valley index
$\tau$. The corresponding eigenstates are $|\Psi_{\tau}\rangle=|-\tau\rangle$,
confirming that the states have no node at the line defect.

Rather than forming isolated Bloch waves, experiments exploiting the valley
filter will likely construct more complex wave patterns. The dimensions of the
system should be chosen such that the mean-free path is longer than the
distance between the source and the line defect. For low-energy
quasiparticles, the wavelength $\lambda=2\pi/q$ is much greater than the
repeating length $2a$ of the line defect. As long as this repeating length is
also much shorter than any spatial features of the waves, to an excellent
approximation the scattering can be treated within ray optics, where rays
travel in straight lines and only scatter at the line defect, where they will
either transmit with a probability approximately given by Eq. (4) or reflect
while obeying the law of specular reflection.

The filter can be used to create a valley-polarized beam of electrons or
holes. By probing the current passing through the line defect at a particular
angle, one could also measure the valley polarization of the incident
quasiparticles. Demonstration of these components should significantly
accelerate research on graphene valleytronics.

###### Acknowledgements.

The authors acknowledge support from the U.S. Office of Naval Research,
directly and through the U.S. Naval Research Laboratory. D.G. thanks C.W.J.
Beenakker for helpful comments.

## References

  * (1) K. I. Bolotin, K. J. Sikes, Z. Jiang, G. Fudenberg, J. Hone, P. Kim, and H. L. Stormer, Solid State Commun. 146, 351 (2008)
  * (2) A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, and C. N. Lau, Nano Lett. 8, 902 (2008)
  * (3) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004)
  * (4) P. R. Wallace, Phys. Rev. 71, 622 (1947)
  * (5) J. Lahiri, Y. Lin, P. Bozkurt, I. I. Oleynik, and M. Batzill, Nat. Nanotech. 5, 326 (2010)
  * (6) P. Y. Huang, C. S. Ruiz-Vargas, A. M. van der Zande, W. S. Whitney, M. P. Levendorf, J. W. Kevek, S. Garg, J. S. Alden, C. J. Hustedt, Y. Zhu, J. Park, P. L. McEuen, and D. A. Muller, Nature 469, 389 (2011)
  * (7) S. Malola, H. Häkkinen, and P. Koskinen, Phys. Rev. B 81, 165447 (2010)
  * (8) O. V. Yazyev and S. G. Louie, Phys. Rev. B 81, 195420 (2010)
  * (9) Y. Liu and B. I. Yakobson, Nano Lett. 10, 2178 (2010)
  * (10) O. V. Yazyev and S. G. Louie, Nat. Mat. 9, 806 (2010)
  * (11) A. Rycerz, J. Tworzydło, and C. W. J. Beenakker, Nat. Phys. 3, 172 (2007)
  * (12) T. Ando and T. Nakanishi, J. Phys. Soc. Jpn. 67, 1704 (1998)
  * (13) J. C. Slonczewski and P. R. Weiss, Phys. Rev. 109, 272 (1958)