# Dynamic Disorder in Supported Metal Nanocatalysts

Nanoparticle materials are ubiquitous in heterogeneous catalytic processes
and there is broad interest in their physical and chemical properties at the
atomic level.
Global experimental probes such as EXAFS, XANES and XPS, however, reveal only
their ensemble
characteristics, obscuring details of their internal structure.
My research [F. D. Vila
*et al.*, Phys. Rev. B **78**, 121404(R) (2008)]
focuses on the nanoscale theoretical modelling of the structure and x-ray
response of these systems. It provides a deep understanding of the
intra-particle heterogeneity of these systems, and their changes under
realistic conditions. [F. D. Vila
*et al.*, J. Phys. Chem. C **117**, 12446 (2013)]
My DFT/MD simulations reveal (see video on right) that both structure and charge
are inhomogeneous, and fluctuate dynamically over several time-scales. These
range from fast (200-400 fs) bond vibrations to slow fluxional bond breaking
(>10 ps). More interestingly, the equilibrium points of the vibrations are not
static, but rather driven by the stochastic motion of the center of mass over 1-4 ps
time-scales. This results in “dynamic structural disorder” (DSD), semi-melted
cluster surfaces, atomic-segregation, and modulation of adsorbate
dynamics. [J. J. Rehr
and F. D. Vila, J. Chem. Phys. **140**,
134701 (2014)] These fluctuations also affect the catalytic properties of
the nanoparticles by changing the reaction landscape (Figure 1).

Figure 1: Reaction paths for the dissociation of
O_{2} on Pt_{10}Sn_{10} as a function of nanoparticle
disorder. These reaction paths show different behaviors, from exothermic to
endothermic thermodynamics. The role of disorder is clear in path 1, where
the incipient formation of a Sn-O bond stabilizes the transition state.

# *Ab Initio* Debye-Waller Factors

