Spin-Diffusions and Diffusive Molecular Dynamics

TL;DR

This paper connects diffusive molecular dynamics to a new spin-diffusion stochastic process, providing a theoretical foundation and exploring model relationships through assumptions, approximations, and test problems.

Contribution

It formulates a spin-diffusion process and demonstrates its connection to diffusive molecular dynamics, clarifying the model's theoretical basis.

Findings

Spin-diffusion can lead to diffusive molecular dynamics models.

Different assumptions yield various related models.

The models are compared through a simple test problem.

Abstract

Metastable condensed matter typically fluctuates about local energy minima at the femtosecond time scale before transitioning between local minima after nanoseconds or microseconds. This vast scale separation limits the applicability of classical molecular dynamics methods and has spurned the development of a host of approximate algorithms. One recently proposed method is diffusive molecular dynamics which aims to integrate a system of ordinary differential equations describing the likelihood of occupancy by one of two species, in the case of a binary alloy, while quasistatically evolving the locations of the atoms. While diffusive molecular dynamics has shown to be efficient and provide agreement with observations, it is fundamentally a model, with unclear connections to classical molecular dynamics. In this work, we formulate a spin-diffusion stochastic process and show how it can be connected to diffusive molecular dynamics. The spin-diffusion model couples a classical overdamped Langevin equation to a kinetic Monte Carlo model for exchange amongst the species of a binary alloy. Under suitable assumptions and approximations, spin-diffusion can be shown to lead to diffusive molecular molecular dynamics type models. The key assumptions and approximations include a well defined time scale separation, a choice of spin exchange rates, a low temperature approximation, and a mean field type approximation. We derive several models from different assumptions and show their relationship to diffusive molecular dynamics. Differences and similarities amongst the models are explored in a simple test problem.

Tables

Table 1. Table 1: Effective mobility for free energy gradient dynamics with β = 160 𝛽 160 \beta=160 .

$N$	$\frac{1}{2}\beta {\tau }^{-1}$ scaling	32	500	1000	2000
$m_c$	80	52	48	47	46
$m_r$	80	37	36	36	36

Table 2. Table 2: Values of the interaction coefficients in ( 2.8 ) for an ABABAB… composition with displacements computed at zero temperature.

Nearest Neighbor	$r_{ij}$	$J_{ij}$
First	2.28	0.0026
Second	4.56	0.0036
Third	6.84	0.0009
Fourth	9.12	0.0001

Equations

$$(\mathbf{X},\mathbf{k})\in \mathrm{argmin}\mathcal{F}(\mathbf{s},\mathbf{X},\mathbf{k}).$$

$$\frac{d{s}_i}{dt}=\sum _{j\in {\mathcal{N}}_i}{m}_{ij}\left(\frac{\partial \mathcal{F}}{\partial s_j}-\frac{\partial \mathcal{F}}{\partial s_i}\right).$$

$$d\mathbf{x}(t)=-{\nabla }_{x}V(\mathbf{x}(t),\mathbf{\sigma })dt+\sqrt{2{\beta }^{-1}}d\mathbf{W}(t),$$

$$r(\mathbf{\sigma }\to {\mathbf{\sigma }}^{\prime };\mathbf{x})e^{-\beta V(\mathbf{x},\mathbf{\sigma })}=r({\mathbf{\sigma }}^{\prime }\to \mathbf{\sigma };\mathbf{x})e^{-\beta V(\mathbf{x},{\mathbf{\sigma }}^{\prime })}.$$

$$V(\mathbf{x},\mathbf{\sigma })=\sum _{i<j}{\varphi }_{{\sigma }_i,{\sigma }_j}(|x_i-x_j|).$$

$$\mu (d\mathbf{x},\mathbf{\sigma })=Z^{-1}\exp (-\beta V(\mathbf{x},\mathbf{\sigma }))d\mathbf{x},Z=\sum _{\mathbf{\sigma }}{}^{{}^{\prime }}{\int }_{{\Omega }^N}\exp (-\beta V(\mathbf{x},\mathbf{\sigma }))d\mathbf{x}.$$

$$\sum _{i=1}^N{\sigma }_i=M,N_{\mathrm{A}}+N_{\mathrm{B}}=N,N_{\mathrm{A}}-N_{\mathrm{B}}=M.$$

$$Z=\sum _{\mathbf{\sigma }}{}^{{}^{\prime }}{\int }_{{\Omega }^N}\exp (-\beta V(\mathbf{x},\mathbf{\sigma }))d\mathbf{x}=\left(\genfrac{}{}{0ex}{}{N}{N_{\mathrm{A}}}\right){\int }_{{\Omega }^N}\exp (-\beta V(\mathbf{x},{\mathbf{\sigma }}_{\star }))d\mathbf{x}.$$

$${\bar{\Omega }}^N={\cup }_{\ell =1}^n{\bar{D}}_{\ell }(\mathbf{\sigma }),D_{\ell }(\mathbf{\sigma })\cap D_{{\ell }^{\prime }}(\mathbf{\sigma })=\varnothing \text{if ℓ=ℓ′}.$$

$$Z=\sum _{\mathbf{\sigma }}{}^{{}^{\prime }}\sum _{\ell =1}^n{\int }_{D_{\ell }(\mathbf{\sigma })}\exp (-\beta V(\mathbf{x},\mathbf{\sigma }))d\mathbf{x}=\sum _{\mathbf{\sigma }}{}^{{}^{\prime }}\sum _{\ell =1}^n{Z}_{\ell }(\mathbf{\sigma }).$$

$$\begin{array}{rl}V(\mathbf{x},\mathbf{\sigma })& =\sum _{i<j}{\varphi }_{{\sigma }_i,{\sigma }_j}(|x_{ij}|)+\sum _i{u}_{{\sigma }_i}(|x_i|)\\ & =\sum _{i<j}\frac{1}{4}(1+{\sigma }_i)(1+{\sigma }_j){\varphi }_{++}(|x_{ij}|)\\ & +\frac{1}{4}[(1+{\sigma }_i)(1-{\sigma }_j)+(1-{\sigma }_i)(1+{\sigma }_j)]{\varphi }_{+-}(|x_{ij}|)\\ & +\frac{1}{4}(1-{\sigma }_i)(1-{\sigma }_j){\varphi }_{--}(|x_{ij}|)\\ & +\sum _i\frac{1}{2}(1+{\sigma }_i)u_{+}(|x_i|)+\frac{1}{2}(1-{\sigma }_i)u_{-}(|x_i|),\end{array}$$

$$\begin{array}{rl}V(\mathbf{x},\mathbf{\sigma })& =\sum _{i<j}-{\sigma }_i{\sigma }_j{J}_{ij}(\mathbf{x})+\sum _i{\sigma }_i{h}_i(\mathbf{x})+f(\mathbf{x}).\end{array}$$

$${\sigma }_k^{ij}=\{\begin{array}{ll}{\sigma }_i& k=j,\\ {\sigma }_j& k=i,\\ {\sigma }_k& k\ne i,j.\end{array}$$

$$\begin{split}\Delta_{ij}V(\mathbf{x},\bm{\sigma})&\equiv V(\mathbf{x},\bm{\sigma}^{ij})-V(\mathbf{x},\bm{\sigma})\\ &\quad=(\sigma_{j}-\sigma_{i})\Big{\{}-\sum_{k\neq i,j}(J_{ik}(\mathbf{x})-J_{jk}(\mathbf{x}))\sigma_{k}+h_{i}(\mathbf{x})-h_{j}(\mathbf{x})\Big{\}}.\end{split}$$

$$\mathbb{E}[{\sigma }_i]=s_i,\sum _{i=1}^N{s}_i=M.$$

$$\begin{array}{r}\nu (d\mathbf{x},\mathbf{\sigma })=Z_{\nu }^{-1}\exp (-\beta V(\mathbf{x},\mathbf{\sigma })+\beta \mathbf{\lambda }\cdot \mathbf{\sigma })d\mathbf{x},\\ Z_{\nu }=\sum _{\mathbf{\sigma }}{\int }_{{\Omega }^N}\exp (-\beta V(\mathbf{x},\mathbf{\sigma })+\beta \mathbf{\lambda }\cdot \mathbf{\sigma })d\mathbf{x}.\end{array}$$

$$\tilde{\nu }(d\mathbf{x},\mathbf{\sigma };\mathbf{X},\mathbf{k})={\tilde{Z}}_{\nu }^{-1}\exp (-\beta \tilde{V}(\mathbf{x};\mathbf{X},\mathbf{k})+\beta \tilde{\mathbf{\lambda }}\cdot \mathbf{\sigma })d\mathbf{x},$$

$$\tilde{V}(\mathbf{x};\mathbf{X},\mathbf{k})=\sum _{i=1}^N\frac{k_i}{2}|x_i-X_i{|}^2.$$

$${\tilde{Z}}_{\nu }$$

$${\tilde{Z}}_i$$

$${\tilde{\lambda }}_i$$

$$\mathcal{R}(\mu ||\nu )=\{\begin{array}{ll}{\mathbb{E}}^{\mu }\left[\log \frac{d\mu }{d\nu }\right]=\int \log \frac{d\mu }{d\nu }d\mu ,& \mu \ll \nu ,\\ \infty ,& \text{otherwise}.\end{array}$$

$$\mathcal{R}(\tilde{\nu }||\nu )=\beta {\mathbb{E}}^{\tilde{\nu }}[\Delta V]-\beta \Delta \mathbf{\lambda }\cdot \mathbf{s}+\log Z_{\nu }-\log {\tilde{Z}}_{\nu },\Delta V=V-\tilde{V},\Delta \mathbf{\lambda }=\mathbf{\lambda }-\tilde{\mathbf{\lambda }}.$$

$$(\mathbf{X},\mathbf{k})\in \mathrm{argmin}\mathcal{R}(\tilde{\nu }||\nu )$$

$$-{\beta }^{-1}\log Z_{\nu }+\mathbf{\lambda }\cdot \mathbf{s}\le {\mathbb{E}}^{\tilde{\nu }}[\Delta V]-{\beta }^{-1}\log {\tilde{Z}}_{\nu }+\tilde{\mathbf{\lambda }}\cdot \mathbf{s}=\mathcal{F}.$$

$$\tilde{\mu }(d\mathbf{x};\mathbf{X},\mathbf{k})=\prod _{i=1}^N{\tilde{Z}}_i^{-1}\exp \left\{-\frac{\beta k_i}{2}|x_i-X_i{|}^2\right\}d{x}_i=\prod _{i=1}^N{\tilde{\mu }}_i(d{x}_i;X_i,k_i),$$

$${\tilde{J}}_{ij}(\mathbf{X},\mathbf{k})={\mathbb{E}}^{\tilde{\mu }}[J_{ij}(\mathbf{x})],{\tilde{h}}_i(\mathbf{X},\mathbf{k})={\mathbb{E}}^{\tilde{\mu }}[h_i(\mathbf{x})],\tilde{f}(\mathbf{X},\mathbf{k})={\mathbb{E}}^{\tilde{\mu }}[f(\mathbf{x})],$$

$$\begin{array}{rl}\mathcal{F}& =\sum _{i<j}-s_i{s}_j{\tilde{J}}_{ij}(\mathbf{X},\mathbf{k})+\sum _i{s}_i{\tilde{h}}_i(\mathbf{X},\mathbf{k})+\tilde{f}(\mathbf{X},\mathbf{k})-\sum _i\frac{k_i}{2}{\mathbb{E}}^{{\tilde{\mu }}_i}[|x_i-X_i{|}^2]\\ & +{\beta }^{-1}\sum _i\frac{1+s_i}{2}\ln \frac{1+s_i}{2}+\frac{1-s_i}{2}\ln \frac{1-s_i}{2}-\ln {\tilde{Z}}_i.\end{array}$$

$$\frac{\partial \mathcal{F}}{\partial s_i}=-\sum _j{s}_j{\tilde{J}}_{ij}(\mathbf{X},\mathbf{k})+{\tilde{h}}_i(\mathbf{X},\mathbf{k})+{\beta }^{-1}\mathrm{arctanh}(s_i).$$

$$m_{ij}=\{\begin{array}{ll}{m}_{\mathrm{c}}& j\in N_i,\\ 0& j\notin N_i.\end{array}$$

$$m_{ij}=\{\begin{array}{ll}{m}_{\mathrm{r}}(1+s_i)(1-s_j),& j\in N_i,\frac{\partial \mathcal{F}}{\partial s_i}>\frac{\partial \mathcal{F}}{\partial s_j},\\ m_{\mathrm{r}}(1-s_i)(1+s_j),& j\in N_i,\frac{\partial \mathcal{F}}{\partial s_i}<\frac{\partial \mathcal{F}}{\partial s_j},\\ 0& j\notin N_i.\end{array}$$

$$\frac{d{s}_i}{dt}=\sum _{j\in N_i}(1+s_j)(1-s_i){\Gamma }_{j\to i}-(1+s_i)(1-s_j){\Gamma }_{i\to j}.$$

Full text

In this work, we formulate a spin-diffusion stochastic process and show how it can be connected to diffusive molecular dynamics. The spin-diffusion model couples a classical overdamped Langevin equation to a kinetic Monte Carlo model for exchange amongst the species of a binary alloy. Under suitable assumptions and approximations, spin-diffusion can be shown to lead to diffusive molecular dynamics type models. The key assumptions and approximations include a well defined time scale separation, a choice of spin-exchange rates, a low temperature approximation, and a mean field type approximation. We derive several models from different assumptions and show their relationship to diffusive molecular dynamics. Differences and similarities amongst the models are explored in a simple test problem.

Version:

1 Introduction

Metals, alloys, and other condensed matter typically exhibit behavior on at least two vastly different time scales making direct numerical simulation by classical molecular dynamics (MD) impractical. The vibrational time scale of the atoms, measured in femtoseconds ($10^{-15}$ s), sets a first time scale, and constrains the time step of MD. Physically relevant phenomena occur on time scales of microseconds ($10^{-6}$ s) or longer, setting a second time scale. Milliseconds of the laboratory time can be achieved with direct MD using special purpose hardware, but typical MD simulations only reach nanoseconds [39].

The source of this scale separation is the metastability present in the physical system. By metastability, we mean the tendency of the system to persist in a particular arrangement, or conformation, of atoms for a relatively long time before rapidly transitioning to some other persistent conformation. In many condensed matter problems, these distinct metastable states correspond to local minimizers of a potential energy landscape and transitions occur through saddle points. The system fluctuates about a local minimum on the fast time scale, before a rare, but rapid, transition occurs. As the feature of physical interest is the sequence of transitions, classical MD will exhaust computational resources resolving the fast, uninteresting, fluctuations. Overcoming this scale separation has led to significant investment in models and algorithms, including transition state theory (TST), kinetic Monte Carlo (kMC), accelerated molecular dynamics, and phase field crystal (PFC)  [30, 40, 42, 9, 41, 22, 19]. All of these methods exploit a “fast/slow” decomposition of the system, averaging or approximating the fast femtosecond time scale, while faithfully resolving details of the infrequent transitions between the metastable regions.

One method that has recently been formulated to overcome the time scale challenge is diffusive molecular dynamics (DMD)  [24, 35, 36]. While originally proposed based on physical grounds, here, we derive a DMD type system of equations through the following steps. We begin with the overdamped Langevin equation, a diffusion, as a model for a binary alloy. Some of the atoms in this system are species A and some are species B, and the type species asserts itself through the potential energy. Metastability also bedevils overdamped Langevin, limiting the usefulness of direct numerical integration. To circumvent the obstacle of metastability, we adjoin a spin-exchange kMC process, allowing for the two species, A and B, to swap. This allows for more rapid evolution of the system than could be observed by direct numerical simulation of overdamped Langevin. We refer to this coupled process as the spin-diffusion model.

Exploiting the time scale separation of the metastability using idealized quasistationary distributions (QSDs), the fast diffusional process averages out, leaving the kMC model with averaged reaction rates. This averaged kMC model then predicts the positions of the atoms, as they are distributed by a distribution conditioned on the current spin configuration. Finally, an approximation of the ensemble average of this scale separated kMC model gives rise to a system of ordinary differential equations. These predict the composition of the alloy and the structure of the distorted lattice. It thus provides a tool for the study of mechanical deformation and compositional evolution, as in the original DMD models. Our derivation, and the main goal of this work, bridges the gap between MD and DMD, and constrains terms in DMD by relating them to our underlying model.

The spin-diffusion model can be seen as a further generalization and approximation of the weakly off-lattice kMC model of [4]. In that work, the authors quenched their system to find the local minimizer, identified saddle points of the associated basin, and then used Arrhenius reaction rates for the associated kMC moves. As we do not insist on quenching, our model allows for finite temperature fluctuations about the local minima. Additionally, while the method in [4] is stochastic, requiring ensembles of simulations to study the evolution of averages, a DMD type approach evolves an approximation of the time dependent ensemble average. This predicts the evolution of observables from a single simulation.

DMD and spin-diffusion also bear resemblance to other models that have appeared in the literature, including time dependent density functional theory (TDDFT), phase field crystals (PFC), mean field kinetic equations (MFKE), and compressible Ising  [1, 13, 27, 26, 38, 12, 34, 29, 11]. All of these models seek to predict scalar quantities governing, for instance, probability of occupancy at a point in space by a particular species. Time dependent quantities evolve under systems of differential equations.

In contrast to TDDFT and PFC, DMD begins with the interatomic potentials which would be used in a traditional MD simulation and augments them to account for additional degrees of freedom. Thus, there is hope for better quantitative agreement. DMD is entirely classical in our formulation, and does not account for quantum effects. Also, unlike PFC, DMD explicitly retains atomistic detail and is not a continuum limit. Indeed, the DMD type models do not correspond to taking a large-volume, fixed-density limit, as in many mean field approximations [32, 14].

Much, but not all, of the work on MFKE, whose ideas we draw from, has been for on-lattice problems; DMD is explicitly off-lattice. The ideas which are closest to DMD include [34, 29] and related works, where the lattice is permitted to fluctuate under harmonic approximations. The distinction with DMD, which also makes harmonic approximations, is that its harmonic approximations are dynamically determined.

1.1 Diffusive Molecular Dynamics

The DMD method, as presented in the materials science literature, begins by reformulating the MD problem in an extended state space. Instead of having $N$ distinct atoms located at positions $\mathbf{x}=\{x_{i}\}_{i=1}^{N}$, one has $N$ distinct atomic sites located at $\mathbf{x}$ with associated scalar parameters $\mathbf{s}=\{s_{i}\}_{i=1}^{N}$. The parameter $s_{i}\in[-1,1]$ captures the likelihood of the site $i$ being occupied by the species A ($s_{i}=1$) or species B ($s_{i}=-1$). In the earliest DMD works, elemental materials were studied, and the scalar parameter was $c_{i}\in[0,1]$, representing probability of occupancy by an atom as opposed to a vacancy. We emphasize the binary alloy problem and work in the $s_{i}$ coordinate.

A free energy, denoted by $\mathcal{F}$, is constructed as a function of $\mathbf{s}$ and the mean atomic positions $\mathbf{X}=\{X_{i}\}_{i=1}^{N}$, along with harmonic parameters, $\mathbf{k}$, defined in Section 3.1. First, $\mathcal{F}$ is minimized over $\mathbf{X}$ and $\mathbf{k}$ at fixed $\mathbf{s}$:

[TABLE]

This corresponds to a variational Gaussian (VG) approximation [23]. Hence, on the fast time scale, the atomic sites, with a fixed composition, are vibrating about some well defined configuration with the mean positions $\mathbf{X}$. On the slow time scale, the composition $\mathbf{s}$ evolves. This minimization aspect of DMD has recently been studied and posed in terms of relative entropy (Kullback-Leibler divergence) [37]. The variational setting using the relative entropy was explored in [15] for constructing a coarse-grained dynamics in a simpler model of spin-diffusion on a fixed lattice describing non-Fickian diffusion under external concentration gradient. However, DMD represents a different approach in which the evolution of the coarse-grained model is deterministic and is defined by a system of ordinary differential equations.

The question is now how to evolve the composition in time. In [35], one proposed model is driven by the free energy gradients,

[TABLE]

Here, $m_{ij}$ is a mobility tensor and $\mathcal{N}_{i}$ is the set of atomic sites neighboring site $i$. Assuming the mobility is positive and symmetric (i.e., $m_{ij}=m_{ji}$ and $\sum_{ij}v_{i}m_{ij}v_{j}>0$ for all vectors $v\neq 0$), flow (1.2) conserves the total mass, $\sum s_{i}$, and has $\mathcal{F}$ as a Lyapunov function. Throughout these dynamics, $(\mathbf{X},\mathbf{k})$ always instantaneously satisfy (1.1), and thus the evolution is quasistatic.

DMD has been used to study vacancy diffusion problems like sintering of nanoparticles, diffusional void growth, and dislocation climb, along with defect migration and binary alloys [7, 8, 24, 36]. More generally, this method appears appropriate for modeling condensed matter systems with an essential coupling between their composition and their mechanical properties.

1.2 Spin-Diffusion Dynamics

System (1.2) is a modeling assumption, without a direct relationship to an underlying MD model. Furthermore, it requires additional modeling of the mobility factor. Here, we develop a stochastic model which bridges MD to DMD, and relates coefficients, like the mobility, to the underlying model.

For a binary alloy, each site has a position, $x_{i}$, and a “spin”, $\sigma_{i}\in\{\pm 1\}$, determining if it is species A ($\sigma_{i}=+1$) or B ($\sigma_{i}=-1$). The ensemble averages of the $\sigma_{i}$ will correspond to the $s_{i}$ of DMD. Fixing the composition at $\bm{\sigma}$, the positions evolve under the diffusion

[TABLE]

for a potential of both $\mathbf{x}$ and $\bm{\sigma}$. With a static value of $\bm{\sigma}$, (1.3) corresponds to a standard MD method.

Fixing the configuration $\mathbf{x}$, we now allow the composition to evolve under a spin-exchange model with Kawasaki dynamics [17, 16, 18]. Recall that such a process is defined by reaction rates, $r(\bm{\sigma}\to\bm{\sigma}^{\prime};\mathbf{x})$, assumed to satisfy detailed balance with respect to the Gibbs distribution, i.e.,

[TABLE]

Both of the above processes have $e^{-\beta V(\mathbf{x},\bm{\sigma})}$ as the invariant measure. If we let $L_{x}$ and $L_{\sigma}$ denote the generators of the diffusion and spin-exchange processes, respectively, then $L_{\epsilon}=\epsilon^{-1}L_{x}+L_{\sigma}$, induces a process that also has $e^{-\beta V(\mathbf{x},\bm{\sigma})}$ as the invariant measure. The parameter $\epsilon$ captures the scale separation between the two processes. In the $\epsilon\to 0$ limit, the positions, $\mathbf{x}$, are instantaneously distributed according to a distribution parametrized by the composition, $\bm{\sigma}$, while the composition evolves under the spin dynamics with rates averaged over $\mathbf{x}$. After then approximating the ensemble average, we obtain a system of equations with many of the same features as DMD. The actual coupling of the processes is more subtle, but this is the overarching idea.

1.3 Outline

This work is organized as follows. In Section 2, we present the state space that underlies these models and the associated potentials. The essential ideas of DMD are then reviewed in Section 3. After that, we present our stochastic spin-diffusion model in Section 4. Numerical examples are given in Section 5 and we conclude with a discussion in Section 6. Additional calculations and details appear in the Appendix.

2 State Space and Potential Energies

The framework of both DMD and the spin-diffusion model begins with the state space $(\mathbf{x},\bm{\sigma})\in\Omega^{N}\times\{\pm 1\}^{N}$. The domain $\Omega\subset\mathbb{R}^{d}$ may be all of $\mathbb{R}^{d}$ or a subset with boundary conditions. The pair $(x_{i},\sigma_{i})$ identifies the position and species of the $i$-th site. Pair potential interactions can then be defined as

[TABLE]

Pair interaction between two atomic sites thus depend on both the species at the sites and the distance between them. The reader can imagine more general potentials of this form, such as EAM. A key assumption that we make on the potential is

Assumption 1** (Permutation Invariance of Potential).**

Given $(\mathbf{x},\bm{\sigma})$, if we permute the indices of both $\mathbf{x}$ and $\bm{\sigma}$ in the same way to obtain $(\mathbf{x}^{\prime},\bm{\sigma}^{\prime})$, then $V(\mathbf{x},\bm{\sigma})=V(\mathbf{x}^{\prime},\bm{\sigma}^{\prime})$.

Given such a potential, we define the distribution

[TABLE]

The summation is primed to indicate that the net composition is constrained:

[TABLE]

$N_{\rm A}$ and $N_{\rm B}$ are fixed beforehand. The partition function $Z$ suffers from both high dimensionality and a combinatorial explosion. Notice that because of the invariance in Assumption 1, for any fixed $\bm{\sigma}_{\star}$ satisfying (2.3),

[TABLE]

Thus, any process which samples $e^{-\beta V(\mathbf{x},\bm{\sigma})}$ over $\mathbb{R}^{d\cdot N}\times\{\pm 1\}^{N}$ can be used to compute ensemble averages over $e^{-\beta V(\mathbf{x},\bm{\sigma}_{\star})}$ of observables which are also invariant to permutations, such as the energy. We also assume:

Assumption 2** (Finitely Many Minima of Potential).**

For a given $\bm{\sigma}$, $V(\mathbf{x},\bm{\sigma})$ has a finite number of minima in $\Omega^{N}$.

By Assumption 1, the number of minima is then invariant to $\bm{\sigma}$. The minima have corresponding basins of attraction, $D_{\ell}(\bm{\sigma})$. Excluding the measure zero set corresponding to points which converge to saddles, which we will not give further consideration to, these open sets disjointly partition the state space, up to the measure zero set. The union of their closures, $\bar{D}_{\ell}(\bm{\sigma})$, cover $\bar{\Omega}^{N}$; for each $\bm{\sigma}$ satisfying (2.3),

[TABLE]

The partition function, from (2.4) can be decomposed as

[TABLE]

The $Z_{\ell}(\bm{\sigma})$ are the partition functions of $1_{D_{\ell}(\bm{\sigma})}(\mathbf{x})e^{-\beta V(\mathbf{x},\bm{\sigma})}$, the restricted Gibbs distributions.

The potential we examine in this work includes both pair potential interactions and self-interaction trapping terms:

[TABLE]

where $x_{ij}\equiv x_{i}-x_{j}$. Expression (2.7) can be recast as

[TABLE]

See Appendix A for details of $J_{ij}$, $h_{i}$ and $f$.

As noted in the introduction, the stochastic process we develop involves exchanges of spins between atomic sites. Given composition $\bm{\sigma}$, a spin-exchange will be denoted $\bm{\sigma}^{ij}$, with components

[TABLE]

Viable exchanges will be limited to a set, $\mathcal{N}_{i}$, of neighboring sites, which may change in time. With fixed configuration $\mathbf{x}$, the change in the potential energy, which drives the exchanges, is

[TABLE]

3 DMD Model

Here we briefly summarize the ingredients of DMD. We refer the reader to [37] for a more detailed treatment.

3.1 Free Energy Approximation

First, an approximate free energy is developed. Towards this, constraint (2.3) is relaxed via a change of ensemble:

Assumption 3** (Change of Ensemble).**

Mass constraint (2.3) only holds on average; the mean spin values satisfy the total mass constraint,

[TABLE]

The mean values, $s_{i}$, are enforced by the Lagrange multipliers, $\lambda_{i}$, for the distribution

[TABLE]

The summation in (3.2) is no longer primed as the $\bm{\sigma}$ varies over all of $\{\pm 1\}^{N}$. To proceed further the aforementioned VG approximation is made [23]. First, we define $\tilde{\nu}$ as

[TABLE]

with harmonic potential

[TABLE]

Note, a variable with a $\tilde{}$ on it indicates it is associated with a VG approximation. For potential (3.4), many of the terms in (3.3) become explicit:

[TABLE]

$\mathbf{X}$ and $\mathbf{k}$ are then chosen to minimize the relative entropy between $\tilde{\nu}$ and $\nu$ [37]. Recall that the relative entropy metric between two probability measures, $\mu$ and $\nu$, on a common state space is given by

[TABLE]

Here, $\mu\ll\nu$ indicates that $\mu$ is absolutely continuous with respect to $\nu$. $\mathcal{R}$ is non-negative and vanishes if and only if $\mu=\nu$. For our problem,

[TABLE]

We thus make the approximation:

Approximation 1** **(Variational Gaussian/Relative Entropy

Approximation).

Replace $\nu$ by $\tilde{\nu}$ in the model by solving

[TABLE]

and then employ the free energy upper bound $\mathcal{F}$ in the dynamics:

[TABLE]

We simplify the energy expression, $\mathcal{F}$, for later comparison. First, let

[TABLE]

corresponding to the configurational portion of the distribution. Next, define

[TABLE]

where $J$, $h$, and $f$ are as in (2.8). With these definitions,

[TABLE]

Thus,

[TABLE]

3.2 Evolution Equations

To evolve $\mathbf{s}$ in time under (1.2), modeling assumptions are needed for the mobility. Simple mobility models include

[TABLE]

The constant mobility model is simplest, but the rate limited mobility, from [35], accounts for limitations to mass transfer due to the amount of A and B at each site.

An alternative to (1.2) is the so-called “Master Equation,” which appeared in [24, 35, 36, 7, 8],

[TABLE]

The $\Gamma_{j\to i}$ are the rates of mass transfer from the site $j$ to site $i$. They are explicitly modulated by steric factors, ensuring that there is both mass available to transfer from $j$ and space available at $i$. In [35], (3.18) was preferred to (1.2) due to reduced numerical stiffness, and it has been used in [24, 35, 36, 7, 8].

A common choice for the jump rate is

[TABLE]

With this choice, (3.18) will also conserve mass and have $\mathcal{F}$ as a Lyapunov function. In the above expression, $\kappa$ is a rate constant, $f_{i}$ is the “formation energy,” and $Q_{\rm m}$ is the activation energy for mass exchange.111 The coefficient of $(f_{i}-f_{j})$ is $\beta$ rather than $\frac{\beta}{2}$ because we work in the $s_{i}$ coordinate system. In principle, $Q_{\rm m}$ should be obtained by finding the saddle point associated with each pair exchange, but many practitioners take it to be a constant for efficiency.

4 Spin-Diffusion Model

In this section, we develop our spin-diffusion model and show how it captures the essential features of DMD. The physical intuition behind our model is that for the basins partitioning up the state space as in (2.5), the diffusion process will capture intrabasin motion while the spin-exchange process will capture interbasin transitions. This idealizes the metastability of the basins with respect to the diffusion, assuming we would never observe an exit by diffusion. Instead, the spin exchange migrates the system between basins. This requires the introduction of the idealized quasistationary distributions, QSDs, which characterize the distribution of diffusions conditioned to never leave a basin [20, 31, 25]. We review the key ideas of the QSD and relate them to the spin-diffusion in the next subsection. Subject to additional assumptions on the time scales, the reaction rates, the temperature, and a mean field approximation, we then obtain our evolution equations.

The reader may question the introduction of the QSD – why is it not sufficient to use the “naive” coupling of the free, unconditioned, diffusion given by (1.3) to the spin-exchange, (1.4)? If we were to proceed with the naive coupling, as sketched in Section 1.2, our joint process would be $L_{\epsilon}=\epsilon^{-1}L_{x}+L_{\sigma}$. This leaves the value of the scale separation, $\epsilon$, ill-defined. Heuristically, it is intended to relate the time for the diffusion to relax to a “local” equilibrium of the system, a basin of attraction $D_{\ell}(\bm{\sigma})$ and the time for the diffusion to escape this basin. But the only time scale associated with the free diffusion is for convergence to the Gibbs distribution on all of $\Omega^{N}$, and we wish to model the time scales of local relaxation and transition. While one might attempt to use restricted Gibbs, it is the QSD that gives both the desired characterization of “locality” and the time scales that define $\epsilon$.

Additionally, the naive joint process allows events that we intend to restrict. Suppose $(\mathbf{x}(0),\bm{\sigma}(0))$ corresponds to Figure 1 (a). There is a nonzero chance the diffusion exits the basin before any spin-exchange occurs, transitioning, at time $t$, to Figure 1 (b), where $(\mathbf{x}(t),\bm{\sigma}(t))=(\mathbf{x}(0)^{7,10},\bm{\sigma}(0)^{7,10})$. Here, the two atomic sites, and the their species, have exchanged through the diffusion. Alternatively, a spin-exchange may occur before the diffusion can exit the basin, and the system may now be as in Figure 1 (c), with $(\mathbf{x}(t),\bm{\sigma}(t))=(\mathbf{x}(0),\bm{\sigma}(0)^{7,10})$. Here, only the spins, but not the sites, have been exchanged. Figures 1 (b) and (c) have the same physical properties, but have been obtained by different mechanisms from (a). Using our modified joint process, where the diffusion is conditioned to never leave the basin, conforms to our design to only allow transitions via spin-exchanges, disambiguating the evolution of the system. This conditioned diffusion converges to a QSD.

4.1 Quasistationary Distributions

Here, we review the essential features of the QSDs for diffusions and apply them to the spin-diffusion model.

4.1.1 Properties

First consider a diffusion independent of the spins,

[TABLE]

Excluding a set of measure zero converging to saddles, at any instant, $\mathbf{x}(t)$ is in a basin of attraction $D$ of the energy $V(\mathbf{x})$, which can be identified by the local minimum to which the system quenches:

[TABLE]

and the basin is the set that quench to the same minimum:

[TABLE]

Given $\mathbf{x}(0)=\mathbf{x}_{0}$, we equip the generator of this process $L_{x}$ with Dirichlet boundary conditions on $D(\mathbf{x}_{0})$. Consider the eigenvalue problem

[TABLE]

The eigenfunctions and eigenvalues depend parametrically on $\mathbf{x}_{0}$ through $D=D(\mathbf{x}_{0})$.

We assume that $D(\mathbf{x}_{0})$ is metastable, which is to say that the relaxation time scale to the QSD is far shorter than the first exit time. This can be quantified in terms of the first two eigenvalues:

[TABLE]

where we use $\approx$ to indicate these are characteristic (i.e. order of magnitude) time scales. Letting $T$ be the first exit time, the QSD $\check{\mu}$ relates to the diffusion via the limit

[TABLE]

and it also has a density with respect to the Lebesgue measure

[TABLE]

The QSD is not an invariant measure for the diffusion, since $\mathbf{x}(t)$ will exit $D(\mathbf{x}_{0})$ in finite time a.s. It is also not the Gibbs distribution restricted to $D(\mathbf{x}_{0})$; by construction, the QSD vanishes on $\partial D(\mathbf{x}_{0})$ and $e^{-\beta V(\mathbf{x})}$ will not. However, in the low temperature limit, these two distributions coincide. The eigenfunction $\varphi_{1}$ in (4.7) converges to a positive constant in the interior of the set $D(\mathbf{x}_{0})$, with a boundary layer ensuring it vanishes on $\partial D(\mathbf{x}_{0})$. This has been analyzed in [6, 21].

4.1.2 QSD Processes

As noted, a diffusion, even one whose initial position is sampled from the QSD, will leave the metastable region in finite time, a.s. Thus, its invariant measure will not be the QSD, and another stochastic process which has this as its invariant measure is needed. We thus assume

Assumption 4** (Ergodicity of a QSD Process).**

Given $\beta$ and a metastable region $D$, there exists a process $\check{\mathbf{x}}(t)$ ergodic with respect to the QSD $\check{\mu}$. Its generator, $\check{L}_{x}$ has only the constants in its kernel, while $\check{L}_{x}^{\dagger}$ has only the density of the QSD, with respect to Lebesgue, in its kernel.

In what follows, the choice of QSD process is not essential. All that is required is that a process with the QSD as its invariant measure and generator $\check{L}_{x}$ exists. An example of such a process is a branching-interacting particle system of $M$ diffusions, each independently obeying (4.1) until one exits $D(\mathbf{x}_{0})$. At that instant, the one which exits is killed and it is resampled from the $M-1$ survivors; see [3] and references therein. The process, denoted $\hat{\mathbf{x}}(t)\in D^{M}$ reflects the entire ensemble,

[TABLE]

with each $\check{\mathbf{x}}^{(j)}(t)\sim\check{\mu}$, in equilibrium. Indeed, $\hat{\mathbf{x}}$ has generator $\hat{L}_{x}$ with invariant measure $\hat{\mu}$. We relate this back to the QSD of interest using a projection, $\hat{P}$, $\check{\mathbf{x}}(t)=\hat{P}(\hat{\mathbf{x}}(t))=\check{\mathbf{x}}^{(1)}(t)$, so that

[TABLE]

which is assumed to hold for all bounded, continuous observable functions $\mathcal{O}$.

Quantities associated with the QSD will be indicated by $\check{}$ marks. This is in contrast to the $\tilde{}{\;}$’s appearing on VG averaged quantities like (3.13).

4.1.3 QSDs for the Spin-Diffusion Model

While the basins of attraction are defined in terms of the $\mathbf{x}$ argument of $V(\mathbf{x},\bm{\sigma})$, they depend parametrically on $\bm{\sigma}$, as in (2.5):

[TABLE]

$\varphi_{1}(\mathbf{x};\bm{\sigma},\ell)$ is the principal eigenfunction of the generator on the set $D_{\ell}(\bm{\sigma})$.

A consequence of the introduction of the QSD framework is that not only do the spins change during the jump process, $\bm{\sigma}\to\bm{\sigma}^{\prime}$, but also the basin of attraction, $D_{\ell}(\bm{\sigma})\to D_{\ell^{\prime}}(\bm{\sigma}^{\prime})$. This is due to the energy landscape $V(\mathbf{x},\bm{\sigma})$ changing to $V(\mathbf{x},\bm{\sigma}^{\prime})$; see Figure 2.

The jump process is now on the joint space $(\bm{\sigma},\ell)$. This is in contrast to the process sketched in the introduction, where the reaction rates were assumed to satisfy (1.4). Instead, the reaction rates satisfy detailed balance with respect to the QSDs: for all $\mathbf{x}\in D_{\ell}(\bm{\sigma})\cap D_{\ell^{\prime}}(\bm{\sigma}^{\prime})$, the rates are

[TABLE]

and defined to be zero outside of this intersection. Thus, if $D_{\ell}(\bm{\sigma})\cap D_{\ell^{\prime}}(\bm{\sigma}^{\prime})=\emptyset$, no exchange is possible. At the moment, we continue to leave the actual reaction rates unspecified, and merely assume rates satisfying (4.11) exist. At fixed $\mathbf{x}$, the generator of the exchange process is

[TABLE]

4.2 Scale Separation

The joint process for our QSD spin-diffusion has generator

[TABLE]

We now characterize $\epsilon$, by assuming:

Assumption 5** (QSD Time Scale Separation).**

There is a typical value, $0<\epsilon\ll 1$, such that in each metastable set, $D_{\ell}(\bm{\sigma})$,

[TABLE]

with $\lambda_{i}(\ell,\bm{\sigma})$ the eigenvalues associated with the QSD on $D_{\ell}(\bm{\sigma})$.

Again, $\approx$ here indicates that all of these ratios are on the same order magnitude. With this choice of $\epsilon$, typical spin exchanges occur on a time scale that is at least as large as a typical diffusion first exit time from any particular basin. Thus, our joint process respects the physical time scales of the diffusion.

Given an observable $v(\mathbf{x},\bm{\sigma},\ell)$, its evolution under the process associated with $\check{L}_{\epsilon}$ is

[TABLE]

Substituting the ansatz $v=v_{0}+\epsilon v_{1}+\epsilon^{2}v_{2}+\ldots$ into this equation leads to

[TABLE]

At $\operatorname{O}(\epsilon^{-1})$, we obtain

[TABLE]

Since the kernel of $\check{L}_{x}$ is assumed to be observables independent of $\mathbf{x}$, $v_{0}=v_{0}(\bm{\sigma},\ell,t)$. At the next order

[TABLE]

For the expansion to be consistent, we must satisfy the solvability condition,

[TABLE]

The leading order averaged equation is then

[TABLE]

We thus define the QSD averaged reaction rates,

[TABLE]

These satisfy detailed balance with respect to

[TABLE]

Therefore,

Approximation 2** (Leading Order Process with QSDs).**

Spins $\bm{\sigma}(t)$ and basins $\ell(t)$ evolve as a jump process with rates (4.20). Atomic sites $\mathbf{x}(t)$ are instantaneously distributed by $\check{\mu}_{\ell(t)}(\bullet\mid\bm{\sigma}(t))$.

4.3 Master Equation

Equation (4.19) and the averaged reaction rates (4.20) enable us to write down a master equation for the evolution of the mean spin. Let $p_{t}(\bm{\sigma},\ell)$ be the time dependent distribution of $(\bm{\sigma},\ell)$. This satisfies

[TABLE]

Since we only allow the spin-exchanges to $\bm{\sigma}^{\prime}=\bm{\sigma}^{jk}$, with neighboring sites, $j$ and $k$, (4.22) becomes

[TABLE]

Applying this to the average of spin at site $i$:

[TABLE]

Exchanging the roles of $(\bm{\sigma},\ell)$ and $(\bm{\sigma}^{jk},\ell^{\prime})$ in the first summation of (4.24),

[TABLE]

Following standard treatments, such as [17, 16, 28, 12], this allows us to write (4.24) as

[TABLE]

To obtain a tractable equation, further approximations are necessary.

4.4 Low Temperature Approximation

Next we make the approximation that as $\beta\to\infty$, $\varphi_{1}$ tends to a constant on $D_{\ell}(\bm{\sigma})$. Thus, in the low temperature regime, the QSD is approximately the restricted Gibbs distribution [6, 21].

Approximation 3** (Low Temperature QSD Approximation).**

$\beta$* is sufficiently large that*

[TABLE]

Hence, we can approximate (4.11) by

[TABLE]

for $\mathbf{x}\in D_{\ell}(\bm{\sigma})\cap D_{\ell^{\prime}}(\bm{\sigma}^{\prime})$ and zero otherwise. Since $\ell$ and $\ell^{\prime}$ only appear in the definition of regions $D_{\ell}(\bm{\sigma})$ and $D_{\ell^{\prime}}(\bm{\sigma}^{\prime})$, it is easy to formulate reaction rates satisfying this approximate detailed balance.

4.5 State Space Approximation

A major challenge in applying this QSD-based formulation is that we must identify all the states neighboring $(\bm{\sigma},\ell)$ to compute the associated reaction rates. While the $\bm{\sigma}^{\prime}$ are easy to identify, finding the sets $D_{\ell^{\prime}}(\bm{\sigma}^{\prime})$ that intersect $D_{\ell}(\bm{\sigma})$ is intractable in high dimension.

A simple approximation is to only consider a single $\ell^{\prime}$ for each $\bm{\sigma}^{\prime}$, given $(\bm{\sigma},\ell)$. This corresponds to finding the set $D_{\ell^{\prime}}(\bm{\sigma}^{\prime})$ containing the local minimum of $D_{\ell}(\bm{\sigma})$. Consequently, if $\bm{\sigma}\to\bm{\sigma}^{\prime}$, the system will undergo a displacive relaxation from the previous local minima to the new one. In the context of Figure 2, if we begin in $D_{2}(\bm{\sigma})$, then, since its local minima is in $D_{2}(\bm{\sigma}^{\prime})$, we go to that metastable region when $\bm{\sigma}\to\bm{\sigma}^{\prime}$. Under this approximation given $\ell(0)$, $\ell(t)=\ell(\bm{\sigma}(t))$:

Approximation 4** (State Space Approximation).**

Given the time series $\bm{\sigma}(t)$ and the initial basin of attraction indexed by $\ell(0)$, we can uniquely infer $\ell(t)$ and the associated basin of attraction for all time.

Thus (4.26) is replaced with

[TABLE]

with the notation $\ell(\bm{\sigma})$ and $\ell^{\prime}(\bm{\sigma},\bm{\sigma}^{ij})$ indicating that these are uniquely determined by the given spin configuration, along with the initial conditions.

4.6 Mean Field Approximations

The next step is to make mean field approximations, given particular reaction rates, to obtain a system of equations entirely in terms of $s_{i}=\mathbb{E}[\sigma_{i}]$. We consider two possible reaction rates here. However, we emphasize that in this approximation we are not taking the large volume, fixed density, limit; discreteness of the system is retained.

4.6.1 $\tanh$ Reaction Rates

A simple choice of reaction rate satisfying (4.28) involves the $\tanh$ function:

[TABLE]

with $\Delta_{ij}V(\mathbf{x},\bm{\sigma})$ given in (2.10), and $\tau$ a rate constant. Note that through the introduction of $\epsilon$ in (4.13), after integrating out to time $t=\operatorname{O}(\epsilon\lambda_{1}^{-1})=\operatorname{O}((\lambda_{2}-\lambda_{1})^{-1})$ an exit event would have occurred in the underlying diffusion, and now a spin-exchange should occur in the joint process. Hence, $\tau=\operatorname{O}((\lambda_{2}-\lambda_{1})^{-1})$. After the scale separation, units of time and $\tau$ can be arbitrarily rescaled, jointly.

Following [28], the conditional average is:

[TABLE]

No approximation has yet been made in (4.31). To obtain equations of motion entirely in terms of $\mathbf{s}(t)$, we now make:

Approximation 5** **(Mean Field Approximation for $\tanh$ Reaction

Rates).

Assume that:

• •

The averages over $\check{\mu}_{\ell}(\bullet\mid\bm{\sigma})$ approximately commute with the $\tanh$ function.

• •

$\check{\mu}_{\ell}(\bullet\mid\bm{\sigma})$* can be approximated by $\check{\mu}_{\ell}(\bullet\mid\mathbf{s})$.*

• •

Spins are approximately uncorrelated, $\mathbb{E}[\sigma_{i}\sigma_{j}]\approx\mathbb{E}[\sigma_{i}]\mathbb{E}[\sigma_{j}]$.

Then (4.29) can be approximated by the system

[TABLE]

with atomic sites $\mathbf{x}(t)\sim\check{\mu}_{\ell(\mathbf{s})}(\bullet\mid\mathbf{s})$.

The QSDs $\check{\mu}_{\ell}(\mathbf{s})$ are defined analogously to $\check{\mu}_{\ell}(\bm{\sigma})$ with $s_{i}$ taking the place of $\sigma_{i}$ in the potential in (2.8). The basins are now defined as $D_{\ell}(\mathbf{s})$, which are implicitly assumed to still satisfy a corresponding version of Assumption 2. The argument of the $\tanh$ function in (4.32) is

[TABLE]

with

[TABLE]

These are the QSD averaged analogs of the VG averages in (3.13). The symbol ${\ell(\mathbf{s})}$ is used in the above expressions to remind the reader that we have made Approximation 4 on the state space. For brevity, we will suppress it in what follows.

4.6.2 Arrhenius Reaction Rates

Another choice is the Arrhenius type reaction rates, subject to a constant saddle height approximation:

[TABLE]

In analogous fashion to Approximation 5,

Approximation 6** **(Mean Field Approximation for Arrhenius Reaction

Rates).

Assume that:

• •

The averages over $\check{\mu}_{\ell}(\bullet\mid\bm{\sigma})$ approximately commute with the $\tanh$ function.

• •

$\check{\mu}_{\ell}(\bullet\mid\bm{\sigma})$* can be approximated by $\check{\mu}_{\ell}(\bullet\mid\mathbf{s})$.*

• •

Spins are approximately uncorrelated, $\mathbb{E}[\sigma_{i}\sigma_{j}]\approx\mathbb{E}[\sigma_{i}]\mathbb{E}[\sigma_{j}]$.

Then (4.29) can be approximated by the system

[TABLE]

with atomic sites $\mathbf{x}(t)\sim\check{\mu}_{\ell(\mathbf{s})}(\bullet\mid\mathbf{s})$.

4.7 Comparison with DMD

In this section we compare the models (4.32) and (4.36), with the DMD models (1.2) and (3.18). Our models are close in spirit to DMD dynamics, with the significant difference being the use of terms like $\tilde{J}_{ij}^{(\ell)}(\mathbf{s},\mathbf{X},\mathbf{k})$ from (3.13) in place of $\check{J}_{ij}^{(\ell)}(\mathbf{s})$. This reflects that, in standard DMD models, a VG approximation of a single mode of the Gibbs distribution is used, while we use a QSD. Both distributions capture local features of $e^{-\beta V(\mathbf{x},\bm{\sigma})}$. Furthermore, in the low temperature limit, the QSD will be well approximated by restricted Gibbs (see Approximation 3), and individual modes of $e^{-\beta V(\mathbf{x},\bm{\sigma})}$, will be well approximated by restricted Gaussians (i.e., with harmonic potentials). Thus, we expect agreement of terms like $\tilde{J}_{ij}^{(\ell)}(\mathbf{s},\mathbf{X},\mathbf{k})$ and $\check{J}_{ij}^{(\ell)}(\mathbf{s})$ in this limit.

4.7.1 $\tanh$ Reaction Rates

For the $\tanh$ reaction rate model, suppose we make a high temperature approximation of (4.32) to expand the $\tanh$ function

[TABLE]

Next, we use (3.15) and (2.10) to write the difference in the free energy gradients as

[TABLE]

Under a constant mobility model, we can write (1.2) as

[TABLE]

Next, neglecting boundary conditions, in the high temperature, near equilibrium limit, we expect, from rigid lattice asymptotics, $s_{i}=M/N+\operatorname{O}(\beta)=s_{\operatorname{eq}}+\operatorname{O}(\beta)$, so that (4.37) becomes

[TABLE]

while (4.39) becomes

[TABLE]

Comparing (4.39) with (4.40), we see what further conditions must hold for the two models to be approximately equivalent. First, the VG approximation must be suitable,

[TABLE]

Next, there is also the additional term, $\beta\tilde{J}_{ij}(s_{j}-s_{i})$, at $\operatorname{O}(\beta)$ which was discussed in [28], in the case of a rigid lattice. In that work, the author commented that such terms motivated the use of the mean field kinetic models, in the spirit of (4.32), rather than the free energy gradient model, (1.2). Finally, the mobility should scale as

[TABLE]

which is consistent with [35]. A similar analysis holds the rate limited mobility model, where, in the high temperature near equilibrium limit,

[TABLE]

resulting in

[TABLE]

A serious weakness of this analysis is the use of high temperature approximations. This specifically voids Approximation 3, and it degrades the quality of the VG approximation. Nevertheless, it provides some basic insight into how the models compare, and offers rough scalings for constants like (4.43) and (4.44).

4.7.2 Arrhenius Reaction Rates

For the Arrhenius reaction rate model, we substitute identity (4.38) into (3.19) and (3.20) to express the DMD model (3.18) as

[TABLE]

Comparing with (4.36), we see a distinction in the multiplicative terms $\exp(\pm\beta\tilde{J}_{ij}(s_{j}-s_{i}))=1+\operatorname{O}(\beta)$. Again, this corresponds to the term identified in [28]. We thus infer

[TABLE]

provided (4.42) holds. Models (3.18) and (4.36) have similar equations, and there is no need for an expansion in $\beta$, excepting the term $\exp(\pm\beta\tilde{J}_{ij}(s_{j}-s_{i}))$. In the high temperature limit, all of these models tend to agree, to leading order.

5 Numerical Simulations

Here, we compare some of the models presented in this paper with a simple test problem.

5.1 Test Problem

Our test problem is a 1D binary alloy, consisting of a chain of atomic sites that are occupied either by Species A or Species B. The pair potential $\phi_{\sigma_{i},\sigma_{j}}(|x_{ij}|)$ is such that is energetically favorable to have like species neighbor one another (A-A or B-B bonds). Additionally, A-A and B-B bonds have longer equilibrium interatomic spacings than A-B bonds resulting in mechanical deformation in response to rearrangement of the species. See Figure 3 for the problem setup. While this problem has certain nonphysical aspects, it captures the essential features of our spin-diffusion model. This allows for comparisons with a number of deterministic models, including DMD free energy gradient dynamics and mean field master equations.

A few additional details on this test problem are of note. Owing to the translation invariance of the problem, the leftmost atom is pinned at $x_{1}=0$, while the rightmost atom is free to move. The species of the atoms at the first and last atomistic sites are fixed, and only $s_{i}$, $i=2,\ldots,N-1$, evolve under the flow.

With respect to energy model (2.7), there is no $u_{\sigma}$ contribution, and the pair potential is given by the smoothed Lennard-Jones type potential (A.1) with cutoff (A.2). The variations to the pair potentials by bond type are controlled by the choice of parameters. We also add the confining pair potential (A.3) between adjacent atoms to ensure that there is no disaggregation of the chain.

5.2 Model Comparison

As discussed throughout this work, there are many choices of dynamics. Here, we compare the following:

• •

The spin-diffusion model, with scale separation, low temperature approximation, and the state space approximation, corresponding to (4.29) with rates (4.30).

• •

Free energy gradient dynamics (1.2), with constant mobility model (3.16).

• •

Free energy gradient dynamics (1.2), with mobilities defined by (3.17).

• •

Mean field Approximation 5, given by equation (4.32), with an additional VG approximation of $\Delta_{ij}\check{V}(\mathbf{s})$, as in (4.42).

• •

Mean field Approximation 5, given by equation (4.32), with an additional point estimate of $\Delta_{ij}\check{V}(\mathbf{s})$:

[TABLE]

with $(\mathbf{X}_{\min},\mathbf{k}_{\min})$ the local minimizer of $\mathcal{F}(\mathbf{s},\mathbf{X},\mathbf{k})$.

To fix a time scale, we choose $\tau=1$ in (4.32). The free energy gradient dynamics with mobilities (3.16) and (3.17) have free parameters $m_{\rm c}$ and $m_{\rm r}$, respectively. These are chosen so that the mean spins evolve on a similar time scales in both the mean field and free energy gradient models. While we calculated scalings (4.43) and (4.44) in the high-temperature, near-equilibrium limits, they might not extend to low temperatures.

Indeed, we choose $m_{\rm c}$ and $m_{\rm r}$ empirically to align a reference feature of the dynamics across the different models. Here, this feature corresponds to the first jump in the minimum of the spatial spin autocorrelation,

[TABLE]

The first minimum of $R(j)$ corresponds to the characteristic width of a spatial domain. The first jump in this minimizer, estimated as a continuous quantity via interpolation, is shown in Figure 4. This is further discussed below. After running the simulations with $m_{\rm c}=4$ and $m_{\rm r}=1$, time is rescaled to match the reference feature. While this fitting is done independently for each $N$ , the values in Table 1 reveal consistency across the values.

5.3 Numerical Implementation

For initial conditions, the $\sigma_{i}$ of the stochastic model alternate $\pm 1$, while in the deterministic models, the $s_{i}$ alternate $\pm 0.9998$. Given this initial spin composition, the initial conditions in the stochastic model are obtained by quenching the energy to a minimizer, while $\mathbf{X}$ and $\mathbf{k}$ are chosen to locally minimize the free energy via gradient descent.

Our simulation of (4.29) is performed by integrating the overdamped Langevin until it samples the QSD, after which a spin-exchange is performed using the time averaged reaction rates. The diffusion is integrated using preconditioned MALA.

In the deterministic models, the $\mathbf{s}$ are evolved by one of the choices of dynamics, and, since each model involves a VG approximation, the $\mathbf{X}$ and $\mathbf{k}$ are evolved by a scaled gradient descent for the free energy. Thus, the system that is actually solved is:

[TABLE]

Here, $f$ comes from one of the models mentioned above. The constants $\eta_{1}$ and $\eta_{2}$ are large numerical parameters that force the system to be close to the local minimizer of $\mathcal{F}(\mathbf{s},\mathbf{X},\mathbf{k})$.

This system of ODEs is solved in Matlab using the ode15s stiff solver. If the magnitude of $\nabla_{\mathbf{X},\mathbf{k}}\mathcal{F}$ exceeds a threshold tolerance, the integrator is halted and gradient descent is applied to obtain a local minimizer, again. After minimization, integration of (5.3) resumes. For efficiency, quadrature is performed using the GSL implementation of QUADPACK routines, http://www.gnu.org/software/gsl/, while gradient descent minimization is done by running Matlab’s ode15s on the gradient dynamics to steady state.

Additional details on our numerical schemes are given in Appendices B and C.

5.4 Results

Here, we present results of the simulation for a small system with $N=32$, allowing us to observe near complete coarsening, along with larger systems with $N=500$, $1000$ and $2000$, to see more physically meaningful behavior.

5.4.1 Small System

In the case of the small system, $N=32$, the composition and configuration profiles are shown in Figure 5. The composition dynamics has three main features:

Mixing: Initially alternating positive and negative spins become almost constant zero. 2. 2.

Nucleation: Phases of positive and negative spins form, initially at the boundaries, and then in the interior. 3. 3.

Coarsening: Phases merge, resulting in larger phases.

Each of these changes in the composition occurs in conjunction with lengthening of the configuration. These features are shown in Figure 5. In, for instance, Figure 5 (c), we see mixing occurs between $t=0$ and $t\approx 1$, nucleation follows between $t\approx 1$ and $t\approx 10$, and then a coarsening event occurs around $t=5\times 10^{3}$.

In this small system, there are a few coarsening events as the system tends towards total segregation of the two species. Near total segregation is visible in the ensemble average, Figure 5 (b), while the deterministic simulations exhibit incomplete coarsening. It is clear that nucleation completes earlier for the mean field master equation model than for free energy gradient dynamics. For $N=32$, $\beta=160$, nucleation completes around $t=10$ for the mean field models, but closer to $t=10^{2}$ for the free energy gradient dynamics. Another distinction is that the mean field models have better defined grains, with sharper interfaces at this temperature. The stages of mixing and nucleation are more difficult to distinguish in the ensemble average, but we do see a coarsening event there around $t=10^{4}$.

Differences between the models are also visible when we examine the strain and fraction of A-B bonds as functions of time; these appear in Figure 6. Observables for the stochastic model evolve without plateauing. We believe the plateaus are a consequence of the mean field approximation and and the use of a system with comparatively short range interactions in the model at hand. Despite this discrepancy, a large fraction of each of the deterministic trajectories is within one standard deviation of the ensemble averages.

After mixing occurs, the system has a more rapidly increasing strain and decreasing number of A-B bonds in the mean field approximations than in the free energy gradient dynamics. The strain fields all agree in the long time limit, and, while there is still some disagreement in the fraction of A-B bonds, as the system tends to full segregation, these too will agree amongst the models.

5.4.2 Large Systems

We also investigated the DMD model in larger systems. The time horizon of the simulated dynamics was limited to $t=5\times 10^{4}$ due to the computational complexity of the stochastic model. While such time horizon is too short to see total segregation, we are still able to observe nucleation, mixing, and multiple coarsening events.

Looking at the autocorrelation in Figure 4 (b), we see coarsening of the system in all cases. Turning to Figures 7 (c) and (d), we see that, as in the smaller system, nucleation completes earlier for the mean field model, around $t=10^{2}$, than for free energy gradient dynamics, around $t=10^{3}$. A large number of coarsening events occur between $t=10^{3}$ and $t=10^{4}$ for the free energy gradient model, about an order of magnitude earlier than in the mean field model.

There are much less well defined phases in the ensemble average, as seen in Figure 7 (b), though there is still coarsening. In particular, the interfaces between the regions of segregated species are less well defined resulting in a “slush.” This can partly be explained by examining Figure 7 (a), where we see there is uncertainty as to where the exact locations of the interfaces are, resulting in an averaging out.

Turning to observables in Figure 8, we see disagreement between the ensemble average and the deterministic models. The smaller $N=32$ system, in comparison, showed better agreement between modes; see Figure 6.

Not only are the plateaus absent, but the deterministic dynamics are far from the distribution of the stochastic ensemble. We posit that this is a consequence of the mean field approximation.

There is consistency across the deterministic models of different sizes $N=500,1000,2000$ in the long time limit of the observables. For strain, there is even agreement with the ensemble average of stochastic trajectories. The long time behavior of the stochastic model is less consistent in the A-B bond fraction observable, but we conjecture we would observe improved agreement as the segregation continues and eventually completes.

6 Discussion

We have formulated and related an underlying stochastic process, the coupled spin-diffusion model, to DMD, and shown that there is qualitative agreement in simulations. Our simulations also show that, depending on the modeling parameters, there are differences in the dynamics. Our contention is that, given that such distinctions occur, practitioners should follow the strategy developed here, and put their key modeling assumptions in the reaction rates and the mean field approximations.

There is clear disagreement between the ensemble average of the stochastic model and the deterministic model, in that the ensemble averages have less sharply defined interfaces between phases, as in Figures 7, and the observables do not plateau as in the deterministic model, as in Figure 8. We believe this is the result of studying a system in one dimension with comparatively short range interactions. If we look at a system with an ABABAB… composition, the values of the $J_{ij}$ in (2.8) are as in Table 2. These are obtained from the equilibrium displacements obtained after quenching. The values indicate that while the system certainly has the first and second nearest neighbor interactions, the third nearest neighbors are comparatively weak, and beyond the third, there is negligible interaction.

In comparison, the mean field approximation of the on lattice Ising model is not expected to be accurate in one dimension, with such short range interactions, [5, 2, 32, 10]. Indeed, Figure 9 plots the time evolution for the on-lattice analog of the test problem presented in Section 5. Here, the energy is,

[TABLE]

We see similar discrepancies to what was observed in Figure 8. We expect such discrepancies will be mitigated in higher dimensional problems where the number of interacting sites increases and longer range potentials are applied.

This serves as a warning for the practitioner who seeks to use DMD type models. If the system in question only has short range interactions, such discrepancies will occur. It is an important problem to formulate and explore 2D and 3D test problems for these models, and to examine the limitations of the mean field approximation. We do expect agreement in the observables in the long time limit. Thus, even if the detailed temporal evolution, e.g., the transient events, is incorrect, the deterministic model is more easily integrated to equilibrium than the corresponding stochastic model.

We summarize the key assumptions and approximations that are at the core of the studied model.

Spin-Diffusion Model

The underlying modeling assumption is that the system is governed by the coupled system of (1.3) and (1.4).

Time Scale Separation

The next assumption is that there is adequate scale separation between the relaxation time and the time scale for barrier crossings. This is captured by Assumption 5, which is the basis for Approximation 2, the rate averaged jump process.

Low Temperature Approximation

To make the computation of the effective reaction rates tractable, the QSD is approximated by the restricted Gibbs distribution (Approximation 3), which in turn may be approximated by a Gaussian distribution using VG. This is sensible in the low temperature limit, and was captured by (4.42) in our comparison of the models.

State Space Approximation

The first major approximation is in Approximation 4, which restricts our attention to a single new basin each time a spin-exchange occurs.

Mean Field Approximation

The final approximation, giving rise to equations like (4.32), is to approximate expectations of nonlinear functions by the nonlinear functions of the expectations.

The time scale separation and low temperature approximations are quite sensible, are consistent with data, and they can be made rigorous. At present, the state space and mean field approximations are somewhat less constrained. The mean field approximation is quite common, and can even be justified in certain limits. The state space approximation merits further investigation.

One interpretation of our model is as a mean field approximation to a master equation for an off-lattice, finite temperature, restricted kMC method. It is clearly off-lattice as the atomic sites evolve in time, and it certainly has finite temperature effects throughout. We view it as a kMC model because the system is driven by the evolution of the composition, i.e., the assignment of the A and B species to different sites.

Indeed, one might arrive at our model by beginning with, for instance, the off-lattice kMC model of [4]. Instead of viewing the positions of the atoms as corresponding to their quenched configurations, one could allow them to sample the local mode of the energy landscape, as is done in accelerated MD approaches [30, 40]. There, the system reaches a local equilibrium, sampling a QSD, and then one of several approaches is used to accelerate the observation of a first exit. As noted, another difference between DMD type methods and the method of [4] is that the latter produces single realizations of the trajectories, requiring ensemble averaging.

This difference is also found when comparing with accelerated MD methods. Accelerated MD can efficiently produce single trajectories, but will require ensemble averaging to predict the evolution of observables. Another difference is that accelerated MD methods allow for all types of transitions amongst basins, while DMD privileges species exchange, which correspond to correlated events. This restricts the underlying kMC model of DMD. Other off-lattice kMC methods, such as those in [42, 43], are constructed from transitions between adjacent basins of $V(\mathbf{x})$ connected by a single saddle point. The exchanges of species would not be captured by single kMC event in such a model, requiring multiple transitions.

There are, of course, many challenges ahead. As noted, an important task is to make direct numerical comparisons between ensemble averages of the spin-diffusion model and the deterministic model in other test problems. In particular, it would be essential to look at problems in two and three dimensions with longer range potentials, where the mean field approximations is expected to be more consistent. This would allow for an assessment of the aforementioned approximations, and regions of validity could be determined. Alternatively, there may be ways to correct the mean field approximation.

Appendix A Details of Potentials

First, we give the details for the coefficients $J_{ij}$, $h_{i}$, and $f$. We have

[TABLE]

For the pair interactions in our test problem, we use a soft-core $4$-$2$ Lennard-Jones type potential:

[TABLE]

We introduce a cutoff of the above pair potential, to obtain a pair potential that is $C^{1}$ with finite support in $[0,r_{\rm c}]$:

[TABLE]

In addition to the pair potential, there is a quartic confining potential to prevent evaporation occurring as a rare event:

[TABLE]

The potential parameters are $\lambda=0.99$, $A_{++}=A_{--}=0.2$, $A_{+-}=0.18$, $r_{\operatorname{eq},++}=r_{\operatorname{eq},--}=2.6$, $r_{\operatorname{eq},+-}=2.55$. The cutoff radius is $r_{c}=10.5$, and the confining potential has parameter $r_{\mathrm{confine}}=2\cdot r_{\operatorname{eq},+-}=5.1$.

Appendix B Details of Numerical Methods for Deterministic

Models

In this Appendix, we provide further details about the numerical computation of the free energy, as well as other computational parameters for the simulations. Though the presentation is in terms of $\mathbf{X}$, $\mathbf{k}$, and $\mathbf{s}$, the simulations were performed in terms of $\mathbf{X}$, $\bm{\alpha}=\beta\mathbf{k}/2$ and $\mathbf{c}=(\mathbf{s}+1)/2$.

B.1 Calculation of free energy and its derivatives

Many of the expressions require computing expectations over a variational Gaussian distribution with parameters $\mathbf{X}$ and $\bm{\alpha}$. Since $\mathbb{E}^{\tilde{\nu}}\left[V(\mathbf{x},\bm{\sigma})\right]=\mathbb{E}^{\tilde{\mu}}\left[V(\mathbf{x},\mathbf{s})\right]$, for pair potential averages, like $\mathbb{E}^{\tilde{\mu}}\left[\phi(|x_{ij}|)\right]$, we can use the coordinate transformation from [24] which makes the integral one dimensional. We also integrate over a truncated region, $(X_{i}-R,X_{i}+R)$, in one dimension, with $R$ chosen sufficiently large that the domain truncation error is within a specified tolerance. This is done with a priori estimates involving the maximum value of the pair potential, and an assumed lower bound $\alpha_{\min}$ on the harmonic constants $\bm{\alpha}$.

The confining potential applies only to nearest-neighbor pairs, so we do not need a cutoff and therefore do not need to truncate the domain of integration for computation of its Gaussian average.

B.2 Computational Parameters

There are several computational parameters for the calculation of the free energy and its gradients. First, there is $R$, the radius of truncation described above. In our simulations, $R\approx 3.7$, which is a function of the a priori lower bound $\alpha_{\min}=5$, the upper bound $M=\max\{\phi_{++}(0),\phi_{+-}(0),\phi_{--}(0)\}$, and the truncation tolerance of $10^{-6}$. The absolute and relative tolerances for the adaptive quadrature have values $10^{-12}$ and $10^{-6}$ respectively.

There are also computational parameters for the minimization of $\mathcal{F}(\mathbf{c},\mathbf{X},\bm{\alpha})$ by gradient descent. The absolute and relative tolerances for the ODE solver are $10^{-6}$ and $10^{-3}$, respectively. The time step is chosen adaptively to meet these error tolerances. The gradient descent is stopped when $\max(\|\nabla_{\mathbf{X}}\mathcal{F}\|_{\infty},\|\nabla_{\bm{\alpha}}\mathcal{F}\|_{\infty})$ falls below the tolerance $10^{-7}$.

Finally, we have several computational parameters for evolution of the $\mathbf{c}$, $\mathbf{X}$, and $\bm{\alpha}$ by the transformed versions of (5.3a) – (5.3c). The absolute and relative tolerances for the ODE solver are $10^{-6}$ and $10^{-3}$, respectively. The evolution is halted when the infinity-norm of the gradient increases above the tolerance $10^{-6}$. The evolution of these equations continues once $\mathbf{X}$ and $\bm{\alpha}$ have been updated by gradient descent minimization. Finally, the prefactors on the evolution of $\mathbf{X}$ and $\bm{\alpha}$ are both $10^{10}$, which is equivalent to $\eta_{1}=10^{10}$ and $\eta_{2}=(2/\beta)^{2}10^{10}$ in (5.3b) and (5.3c).

Appendix C Details of Numerical Methods for Spin-Diffusion

Models

In this Appendix, we provide details of the numerical computation of spin-diffusion trajectories. We simulate the leading order process in Approximation 2, with exchange rates (4.30), and the State Space Approximation 4. The QSD averages of exchange rates and observables are estimated with time averages. These are computed by integrating the diffusion with preconditioned Metropolis Adjusted Langevin (MALA), [33]. While not strictly applicable to QSD sampling, in our low temperature regime, this efficiently produced consistent results. We obtained one hundred trajectories of the spin-diffusion process, from which we can estimate ensemble averages and confidence intervals.

C.1 Basin Sampling via Preconditioned MALA

For a given spin state, the basin is determined by finding the local energy minimizer from the last known minimizer. In this way, we capture Approximation 4. Energy minimization is performed using the Matlab function fminunc. Given the minimizer $\mathbf{x}_{\min}$ of the new spin state and basin, we use it to initialize preconditioned MALA with the inverse of the Hessian $H(\mathbf{x}_{\min})$ used as the preconditioner. The time step is chosen to maintain an acceptance rate between 65% and 85%. Estimates of the exchange rates and observables are obtained by averaging over the samples produced by preconditioned MALA.

C.2 Computational Parameters

The tolerances for minimization of $V(\mathbf{x})$ by the trust region method are $10^{-6}$ for the norm of the gradient, $10^{-6}$ for the relative change in the energy, and $10^{-8}$ for the norm of the change in $\mathbf{x}$.

The time step for MALA is $\Delta t=0.2$ for $N=32$ and $\Delta t=0.1$ for $N=500$. The interval of time averaging is $t_{\max}=1000$ for $N=32$ and $t_{\max}=500$ for $N=500$.