Thermal and structural disorder effects are of crucial importance in X-ray
spectroscopy since they lead to strong exponential damping of the fine
structure. This damping is dominated by the XAFS Debye-Waller factor
\( W = \mathrm{exp}\left(-2\sigma^{2}_{R} k^2 \right)\)
where \( \sigma^{2}_{R} \)
is the mean square relative displacement (MSRD) for path \(R\).
The MSRD is defined by the thermal average
$$
\sigma^{2}_{R} = \left< \left( R - \bar{R} \right)^2 \right>,
$$
where \( R \) and \( \bar{R} \) are the instantaneous and mean path
lengths, respectively.
In the quasi-harmonic approximation, the average can be computed as a sum over
vibrational or phonon eigenmodes or, alternatively, in terms of the projected
vibrational density of states \( \rho_R \left( \omega \right) \) or VDOS, and
the Debye integral:
$$
\sigma^{2}_{R} = \frac{\hbar}{2\mu_R} \int_{0}^{\infty} \frac{1}{\omega}
\mathrm{coth} \left( \frac{\beta \hbar \omega}{2} \right) \rho_R \left(
\omega \right) d\omega,
$$
where \( \mu_R \) is the reduced mass for path \( R \).
Following on work by Poiarkova and Rehr, I have implemented
[F. D. Vila
*et al.*, Phys. Rev. B **76**, 014301 (2007)]
an efficient, *ab initio* way to compute the PDOS from the imaginary
part of lattice dynamical Green's function
$$
\rho_R \left( \omega \right) = -\frac{2\omega}{\pi} \mathrm{Im}
\left< 0 \left| \frac{1}{\omega^2 - \mathbf{D} + i\epsilon} \right| 0 \right>
\simeq \sum_{\nu=1}^{N} w_{\nu} ~ \delta \left(\omega - \omega_{\nu} \right).
$$
Here \( \left| 0 \right> \) is a vector representing a normalized mass-weighted
displacement of the atoms along the multiple-scattering path \( R \). The
lattice dynamical Green's function is represented in the form of
a continued fraction obtained from the iterative Lanczos
algorithm. This yields a \(N\)-pole representation for the VDOS where
\( \omega_{\nu} \) and \( w_{\nu} \)
are, respectively, the pole frequencies and weights of pole \( \nu \).
The dynamical matrix of force constants or Hessian, \(\mathbf{D} \), is the
main ingredient required for the calculation of theoretical MSRDs and is
defined as
$$
D_{jl\alpha,j'l'\beta} = \left( M_j M_{j'} \right)^{-1/2}
\frac{\partial^2 E}{\partial u_{jl\alpha} \partial u_{j'l'\beta}},
$$
where \( u_{jl\alpha} \) and \( M_j \) are the \( \alpha = \left\{x,y,z\right\} \)
Cartesian displacement and mass of atom \(j\) in unit cell \(l\),
respectively, and \(E\) is the internal energy of the
system.
I have implemented this method in the DMDW module of FEFF.
For molecular systems, \( \mathbf{D} \) can be computed and imported from
quantum chemistry packages like
Gaussian,
NWChem,
Orca, etc.
For solids,
the dynamical matrix can be obtained using
ABINIT,
VASP and
Quantum ESPRESSO.
Besides the MSRD, the DMDW module can be used to compute thermal expansion
coefficients and vibrational free energies.

Figure 2: Left: Comparison of the DMDW, Correlated Debye, and experimental near-neighbor MSRD for Cu and Ge. Right: Comparison of the DMDW and experimental thermal expansion for GaAs.

# Real-Time Electron Dynamics

In the past decade, Time-Dpendent Density Funcitonal Theory (TDDFT) has become
a versatile method for computing electronic response. Real time (RT) variants
of this method are particularly useful in the linear and non-linear regimes
when arbitrary electromagnetic perturbations such as pulses are used.
I am part of the development of
RT-Siesta,[Y. Takimoto
*et al.*, J. Chem. Phys. **127**, 154114 (2007)]
an efficient implementation of RT-TDDFT based on the Siesta LCAO approach.
RT-Siesta solves the time-dependent Kohn-Sham equation by time-evolving the
Kohn-Sham states in the presence of an external electric field:

$$ i \frac{\partial\Psi}{\partial t} = H\left(t\right) \Psi \qquad \longrightarrow \qquad \Psi\left(t\right) = T \mathrm{exp}\left( -i \int_{0}^{t} H\left(t'\right) dt' \right) \Psi\left(0\right) \\ H = -\frac{1}{2}\nabla^2 + V_{ext}\left(\mathbf{r},t\right) + V_H\left[\rho\right]\left(\mathbf{r},t\right) + V_{xc}\left[\rho\right]\left(\mathbf{r},t\right) $$

The optical properties are determined from the total dipole moment:

$$ \mathbf{p}\left(t\right) = \int \rho\left(\mathbf{r},t\right)\mathbf{r} d^3\mathbf{r} \\ \underbrace{\chi_{ij}^{(1)}\left(\omega\right)=\delta p_i\left(\omega\right)/E_j\left(\omega\right)=\alpha_{ij}\left(\omega\right)}_\text{Linear Response} \qquad \underbrace{\sigma\left(\omega\right) \sim \omega ~ \mathrm{Im} \left<\alpha\left(\omega\right)\right>}_\text{Absorption} $$

In RT-Siesta the states are represented in an LCAO basis, thus we just propagate the coefficients \(c\):

$$ i \frac{\partial c\left(t\right)}{\partial t} = S^{-1}H\left(t\right) c\left(t\right) $$

where \(S\) is the basis set overlap matrix. The evolution is done with a Crank-Nicolson propagator:

$$ c\left(t+\Delta t\right) = \frac{1-iS^{-1}H\left(t+\Delta t/2\right)\Delta t/2}{1+iS^{-1}H\left(t+\Delta t/2\right)\Delta t/2}c\left(t\right) + \mathcal{O} \left(\Delta t^2 \right) $$

I have implemented a predictor-corrector scheme for the propagator that ensures time-reversibility, thus enabling very long time steps. The method has been applied to the simulation of full optical spectra in complex chromophores, and the computation of their NLO response over a range of frequencies (Figure 3). Recent developments have extended RT-SIESTA to compute core spectroscopies such as XAS, XES, and response to core perturbations.

Figure 3: NLO optical response of para-nitroaniline.