# Combining phase field crystal methods with a Cahn-Hilliard model for binary
alloys

Ananya Renuka Balakrishna  Department of Materials Science and Engineering,
Massachusetts Institute of Technology, Cambridge, MA02139, USA  W Craig Carter
Department of Materials Science and Engineering, Massachusetts Institute of
Technology, Cambridge, MA02139, USA

###### Abstract

During phase transitions certain properties of a material change, such as
composition field and lattice-symmetry distortions. These changes are
typically coupled, and affect the microstructures that form in materials.
Here, we propose a 2D theoretical framework that couples a Cahn-Hilliard (CH)
model describing the composition field of a material system, with a phase
field crystal (PFC) model describing its underlying microscopic
configurations. We couple the two continuum models via coordinate
transformation coefficients. We introduce the transformation coefficients in
the PFC method, to describe affine lattice deformations. These transformation
coefficients are modeled as functions of the composition field. Using this
coupled approach, we explore the effects of coarse-grained lattice symmetry
and distortions on a phase transition process. In this paper, we demonstrate
the working of the CH-PFC model through three representative examples: First,
we describe base cases with hexagonal and square lattice symmetries for two
composition fields. Next, we illustrate how the CH-PFC method interpolates
lattice symmetry across a diffuse composition phase boundary. Finally, we
compute a Cahn-Hilliard type of diffusion and model the accompanying changes
to lattice symmetry during a phase transition process.

## Introduction

Phase transitions in materials are typically accompanied by structural changes
to the lattice symmetry [1, 2, 3]. These changes include individual lattice
distortions [1, 3], grain rotations [2] and lattice defect formations [4]. All
of these influence the microstructures that form in a material [4, 5]. In this
paper, we introduce a modeling approach that couples the lattice symmetry with
the composition field. Two continuum methods, namely a Cahn-Hilliard model and
a phase field crystal model, are coupled and solved to describe a phase
transition process.

The lattice symmetry and grain orientations in microstructures are known to
affect physical properties of materials. For example, lithium diffusion in a
battery electrode induces lattice deformation [6, 7, 8, 9], which affects
lithium ion kinetics [6] and causes anisotropic expansion of electrodes [7].
Likewise, in a paraelectric to ferroelectric phase transitions, lattices
transform from centrosymmetric to other point groups lacking an inversion
centre. This transformation introduces stress-free spontaneous strains in the
ferroelectric system [1]. At present, theoretical models like phase field
methods describe complex microstructures in electrode/ferroelectric systems as
a function of the composition field (lithium-ion concentration, temperature or
polarization) [10, 11, 12, 13]. The Kobayashi-Warren-Carter phase field model
[9] further accounts for crystallographic misorientation at grain boundaries
during a phase transition process. While these modeling approaches provide
insights on the position of phase and grain boundaries, they only account for
grain orientations as an empirical parameter [9]. The current phase field
approaches do not allow for lattices to distort independently. Consequently,
the local strain fields arising from individual lattice distortions and the
presence of defects in a material system are not explored.

Alternatively, a phase field crystal (PFC) method proposed by Elder and Grant
[14, 15] describes atomistic details of material systems with periodic
solutions. This modeling technique describes coarse-grained symmetry of a
periodic system [16, 17], and is computed at faster time scales than the
molecular dynamics simulations [18]. The PFC model has been applied to explore
lattice defects in graphene [19] and nucleation problems in colloidal systems
[20]. Binary alloy models [21, 22, 23, 24], an extension to the PFC formalism,
was demonstrated to describe solidification [23, 25], crystallization [26, 27]
and phase segregation processes [28, 29]. The PFC approach is a useful tool
for multiscale modeling to describe the lattice symmetry of a material system.

In the current work, we combine the phase field crystal methods with a Cahn-
Hilliard model in a 2D theoretical framework to model a phase transition
process. The modeling approach couples two field parameters of a model system,
namely the composition field and the coarse-grained lattice-symmetry
distortions. The Cahn-Hilliard equation describes microstructures with a
composition order-parameter field. The phase field crystal equation models a
coarse-grained representation of lattice symmetry with peak density field as
its order parameter. In the PFC equation, we introduce coordinate
transformation coefficients to relate lattice symmetries in 2D point groups
via affine transformations. These transformation coefficients are coupled with
the composition field and influence the underlying lattice symmetry of a
material system. As the composition field evolves following the Cahn-Hilliard
equation, the transformation coefficients are updated in the PFC model. The
PFC model computes the equilibrium lattice arrangements of the material system
during composition evolution. Here, an assumption is that the dynamics of the
PFC model is fast relative to the composition field dynamics. Using this
coupled approach, we model the structural evolution of lattice distortions and
defects during a phase transition process.

In this paper, we investigate the nature of the coupled CH-PFC methods by
modeling three representative examples. First, we consider base cases to
understand how transformation coefficients affect the coarse-grained lattice
symmetries of the material system. Here, we stabilize hexagonal and square
symmetries as representative lattice structures corresponding to two
composition-field values. Second, we extend these base cases to investigate
how the model interpolates the peak density field across a diffuse composition
phase boundary. We model a representative binary alloy with hexagonal and
square symmetry phases and explore the lattice distortions across diffuse
interfaces. Finally, we model the composition field in the binary alloy to
follow a Cahn-Hilliard type of diffusion and study the accompanying
equilibrium lattice arrangements described by the peak density field. The
simulations show lattice distortions at coherent interfaces and demonstrate
structural evolution of lattice arrangements during a phase transition.

## Coupled Cahn-Hilliard – phase-field crystal model

The aim is to couple the Cahn-Hilliard (CH) and the phase field crystal (PFC)
methods, to explore structural changes to lattice symmetry during diffusion
induced phase transition. In this section, we first introduce the two
continuum models and explain how these methods are coupled in a 2D theoretical
framework. Next, we describe the evolution of the two order parameters, namely
the composition field and the peak density field. Finally, we discuss the
numerical procedure followed to compute the coupled CH-PFC methods.

The first model is a Cahn-Hilliard method that describes the continuum
composition field of a model system. This method utilizes a double-well free-
energy function in terms of a composition field, $\overline{c}$, which is its
order parameter. The second model is a phase field crystal (PFC) method that
describes the coarse-grained lattice symmetry for the model system, and
statistically illustrates lattice orientation, distortion and defect density.
This approach describes a free energy functional that is minimized by a
spatially-periodic order parameter, $\phi$. In the current work, we couple the
two models by using the composition field to influence the underlying lattice
symmetry of the model system. The composition is not coupled to the peak
density field $\phi$ via a homogeneous free energy, but rather as the
coordinate transformation coefficients of the composition-dependent Laplacian
$\nabla_{c}^{2}$, relative to a Cartesian basis. That is, each of the 5
Bravais lattices in the 2-dimensional space are stabilized by computing the
Laplace operator in a transformed space on a coordinate plane. The
transformation coefficients which control lattice deformations are described
as functions of the composition field. These coefficients are updated during
the evolution of the composition field.

The total free energy functional for the CH-PFC model is given by:

$\displaystyle F$
$\displaystyle=\int[g(\overline{c})+\kappa|\nabla\overline{c}|^{2}+f(\phi)+\frac{\phi}{2}G(\nabla_{c}^{2})\phi]d\vec{r}$
$\displaystyle=\int[\frac{F_{0}}{(\overline{c}_{a}-\overline{c}_{b})^{4}}(\overline{c}-\overline{c}_{a})^{2}(\overline{c}-\overline{c}_{b})^{2}+\kappa|\nabla\overline{c}|^{2}$
$\displaystyle\hskip 2.84544pt\hskip 2.84544pt\;+a\Delta
T\frac{\phi^{2}}{2}+u\frac{\phi^{4}}{2}+\frac{\phi}{2}\hskip
2.84544pt(\lambda(q_{0}^{2}+\nabla_{c}^{2})^{2})\hskip 2.84544pt\phi]d\vec{r}$
(1)

Here, $g(\overline{c})$ and $f(\phi)$ describe the homogeneous energy
contributions from the Cahn-Hilliard and PFC equations respectively. The
composition gradient-energy coefficient is given by $\kappa$. The operator
$G(\nabla_{c}^{2})$ controls the coarse-grained lattice symmetry described by
the particle density field. This operator is modeled as a function of the
composition field and is discussed in detail later on in this section. The
coefficients, $F_{0}$ and $(\overline{c}_{a},\overline{c}_{b})$, correspond to
the energy barrier height and to the local equilibrium states of
$g(\overline{c})$ respectively. The parameter $a\Delta T$ , controls the
second-order phase transition of the PFC model. In this paper, we model
$a\Delta T$ as a constant to always describe a crystalline-solid state. The
parameters $\lambda,q_{0},u$, relate the PFC equation to the first-order peak
in an experimental structure factor. Further details on these coefficients are
explained in the work by Elder and Grant [15].

Before proceeding with the model description, we first normalize the free
energy functional:

$\displaystyle\mathcal{F}$ $\displaystyle=\frac{F}{F_{0}}$
$\displaystyle=\int\\{c^{2}(c-1)^{2}+|\nabla
c|^{2}+\gamma[\frac{\psi}{2}(r+(1+\nabla_{c}^{2})^{2})\psi+\frac{\psi^{4}}{4}]\\}d\vec{x}.$
(2)

The composition field, $\overline{c}$ is normalized as
$c=\frac{\overline{c}_{a}-\overline{c}}{\overline{c}_{a}-\overline{c}_{b}}$,
with local equilibrium states at $c=0$ and $c=1$. The dimensionless peak
density field, $\psi$ is given by $\psi=\phi\sqrt{\frac{u}{\lambda
q_{0}^{4}}}$. We set, $r=\frac{a\Delta T}{\lambda q_{0}^{4}}=-0.2$, as a
constant such that Eq. 2 always models a stable crystalline-solid phase for
the peak density field [15]. With $c=0$, the peak density field in Eq. 2
describes a hexagonal symmetry with a periodic spacing of
$\frac{4\pi}{q_{0}\sqrt{3}}$ at equilibrium. Note, $\frac{1}{q_{0}}$ is the
length scale of the PFC model and $\vec{x}=q_{0}\vec{r}$. The gradient energy
coefficient
$\kappa=\frac{F_{0}}{(\overline{c}_{a}-\overline{c}_{b})^{2}}\left(\frac{16\pi}{q_{0}\sqrt{3}}\right)^{2}$,
is numerically calibrated such that the width of the diffuse composition
interface spans over $\sim 4$ peaks described by the peak density field,
$\psi$. We introduce a constant, $\gamma=\frac{\lambda^{2}q_{0}^{5}}{uF_{0}}$
that relates the free energy normalizations of the Cahn-Hilliard and the PFC
model. For simulations in this paper, we set $\gamma=1$.

The composition-dependent Laplacian $\nabla_{c}^{2}$ in Eq. 1, introduces the
composition-lattice symmetry coupling. Here, the composition terms enter the
Laplacian via its coordinate transformation coefficients. The Laplace operator
is written in terms of its second partial derivatives:

$\nabla_{c}^{2}=(\mathrm{A}_{1,1}^{2}+\mathrm{A}_{1,2}^{2})\frac{\partial}{\partial
x^{2}}+\mathrm{A}_{2,2}^{2}\frac{\partial}{\partial
y^{2}}+2\mathrm{A}_{1,2}\mathrm{A}_{2,2}\frac{\partial}{\partial x\partial
y},$ (3)

where $\mathrm{A}_{k,l}$ are the coordinate transformation coefficients. These
coefficients are described as functions of the dimensionless composition field
$c$ and correspond to the elements of a $k\times l$ transformation matrix:

$\mathbf{A(\mathit{c})}=\left[{\begin{array}[]{cc}\alpha(c)&\frac{2\alpha(c)}{\sqrt{3}}\hskip
2.84544pt\mathrm{cos}[\theta(c)]-\frac{\alpha(c)}{\sqrt{3}}\\\
0&\frac{2\beta(c)}{\sqrt{3}}\hskip
2.84544pt\mathrm{sin}[\theta(c)]\end{array}}\right].$ (4)

The matrix $\mathbf{A(\mathit{c})}$, describes affine lattice transformations
using hexagonal symmetry as the reference structure [30], [31]. With $c=0$,
the transformation matrix is an identity matrix and Eq. 2 describes a
hexagonal symmetry in 2D [15]. In the current work, we choose the hexagonal
and square symmetries to represent phases with compositions $c=0$ and $c=1$
respectively. These symmetries are chosen to illustrate exaggerated symmetry
deformations during phase transition. The transformation coefficients in Eq.
4, $(X=\alpha,\beta,\theta)$, are modeled as linear functions of the
dimensionless composition field, $X(c)=X_{0}+c\Delta X$. We define $X_{0}$ to
be the transformation coefficients corresponding to the hexagonal lattice
$(\alpha_{0}=\beta_{0}=1,\theta_{0}=\frac{\pi}{3})$, and $\Delta X$ is the
deformation required to transform the lattice with a hexagonal symmetry to a
square symmetry $(\Delta\alpha=\Delta\beta=0,\Delta\theta=\frac{\pi}{6})$.
Note, in both the hexagonal and square lattice symmetries, the transformation
matrix encourages a periodic lattic- symmetry spacing of
$\frac{4\pi}{q_{0}\sqrt{3}}$.

(a)

(b)

Figure 1: Schematic representations of the Cahn-Hilliard – phase field crystal
(CH-PFC) method. (a) The reference hexagonal symmetry (in blue) deforms to a
square symmetry (in red) under the transformation matrix described by Eq. 4,
with composition field $c=1$. (b) The lattice symmetry transforms from a
square to a hexagonal symmetry as a function of the composition field. The
dashed quadrilaterals across the diffuse phase-boundary illustrate the
intermediate lattice symmetries.

Fig. 1 shows a schematic illustration of the Cahn-Hilliard – phase field
crystal concept. In Fig. 1a, the transformation matrix describes lattice
symmetry as a function of the composition field. For $c=0$, the transformation
matrix $\mathbf{A}(c=0)$ is an identity matrix, which describes the
composition-dependent Laplacian $\nabla_{c}^{2}$ (in Eq. 3) in an isotropic
coordinate space. With $\mathbf{A}(c=0)$ the CH-PFC model stabilizes a
hexagonal lattice symmetry at equilibrium, see blue hexagonal symmetry in Fig.
1a. However, for a system with $c=1$, the transformation matrix
$\mathrm{\mathbf{A}}(c=1)$ introduces anisotropy in the transformation
coefficients (in Eq. 3), which models the composition-dependent Laplacian in a
transformed coordinate space. With $\mathrm{\mathbf{A}}(c=1)$ the CH-PFC model
results in a square symmetry at equilibrium, see the red square in Fig. 1a.
Next, Fig. 1b schematically illustrates how the CH-PFC model interpolates the
lattice symmetry across a diffuse phase boundary. Here, the transformation
matrix is locally defined in space as a function of the composition field. For
$0<c<1$, the transformation matrix $\mathbf{A}(0<c<1)$ interpolates the peak
density field to describe intermediate lattice symmetries between the square
and the hexagonal, see the dashed quadrilaterals in Fig. 1b.

Next, we describe the evolution of the two order parameters during phase
transition. Here, we assume that the elastic relaxation of the dimensionless
peak density field, $\psi$, is achieved instantaneously in comparison to the
evolution of the composition field. Consequently, we model
$\frac{\delta\mathcal{F}}{\delta\psi}\approx 0$ to be maintained throughout
the phase transition process.

The composition field evolves using a generalized Cahn-Hilliard equation:

$\displaystyle\frac{\partial c}{\partial\tau}$
$\displaystyle=-\nabla^{2}\frac{\delta\mathcal{F}}{\delta c}$
$\displaystyle=-\nabla^{2}(\gamma\frac{\psi}{2}\frac{\partial(\nabla_{c}^{4}+2\nabla_{c}^{2})\psi}{\partial
c}+4c^{3}-5c^{2}+2c-\nabla^{2}c).$ (5)

Here, $\gamma=1$ and $\tau$ is the dimensionless time variable
$\tau=t\frac{D}{L^{2}}$. $D$ is the isotropic diffusion coefficient in Eq. 5
and $L$ is the size of the simulation grid. The variational derivative in Eq.
5, produces coupled terms connecting the peak density field and the
composition field. In Eq. 5, it is of interest to note the two types of
Laplace operators, $\nabla^{2}\mathrm{\thinspace
and\thinspace}\nabla_{c}^{2}$, respectively. The Laplace operator $\nabla^{2}$
is $\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$.
This Laplacian computes the Cahn-Hilliard diffusion isotropically. The
composition-dependent Laplacian $\nabla_{c}^{2}$ describes its partial
derivatives in a transformed-coordinate space, see Eq. 3. The transformation
coefficients are influenced by the local composition field values and
$\nabla_{c}^{2}$ computes the derivatives of $\psi$ in a transformed-
coordinate space. The propagation of the composition diffusion front given by
Eq. 5 is affected by both the coarse-grained lattice arrangements and the
local-composition of the model system. As the composition field evolves, the
transformation coefficients in the composition-dependent Laplacian
$\nabla_{c}^{2}$, are updated accordingly.

As the elastic relaxation is much faster than composition evolution, we
introduce a time-like fictive variable $n$ to compute
$\frac{\delta\mathcal{F}}{\delta\psi}\approx 0$. The variable $n$ is treated
as a rapidly changing parameter in comparison to the dimensionless time,
$\tau$. This variable $n$ is used as a relaxation parameter to approximate
equilibrium of $\psi$ at each $c(\tau)$:

$\displaystyle\frac{\partial\psi}{\partial n}$
$\displaystyle=-\frac{\delta\mathcal{F}}{\delta\psi}+\frac{1}{n_{x}n_{y}}\int\frac{\delta\mathcal{F}}{\delta\psi}d\vec{x}$
$\displaystyle=\\{-\gamma[(r+(1+\nabla_{c}^{2})^{2})\hskip
2.84544pt\psi+\psi^{3}]$ $\displaystyle\ \
+\frac{1}{n_{x}n_{y}}\int\gamma[(r+(1+\nabla_{c}^{2})^{2})\hskip
2.84544pt\psi+\psi^{3}]\\}d\vec{x}.$ (6)

Here, $n_{x}$ and $n_{y}$ are the sides of a rectangular simulation domain,
and $\gamma=1$. Eq. 6 follows from the numerical scheme introduced by
Melenthin et al. [32] that allows equilibrium states to be attained faster in
comparison to the standard equation of motion of the PFC model [15]. Here,
$\psi$, is treated as a locally nonconserved order parameter, while the mass
conservation, $\int\psi d\vec{x}$, is ensured globally. Other approaches to
model faster dynamics for the peak density field can be found in the work by
Heinonen et al. [33]. Note, the variational derivative in Eq. 6 introduces
coupled composition-lattice symmetry terms. These coupled terms affect the
symmetry of the periodic system.

Eqs. 5 and 6 are computed using an Euler discretization scheme in a 2D finite-
difference framework. Simulation grids of size $n_{x}\times n_{y}$ are modeled
with periodic boundary conditions and with grid spacings of $\delta x=\delta
y=\frac{4\pi}{q_{0}6\sqrt{3}}$. At each grid point, the dimensionless
composition and peak density fields are represented in their discrete forms as
$c_{ij}$ and $\psi_{ij}$ respectively. The dimensionless composition time
derivative in Eq. 4 is computed at regular time steps of $\Delta\tau$, to
track the evolving composition field. At each time step, $\tau+\Delta\tau$,
the transformation coefficients of the Laplace operator $\nabla_{c}^{2}$ , are
updated to correspond with the evolving composition field. Next, the
equilibrium lattice symmetry at time $\tau+\Delta\tau$, is identified by
maintaining $\frac{\delta\mathcal{F}}{\delta\psi}\approx 0$. This general
numerical procedure is iterated. In other work, we apply the CH-PFC method to
model Li-ion diffusion in electrode materials [34].

## CH-PFC simulations

In this section we investigate the nature of the CH-PFC methods by simulating
a few representative examples. First, we explore how the transformation
coefficients stabilize hexagonal and square symmetries as a function of the
composition field. Using the hexagonal and square symmetries as base cases, we
next model a representative binary alloy with diffuse interfaces. Here, we
study how the model interpolates the peak density field across a diffuse phase
boundary. Finally, we simulate a Cahn-Hilliard type of diffusion for the
composition field and model the accompanying structural changes to the
underlying lattice symmetry during a phase transition.

### Lattice symmetry

At first, we describe two representative systems (not necessarily a physical
system) with homogeneous composition fields, $c_{ij}=0$ and $c_{ij}=1,$
respectively. The composition fields are treated as fixed. Using the
composition fields as input, we compute the peak density fields for the
periodic systems. These representative systems will be generalized
subsequently in the following subsections. Note, the peak density field is
rapidly evolving with reference to the composition field dynamics, and is
modeled with a fictive time in the subsequent computations, see Eq. 6.

Two simulation grids of size $100\times 100$ are modeled with periodic
boundary conditions. The transformation matrices at each grid point, for the
two representative systems with $c_{ij}=0$ and $c_{ij}=1$ are computed
following Eq. 4:

$\displaystyle\mathbf{\mathbf{A}_{\mathbf{\mathrm{H}}}=A}(c_{ij}=0)$
$\displaystyle=\left[\begin{array}[]{cc}1&0\\\ 0&1\end{array}\right],$ (9)
$\displaystyle\mathbf{A}_{\mathbf{\mathrm{S}}}=\mathbf{A}(c_{ij}=1)$
$\displaystyle=\left[\begin{array}[]{cc}1&-1/\sqrt{3}\\\
0&2/\sqrt{3}\end{array}\right].$ (12)

Matrices, $\mathbf{A}_{\mathbf{\mathrm{H}}}$ and
$\mathbf{A}_{\mathbf{\mathrm{S}}}$ describe the transformation coefficients to
model the hexagonal and square lattice symmetries respectively. Note, the
determinant of the matrices in Eq. 12 are
$\mathrm{det}(\mathbf{A}_{\mathbf{\mathrm{H}}})=1$ and
$\mathrm{det}(\mathbf{A}_{\mathbf{\mathrm{S}}})=1.15$ respectively. The
difference in the determinants
$\mathrm{det}(\mathbf{A}_{\mathbf{\mathrm{S}}})-\mathrm{det}(\mathbf{A}_{\mathrm{H}})=0.15$,
indicates an area change between the square and hexagonal lattices. This is
because, in the current work we model hexagonal and square symmetries to
assume equal lattice spacing of $\frac{4\pi}{q_{0}\sqrt{3}}$. Therefore the
number density of peaks changes with lattice symmetry.

(a)

(b)

Figure 2: Evolution of the peak density fields in representative systems
forming (a) hexagonal and (b) square symmetries with homogeneous composition
fields $c_{ij}=0$ and $c_{ij}=1$ respectively. Subfigures illustrate the peak
density fields starting from a randomized initial state (far left), during
evolution (centre) and at the final equilibrium state (far right).

Using the transformation matrices in Eq. 12 we next compute the peak density
fields of the periodic systems. The simulation grids are initialized with
random peak-density field values, $-0.1\leq\psi_{ij}\leq 0.5$ – a condition
that we will refer to as the “random initial seed”. Starting from this random
state and average density, $\overline{\psi_{ij}}=0.2$, the evolution of the
peak density field, Eq. 6, is iterated until equilibrium is reached.

Fig. 2 shows the evolution of density fields from randomized initial states,
for the two homogeneous composition fields, $c_{ij}=0$ and $c_{ij}=1$,
respectively. During evolution, individual grains with hexagonal and square
lattice symmetries nucleate in Fig. 2a and Fig. 2b respectively. Note, grains
of different sizes and lattice orientations form during a CH-PFC simulation,
see ’During evolution’ in Fig. 2. At the grain boundaries, lattice symmetries
distort to form coherent interfaces. At equilibrium, individual grains arrange
to minimize lattice misfits at the grain boundaries. A coarse-grained
representation of hexagonal and square lattice symmetries are formed in Fig.
2a and Fig. 2b respectively.

Fig. 2 shows the formation of multiple grains in homogeneous composition
fields and identifies the position/orientation of the grain boundaries in the
model system. In Fig. 2a, the density peaks that model the hexagonal symmetry
are of circular shape. However, for the square symmetry in Fig. 2b, the
density peaks are ellipsoidal in shape. This difference in the density peak
shapes is explained from the use of transformation matrices
$\mathbf{A_{\mathrm{H}}}$ and $\mathbf{A_{\mathrm{S}}}$ in Eq. 12. The
transformation matrix for hexagonal symmetry, $\mathbf{A_{\mathrm{H}}}$
describes an isotropic composition-dependent Laplacian, $\nabla_{c}^{2}$. This
computes the density peaks to be of circular shape. While, the transformation
matrix for a square symmetry, $\mathbf{A_{\mathrm{S}}}$ introduces
transformation coefficients in the composition-dependent Laplace operator, see
Eq. 12 and Eq. 3. These transformation coefficients shear the density peaks to
an ellipsoidal shape. Similar ellipsoidal density peaks are observed in the
anisotropic PFC simulations [23, 35]. Furthermore, the density peaks near
grain boundaries in both Fig. 2a and Fig. 2b, appear smeared and deviate from
the regular ellipsoidal/circular shapes. Here, an interpretation is that the
smeared appearance indicates lattice distortion at the interfaces to maintain
coherency between neighboring grains.

### Diffuse interface

Next, we investigate the model behaviour to interpolate the peak density field
across a diffuse interface in a representative binary alloy. Here, the
hexagonal and square lattice symmetries at compositions $c_{ij}=0$ and
$c_{ij}=1$ are used as base cases, and correspond to the two phases of the
binary alloy. A representative binary alloy with diffuse phase boundaries is
modeled and its composition field is treated to be fixed. The equilibrium
lattice symmetry for this system with heterogeneous composition field is
computed.

A periodic simulation grid of size $200\times 30$ is modeled. Here, two phases
with $c_{ij}=0$ and $c_{ij}=1$ separated by a sharp interface is assumed in
the initial state:

$c_{ij}=\left\\{\begin{array}[]{c}1\\\
0\end{array}\right.for\begin{array}[]{c}i<20,i>180\\\ 20\leq i\leq
180\end{array}.$ (13)

Next, the composition field is evolved following Eq. 5, without any influence
from the peak density field. That is, $\frac{\partial
c}{\partial\tau}=-\nabla^{2}\frac{\delta\mathcal{F}(c,\psi=0)}{\delta c}$. The
composition time derivative is iterated until the phase boundary begins to
smooth and is then held fixed. This is to explore the coupling of the fast
kinetics of $\psi$ for a single interation of $c$. Fig. 3(a) illustrates the
composition of a binary alloy with diffuse phase boundaries. Fig. 3(b) shows
the composition variation across the simulation grid at $j=15$.

Figure 3: (a) The composition field of a representative binary alloy. (b)
Variation of composition across AA in the simulation grid at $j=15$. (c) The
equilibrium peak density field describing the underlying symmetry of the
binary alloy corresponding to the composition field in subfigure 3(a). Square
and hexagonal symmetries are described for phases with $c_{ij}=1$ and
$c_{ij}=0$ respectively. The width of the composition phase boundary
(illustrated by vertical dashed-lines) is numerically calibrated to span
across $\sim 4$ peaks.

Following Eq. 3-4, the transformation matrix, $\mathbf{A}(c_{ij})$, is next
computed with $c_{ij}$ describing the discrete composition field shown in Fig.
3(a):

$\mathbf{A\mathrm{(\mathit{c_{ij}\mathrm{)}}}}=\left[{\begin{array}[]{cc}1&\frac{2}{\sqrt{3}}\hskip
2.84544pt\mathrm{cos}[\frac{\pi}{3}+\frac{\pi}{6}c_{ij}]-\frac{1}{\sqrt{3}}\\\
0&\frac{2}{\sqrt{3}}\hskip
2.84544pt\mathrm{sin}[\frac{\pi}{3}+\frac{\pi}{6}c_{ij}]\end{array}}\right].$
(14)

Here, $\mathbf{A}(c_{ij})$ defines the transformed space for the composition-
dependent Laplace operator at each grid point. Using this transformation
matrix as an input, the equilibrium peak density field is next computed, Eq.
6.

To model the lattice symmetry of the binary alloy shown in Fig. 3(a), the
simulation grid is initialized with random peak density field values,
$-0.1\leq\psi_{ij}\leq 0.5$. Using $\mathbf{A}(c_{ij})$ from Eq. 14, the
evolution of the peak density field, Eq. 6, is iterated to find the
equilibrium lattice-symmetry for the model system. Fig. 3(c) shows the
equilibrium lattice-arrangements described for the heterogeneous composition
field (shown in Fig. 3(a)). Lattices with square symmetry are stabilized in
the phase with $c_{ij}=1$, and hexagonal symmetry is observed in the phase
with $c_{ij}=0$. At the phase boundaries, $0<c_{ij}<1$, the coupled CH-PFC
model describes a coarse-grained representation of deformed lattices. Here,
the density peaks are smeared to illustrate the lattice distortion at the
phase boundaries, see Fig. 3(c). Note, the composition phase boundary is
numerically calibrated to span over $\sim 4$ density peaks (about 25 grid
spacings). Fig. 3 provides an atomistic insight into the coarse-grained
lattice arrangements across a diffuse phase boundary.

### Phase transition

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Phase transition showing structural evolution of lattices from a
hexagonal (in blue), to a square symmetry (in red). Subfigures illustrate
lattice transformations as a function of the dimensionless composition field,
$c_{ij}$. The green arrows indicate the orientation of the square symmetry
during phase transition. The dashed line in subfigures 4(e-f) indicate a grain
boundary in the square symmetry phase.

Up to this point, we only modeled the microscopic configurations at fixed
compositions. However, to model phase transition with microscopic insights on
the coarse-grained lattice symmetry, we need to simulate the evolution of the
composition field. The binary alloy in Fig. 3 is considered as the initial
state, and we next extend the simulation to describe the propagation of the
diffusion front. A representative Cahn-Hilliard type of diffusion for the
composition field is modeled. During the phase transition, the equilibrium
lattice arrangements of the underlying system is computed. An assumption made
in this simulation is that the dynamics of elastic relaxation (equilibrating
the peak density field) is several times faster than the diffusion of the
composition field. Using this CH-PFC approach we investigate how composition
field influences the lattice arrangements in a model system during phase
transitions.

Taking as an initial state, the lattice arrangements described for a binary
alloy from Fig. 3, the phase transition is modeled by allowing composition to
diffuse into the simulation domain. The composition field is held fixed at
$c_{ij}=1$, for $i<20$ and $i>180$ throughout the simulation. This boundary
condition is a proxy for having a consistent composition reservoir. The
composition field on the remaining part of the simulation grid, $20\leq i\leq
180$, is allowed to evolve with time. The composition time derivative, Eq. 5,
is iterated from $\tau=0$ to $\tau=2500$, in dimensionless time intervals of
$\Delta\tau=25$. Note, the composition evolution at $\tau$, receives input
from the equilibrium peak density field calculated for the $(\tau-1)$ time
step. The composition field is tracked as
$c_{ij}(\tau+\Delta\tau)=c_{ij}(\tau)+\Delta\tau\frac{\partial
c_{ij}}{\partial\tau}$, until a homogenous phase is obtained.

The composition field at each evolution step, $c_{ij}(\tau)$, is used as an
input to compute the transformation matrix in Eq. 14. At a given time step,
$\tau$, the transformation matrix $\mathbf{A}(c_{ij}(\tau))$, is used to
calculate the equilibrium peak-density-field following Eq. 6. The composition
and peak density fields are iterated until the phase transition is complete.

Fig. 4 shows the structural evolution of the coarse-grained lattice
arrangements during the phase transition. At the intial state $\tau=0$, the
coarse-grained lattice symmetry for the heterogeneous composition,
$c_{ij}(\tau=0)$ is described, see Fig. 4a. Here, two coherent phases with
square and hexagonal symmetries are formed in domains with $c_{ij}=1$ and
$c_{ij}=0$ respectively. Note in Fig. 4a, the edges of the square lattices are
mostly aligned with the axes of the simulation grid. A pair of green arrows in
Fig. 4a illustrates the orientation of square lattices in the simulation grid.
Across the diffuse phase boundary, hexagonal and square lattices are distorted
to maintain coherency, see Fig. 4a. Next, in Fig. 4(b-e), as the composition
field diffuses into the simulation domain, the hexagonal lattice symmetry is
transformed to a square symmetry.

In Fig. 4b, the phase with square symmetry occupies $\sim 50\%$ of the
simulation grid. Here, it is interesting to note that square lattices begin to
rotate uniformly as the diffusion front propagates through the simulation
grid. In Fig. 4(c-d), the square lattice symmetry is observed to rotate
further (e.g., orientation of the green arrows in Fig. 4(c-d)). We interpret
that the square lattices rotate to maintain coherency with the neighboring
hexagonal phase. Note, the periodic boundary conditions on the simulation
domain further enforce an additional strain on the peak density field. This is
discussed in detail in the next section of this paper. In Fig. 4e, a grain
boundary (as indicated by the dashed line) is formed in the square symmetry
phase. This grain boundary migrates in the square symmetry phase and remains
in the homogeneous phase, see Fig. 4f. At $\tau=2500$, the phase transition is
complete with a homogenous composition field and a phase with square symmetry
is described at equilibrium, see Fig. 4f.

## Discussion of the CH-PFC model

The coupled Cahn-Hilliard – phase field crystal model provides a theoretical
framework to describe continuum phase transition with microscopic insights.
There are several issues we feel remain to be clarified in interpreting the
simulations. Among these issues are three questions: Do the peaks in the CH-
PFC simulations represent atomic sites or illustrate the underlying lattice
symmetry? Are the total number of peaks in a simulation grid conserved? In Eq.
1 why was the composition field coupled with the peak density field only via
the Laplace operator? In this section, we discuss these key details of the
coupled CH-PFC model and explore potential further work.

First, the peak density field in CH-PFC simulation describes the coarse-
grained lattice symmetry of the underlying atomic arrangements. Individual
peaks do not represent atomic sites, however the arrangement of peaks
indicates the unit cell symmetry of the model system. Similarly, a grain
boundary in a CH-PFC simulation is a coarse-grained approximation of the
underlying lattice orientations, distortions and defects. Fig. 5 provides a
schematic illustration of the difference between atomic sites, peak positions
and coarse-grained lattice symmetry. In Fig. 5, the small-black dots indicate
atomic sites, which correspond to the deterministic positions of atoms in the
unit cell. The big-green dots highlight representative peak positions modeled
by a CH-PFC method. The dashed-red lines connecting the peaks in Fig. 5,
indicate an example of a coarse-grained lattice symmetry. In Fig. 5, the side
of the coarse-grained lattice is four times that of the unit cell. However, in
our CH-PFC simulations, the coarse-grained lattice is several times larger
than a unit cell.

Figure 5: Schematic illustration of the atomic sites, peak positions and
coarse-grained lattice symmetry in a square phase system. The small-black dots
represent atomic sites of the unit cell. The big-green dots schematically
indicate peak positions in the CH-PFC simulation. The dashed-red lines
connecting the peaks highlight an example of a coarse-grained lattice symmetry
of the underlying unit cells. Note, in this figure the coarse grained symmetry
is four times the unit cell size. However, in our CH-PFC simulations the
coarse-grained symmetry is multiple times larger than a unit cell.

Second, the number density of peaks in the CH-PFC simulations are not
necessarily conserved. Let us consider Fig. 2, where hexagonal and square
symmetries are described on identitical computational grids of size $100\times
100$. The total number of peaks in both these symmetry systems are not
necessarily the same. This can be explained from two reasons: First, the
transformation matrices in Eq. 12, $\mathbf{A_{\mathrm{H}}}$ and
$\mathbf{A_{\mathrm{S}}}$, describe lattice symmetries with an area difference
$(\sim 0.15)$. Second, the periodic boundary conditions on the computational
grid enforces a strain on the peak density field. That is, on an infinitely
large grid size, the peak density field would assume the fundamental length
scale of $\frac{4\pi}{q_{0}\sqrt{3}}$ specified by the transformation
coefficients in Eq. 3. However, by modeling this density field on a periodic
grid with finite dimensions, we strain the lattice symmetry spacing and force
the peak density field to satisfy periodic boundary conditions. To minimize
these imposed strains and to simultaneously maintain periodicity, the CH-PFC
model introduces (or removes) peaks to (or from) the simulation grid. For
future applications of the CH-PFC model, the computational grid size is to be
calibrated to correspond to the closest fundamental length scale of the peak
density field. Alternatively, numerical correction terms to Eq. 6 to conserve
the number of peaks can be used [36].

Finally, in Eq. (1), the composition and the peak density fields are coupled
only via the Laplace operator. That is, the coarse-grained lattice symmetry
described by Eq. 1 is solely determined by the coordinate transformation
coefficients of the Laplace operator. These transformation coefficients (which
are functions of the composition field) describe lattice transformations with
hexagonal symmetry as the reference structure. In this paper, we assumed the
ideal free energy contribution from other non-linear terms
$(\psi^{2},\psi^{4})$ in Eq. 1, to be independent of the composition field for
a couple of reasons: First, this assumption allows the CH-PFC model to
stabilize a reference hexagonal lattice symmetry for composition field,
$c_{ij}=0$. Second, Eq. 1 will always describe a crystalline/ordered state for
the model system. This is because the driving force for the peak density field
towards the disordered state (controlled by term $\psi^{4}$) is not a function
of the composition field.

## Summary

We introduced a 2D theoretical framework, which combined a Cahn-Hilliard (CH)
model and a phase field crystal (PFC) model, to describe a phase transition
process. In this CH-PFC method, the composition field was coupled to the
coarse-grained lattice symmetry (peak density field) of the periodic system.
The CH-PFC modeling approach captured the effects of microscopic
configurations, such as lattice orientations, distortions and presence of
defects, on the phase-transition process. Furthermore, the model described the
structural evolution of the coarse-grained lattice symmetry during a phase
change.

Using the CH-PFC approach, we stabilized representative lattice symmetries
(hexagonal and square) as a function of the composition field. Here, we found
that multiple grains formed in a single phase, and identified the position and
orientation of grain boundaries. Next, in a binary alloy, we described the
coarse-grained distortion of lattice symmetry across a diffuse phase boundary.
Finally, we modelled a representative phase transition process – here, the CH-
PFC simulations modeled grain rotations and grain boundary migrations during
phase change.

## Acknowledgements

A.R.B acknowledges the support of the Lindemann trust fellowship. The authors
would like to acknowledge the support by the grant DE-SC0002633 funded by the
U.S. Department of Energy, Office of Science, in carrying out this work.
Further, the authors wish to thank Dr. Rachel Zucker for useful discussions on
phase field crystal modeling methods.

## References

  * [1] Glass A. M. Lines M.E.  Principles and applications of ferroelectrics and related materials.  Oxford university press, 1977.
  * [2] R.D. James.  Displacive phase transformations in solids.  J. Mech. Phys. Solids, 34:359–394, 1986.
  * [3] N. Meethong, Y. H. Kao, M. Tang, H. Y. Huang, W. C. Carter, and Y.-M. Chiang.  Electrochemically Induced Phase Transformation in Nanoscale Olivines Li1-xMPO4 (M = Fe, Mn).  Chem. Mater., 20:6189–6198, 2008.
  * [4] A. S. Levanyuk, A. P., & Sigov.  Defects and structural phase transitions.  Gordon and Breach Science Publishers, 1988.
  * [5] R. D. Song, Y., Chen, X., Dabade, V., Shield, T. W., & James.  Enhanced reversibility and unusual microstructure of a phase-transforming material.  Nature, 502(7469):85–88, 2013.
  * [6] A. Nie, L.-Y. Gan, Y. Cheng, Q. Li, Y. Yuan, F. Mashayek, H. Wang, R. Klie, U. Schwingenschlogl, and R. Shahbazian-Yassar.  Twin boundary-assisted lithium ion transport.  Nano Lett., 15:610–615, 2014.
  * [7] Y. Yuan, A. Nie, G.M. Odegard, R. Xu, D. Zhou, S. Santhanagopalan, K. He, H. Asayesh-Ardakani, D.D. Meng, R.F. Klie, and C. Johnson.  Asynchronous crystal cell expansion during lithiation of K+-stabilized $\alpha$-MnO2.  Nano Lett., 15:2998–3007, 2015.
  * [8] M. Tang, W.C. Carter, and R.M. Cannon.  Diffuse interface model for structural transitions of grain boundaries.  Phys. Rev. B, 73:024102, 2006.
  * [9] J.A. Warren, R. Kobayashi, A.E. Lobkovsky, and W.C. Carter.  Extending phase field models of solidification to polycrystalline materials.  Acta Mater., 51:6035–6058, 2003.
  * [10] M. Tang, H.Y. Huang, N. Meethong, Y.H. Kao, W.C. Carter, and Y.M. Chiang.  Model for the particle size, overpotential, and strain dependence of phase transition pathways in storage electrodes: application to nanoscale olivines.  Chem. Mater., 21:1557–1571, 2009.
  * [11] D.A. Cogswell and M.Z. Bazant.  Coherency strain and the kinetics of phase separation in LiFePO4 nanoparticles.  ACS Nano, 6:2215–2225, 2012.
  * [12] L.Q. Chen.  Phase-field models for microstructure evolution.  Annu. Rev. Mater. Res., 32:113–140, 2002.
  * [13] A.R. Balakrishna, J.E. Huber, and I. Münch.  Nanoscale periodic domain patterns in tetragonal ferroelectrics: A phase-field study.  Phys. Rev. B, 93:174120, 2016.
  * [14] K.R. Elder, M. Katakowski, M. Haataja, and M. Grant.  Modeling elasticity in crystal growth.  Phys. Rev. Lett., 88:245701, 2002.
  * [15] K.R. Elder and M. Grant.  Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals.  Phys. Rev. E, 70:051605, 2004.
  * [16] H. Emmerich, H. Löwen, R. Wittkowski, T. Gruhn, G.I. Tóth, G. Tegze, and L. Gránásy.  Phase-field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: an overview.  Adv. Phys., 61:665–743, 2012.
  * [17] N. Provatas, J.A. Dantzig, B. Athreya, P. Chan, P. Stefanovic, N. Goldenfeld, and K.R. Elder.  Using the phase-field crystal method in the multi-scale modeling of microstructure evolution.  JOM J. Miner. Met. Mater. Soc., 59:83–90, 2007.
  * [18] P.F. Tupper and M. Grant.  Phase field crystals as a coarse-graining in time of molecular dynamics.  EPL (Europhysics Lett., 81:40007, 2008.
  * [19] M. Seymour and N. Provatas.  Structural phase field crystal approach for modeling graphene and other two-dimensional structures.  Phys. Rev. B, 93:035447, 2016.
  * [20] T. Gránásy, L., Tegze, G., Tóth, G.I. and Pusztai.  Phase-field crystal modelling of crystal nucleation, heteroepitaxy and patterning.,.  Philos. Mag., 91:123–149, 2011.
  * [21] M. Elder, K.R., Provatas, N., Berry, J., Stefanovic, P. and Grant.  Phase-field crystal modeling and classical density functional theory of freezing.  Phys. Rev. B, 75:064107, 2007.
  * [22] N. Elder, K.R., Huang, Z.F. and Provatas.  Amplitude expansion of the binary phase-field-crystal model.  Phys. Rev. E, 81:011602, 2010.
  * [23] H. Kundin, J., Choudhary, M.A. and Emmerich.  Bridging the phase-field and phase-field crystal approaches for anisotropic material systems.  Eur. Phys. J. Spec. Top., 223:363–372, 2014.
  * [24] P.W. Alster, E., Elder, K.R., Hoyt, J.J. and Voorhees.  Phase-field-crystal model for ordered crystals.  Phys. Rev. E, 95:022105, 2017.
  * [25] T. Heinonen, V., Achim, C.V., Elder, K.R., Buyukdagli, S. and Ala-Nissila.  Phase-field-crystal models and mechanical equilibrium.  Phys. Rev. E, 89:032411, 2014.
  * [26] M. Elder, K.R., Provatas, N., Berry, J., Stefanovic, P. and Grant.  Phase-field crystal modeling and classical density functional theory of freezing.  Phys. Rev. B, 75:064107, 2007.
  * [27] J. Zhang, W. and Mi.  Phase field crystal modelling of the order-to-disordered atomistic structure transition of metallic glasses.  IOP Conf. Ser. Mater. Sci. Eng., 117:012056, 2016.
  * [28] W. Tao, Y., Zheng, C., Jing, Z., Wei-Ping, D. and Lin.  Effect of grain boundary on spinodal decomposition using the phase field crystal method.  Chinese Phys. Lett., 29:078103, 2012.
  * [29] M. Greenwood, M., Sinclair, C. and Militzer.  Phase field crystal model of solute drag.  Acta Mater., 60:5752–5761, 2012.
  * [30] R.D. Ball, J.M. and James.  Fine phase mixtures as minimizers of energy.  Arch. Ration. Mech. Anal., 100:13–52, 1987.
  * [31] K. Bhattacharya.  Microstructure of martensite: why it forms and how it gives rise to the shape-memory effect.  Oxford university press, 2003.
  * [32] M. Mellenthin, J., Karma, A. and Plapp.  Phase-field crystal study of grain-boundary premelting.  Phys. Rev. B, 78:184110, 2008.
  * [33] T. Heinonen, V., Achim, C.V., Kosterlitz, J.M., Ying, S.C., Lowengrub, J. and Ala-Nissila.  Consistent hydrodynamics for phase field crystals.  Phys. Rev. Lett., 116:024303, 2016.
  * [34] A.R. Balakrishna and W.C. Carter.  Modeling phase transition in nanoscale electrodes using the coupled Cahn-Hilliard - phase field crystal methods.  Prep., 2017.
  * [35] H. Prieler, R., Hubert, J., Li, D., Verleye, B., Haberkern, R. and Emmerich.  An anisotropic phase-field crystal model for heterogeneous nucleation of ellipsoidal colloids.  J. Phys. Condens. Matter, 21:464110, 2009.
  * [36] A. Aland, S., Hatzikirou, H., Lowengrub, J. and Voigt.  A mechanistic collective cell model for epithelial colony growth and contact inhibition.  Biophys. J., 109:1347–1357, 2015.


# Motion Planning Transformers: One Model to Plan Them All

Jacob J. Johnson, Linjun Li, Ahmed H. Qureshi, and Michael C. Yip
University of California San Diego
La Jolla, CA 92093
<EMAIL_ADDRESS>

###### Abstract

Transformers have become the powerhouse of natural language processing and
recently found use in computer vision tasks. Their effective use of attention
can be used in other contexts as well, and in this paper, we propose a
transformer-based approach for efficiently solving the complex motion planning
problems. Traditional neural network-based motion planning uses convolutional
networks to encode the planning space, but these methods are limited to fixed
map sizes, which is often not realistic in the real-world. Our approach first
identifies regions on the map using transformers to provide attention to map
areas likely to include the best path, and then applies local planners to
generate the final collision-free path. We validate our method on a variety of
randomly generated environments with different map sizes, demonstrating
reduction in planning complexity and achieving comparable accuracy to
traditional planners.

##  1 Introduction

Motion planning is an integral component of any navigation stack. It is an
age-old problem with more than a few decades of innovative solutions. By far
the most popular recent methods involve constructing trees or graphs in the
planning space by randomly sampling points [LaValle and James J. Kuffner,
2001, Kavraki et al., 1996] or discretizing the space into grids [Hart et al.,
1968].

Figure 1: Given a planning problem (left), we propose a hierarchical approach
to planning. In our method, the MPT proposes a region where a valid path might
exist (center). A sampling-based planner is used to find a collision-free path
(right) given the proposed region.

Although these methods come with theoretical guarantees, they do not scale
well to planning problems in larger planning spaces where there may be
significant obstacles. The problems come up frequently in the real world,
where map sizes may vary significantly, and obstacles and impedances may vary
across tasks and domains. Hence, in our ever-expanding pursuit of automation,
we need planning algorithms that can scale up while effectively keeping the
planning problem’s complexity and computation time in check.

In the robotics community, learning-based methods to solve planning problems
have been proposed to offer speed improvements to planning and typically
revolve around using convolution neural networks (CNN) and multi-layer
perceptrons (MLP) to predict where samples should be generated to construct
trees [Qureshi et al., 2019, Johnson et al., 2020]. Other works employ
different forms of latent representation such as gird cells [Banino et al.,
2018] and Sparse Graphical Memory (SGM) [Emmons et al., 2020] to generate a
plan. Chaplot et al. [2021] proposed one of the few works that use
transformers to learn a value function similar to Value Iteration Networks
[Tamar et al., 2016], but the extension of this method to larger maps still
need to be explored. Overall, the utilization of Transformers for motion
planning remains unexplored at large.

In this work, we propose a hierarchical approach using transformers. Our work
is motivated by the recent progress in expanding the language transformer
models for vision tasks and taking advantage of the power of self-attention
[Dosovitskiy et al., 2021, Liu et al., 2021]. Our model plans in a two-stage
approach (Fig. 1), first a transformer-based region proposal network
recommends a region where a potential path can exist. We refer to this module
as the Motion Planning Transformer (MPT). Next, a sampling-based planner tries
to plan a path given the proposed region. Unlike previous neural motion
planning methods, the objective of the MPT is to reduce the search space for
the underlying motion planners. Hence MPT can be integrated with any
traditional path planner.

We assess the performance of MPT augmented Sampling-based Motion Planners
(SMPs) [Karaman and Frazzoli, 2011, Gammell et al., 2014] on a synthetic and
real-world dataset comprising different map sizes in cluttered environments
for point-mass and non-holonomic car models. Our results indicate that MPT
augmented planners achieve sufficiently high accuracy in planning while
substantially reducing planning time and the vertices on the planning tree
compared to traditional methods. We also show that MPT maintains the same
qualities to maps of different sizes without fine-tuning or additional
training and generalizes to real-world environments. The diversity of
environments represented, from block obstacles to narrow passage mazes speak
to the potential general application of the MPT approach across multiple
environment types and domains.

##  2 Related Work

The most relevant work to our transformers-based region proposal network for
motion planning is perhaps the guided sampling-based motion planning methods.
They analytically or through learned heuristics determine a subset in robot
space that probably contains a path solution. For instance, [Qureshi and Ayaz,
2016, Tahir et al., 2018] employ Artificial Potential Fields (APF) within
sampling-based methods such as RRT* [Karaman and Frazzoli, 2011] and
Bidirectional RRT* [Qureshi and Ayaz, 2015] to guide a subset of random
samples towards promising regions that possibly contain an optimal path
solution. In contrast, Informed-RRT* (IRRT*) [Gammell et al., 2014] depends on
an initial path from an RRT* algorithm to compute an ellipsoidal region
probably containing an optimal path solution. However, in most planning
problems finding an initial path solution is itself challenging. In a similar
vein, Batch Informed Trees (BIT*) [Gammell et al., 2015] begins from an
elliptical region formed by a straight line path ignoring all obstacles and
incrementally expand it until an initial path solution is found. Once an
initial path is determined, it is further optimized by adapting the
precomputed ellipsoid and generating new samples within that space.

Another approach for region selection is by Zucker et al. [2008], which uses
REINFORCE algorithm [Williams, 1992] on discretized workspaces to guide
samples for underlying SMP methods. Similarly, Value Iteration Networks (VIN)
[Tamar et al., 2016] also discretizes the space and learns a value map to
guide path planning. Universal Planning Networks (UPN) [Srinivas et al., 2018]
extends VIN to continuous control spaces. Despite all advancements from APF-
based directed trees to learning value maps, these methods only consider
geometric planning, i.e., path planning under collision-avoidance constraints.
They are yet to be evaluated in practical problems with complex constraints
such as kinodynamic or non-holonomic constraints in robot manipulation and
autonomous car navigation tasks.

Although there exist random sampling-based approaches such as RRT* [Karaman
and Frazzoli, 2011, Arslan et al., 2017] and SST [Li et al., 2016] that
randomly explore the robot state space and satisfy advance constraints, they
suffer from the curse of dimensionality and take a considerable computational
time in cluttered spaces. Neural Motion Planning [Qureshi et al., 2020, Ichter
et al., 2018] has recently emerged as a promising tool for solving a wide
range of planning problems under various task constraints, ranging from non-
holonomic [Johnson et al., 2020, Li et al., 2021] to advanced manifold
kinematic constraints [Qureshi et al., 2020], with high computational speed.
These scalable methods learn sampling distributions from expert demonstrations
and, on execution, generate samples for an underlying planner forming a subset
that potentially contains a path solution. However, these approaches assume a
fixed size input environment map and often require redefining network
architectures and retraining for different map sizes. However, recent
developments in deep learning, primarily through Transformers [Dosovitskiy et
al., 2021, Liu et al., 2021], have provided us with ways to relax such
assumptions. Our proposed approach leverages these developments and introduces
a region proposal framework that can work with variable map sizes and enhances
underlying motion planners to solve complex problems in cluttered
environments.

##  3 Method

This section will go over our planning pipeline, which includes the MPT module
and the traditional SMPs used to generate the final plan. Our framework
proposes promising regions from large maps to the underlying planners for
quickly computing path solutions.

###  3.1 Motion Planning Transformer (MPT)

Figure 2: Overview of MPT Module. (Start from the top left and move clockwise)
The input map, and the start and goal encoded map are concatenated together
and passed through a FCN (Fully Connected Networks). The output of which is
reshaped and passed through the standard Transformer module. The classifier
predicts the probability of each of the token from the encoder output to
create the mask superimposed on the output map. The light green patches
throughout the pipeline represents how a patch $p\times p$ remains connected
with output of the classifier.

The MPT module, a region proposal network, is similar to that proposed in Ren
et al. [2015], but apart from having CNN and linear layers, we also include a
Transformer model. The Transformer block we use is similar to Dosovitskiy et
al. [2021], which is derived from Vaswani et al. [2017]. An overview of the
model is shown in Figure 2. For a given input costmap
$\bm{m_{c}}\in\mathbb{R}^{H\times W}$, start position
$\bm{x_{s}}\in\mathbb{R}^{2}$, and goal position $\bm{x_{g}}\in\mathbb{R}^{2}$
we encode the start and goal position onto an array of size $H\times W$. The
goal and start location are represented by patches of size $p\times p$ with
values -1 and 1 respectively. The encoded map,
$\bm{m_{e}}\in\mathbb{R}^{H\times W}$, is concatenated with $\bm{m_{c}}$ to
form the input, $\bm{m}\in\mathbb{R}^{H\times W\times 2}$, to the feature
extractor consisting of fully convolutional networks (FCN).

To extract the features from the encoder, we use a series of convolution,
ReLU, and MaxPool layers, which have a combined receptive area of $p\times p$.
We refer to this module as the feature extractor. The feature extractor is
similar to the CNN used by the hybrid architecture proposed in Dosovitskiy et
al. [2021], but instead of using a pre-trained network, we learn the filter
parameters from scratch. The feature extractor is similar to the convolutional
layers used in Faster R-CNN’s [Ren et al., 2015].

The output of the feature extractor is of size $H_{l}\times W_{l}\times d$,
where $H_{l}$ and $W_{l}$ is determined by the size of $\bm{m_{c}}$ and the
convolution layers, and $d$ is the latent dimension of the Transformer. The
output is reshaped to size $(H_{l}\cdot W_{l})\times d$ in order to be
compatible with the Transformer model. Each row vector is analogous to the
word token’s in NLP.

Each row of the flattened output of the feature extractor can be mapped to a
patch of size $p\times p$ on $\bm{m_{c}}$ centered at position
$\hat{x}\in\mathbb{R}^{2}$. We call all such points as anchor points. There
are $H_{l}\times W_{l}$ anchor points for a given input map. To encode the
positions of each patch, we add fixed sinusoidal position embeddings to the 1D
features similar to [Dosovitskiy et al., 2021, Vaswani et al., 2017]. Since
the positional encoding used is of a fixed size, there is an upper bound on
the size of the input map.

The modified patch embedding is passed through a series of multi-headed self-
attention (MSA), and multi-layer perceptron (MLP) blocks, identical to the
architecture used by Vaswani et al. [2017]. The MLP block consists of two
linear layers with RELU non-linearity. We also apply Dropout, Layer Norm and
residual connections after every MSA and MLP block. To save on memory, we add
gradient checkpoint’s [Chen et al., 2016] to all MSA blocks. Since the feature
extractor and the Transformer block are both invariant to the size of the
input map, this architecture is scalable to maps of any size.

Each of the output token from the transformer encoder is passed through a
classifier to identify regions that might represent a path. Similar to Ren et
al. [2015], this is implemented using a $1\times 1$ convolution layer. Once
all the tokens have been identified, a mask of dimension $H\times W$ is
generated which highlights the region with a potential path.

###  3.2 Path Planning

The masked map is used as the reference map for finding a path and checking
its feasibility. Any traditional or learning-based planners can be used to
solve the given problem. In this work, we use variations of the Rapidly
Exploring Random Trees (RRT) algorithm that guarantee optimality [Karaman and
Frazzoli, 2011, Li et al., 2016, Gammell et al., 2014] to find the path under
different settings. Note that MPT can enhance any existing path planner. It is
robust to different map sizes, making it ideal for real-world scenarios,
especially for robot navigation indoors (e.g., hospitals and homes) and
outdoors (e.g., streets, highways) environments, where there are areas of
variable sizes. MPT selects a subset in a larger map with fixed local
receptive fields. The selected local patches become the new input map for the
path planner, which would otherwise explore the entire space for pathfinding.

##  4 Experiments

We evaluate the planning capabilities of SMPs and Motion Planning Transformers
(MPT) aided SMPs. To test the efficiency of the planner, we train the model on
different types of environments of fixed size and evaluate on unseen maps of
various sizes. To evaluate the applicability of the planner to other robotic
systems, we train MPT on different types of robot. Lastly, we evaluate our
algorithm on real world map without training on any additional data.

###  4.1 Setup

####  4.1.1 Environments

To test the planning capabilities of our method, we evaluated the model on two
different kinds of synthetic environments. The first environment is called the
Random Forest, where 100 circular and square objects are randomly placed in
different positions and orientations on the map. The second environment is
called the Maze environment. A perfect maze is generated using the randomized
depth-first search. One characteristic of a perfect maze is that any start and
goal pairs on this map are reachable by a collision-free path. The distance
resolution of each map was set at 5cm per pixel.

####  4.1.2 Robot Models and SMP

We test our algorithm on two robotic systems. The first is a simple Point
Robot Model that can move in any direction in $\mathbb{R}^{2}$. Planning for
this robot is done in the $\mathbb{R}^{2}$ space using the RRT* [Arslan et
al., 2017] and Informed-RRT* (IRRT*) [Gammell et al., 2014] algorithms. The
other robot system is a Dubins Car Model, controlled using wheel velocity and
steering angle. The planning for this robot is done in the control space using
the Stable Sparse RRT (SST) algorithm [Li et al., 2016]. We used the
implementation from the Open Motion Planning Library (OMPL) [Şucan et al.,
2012] for all SMPs.

####  4.1.3 Dataset

We created a simulated dataset for the environments and robot models mentioned
above. We collected 25 paths for random start and goal pairs for 1750
environments using RRT* for both the Random Forest and Maze environments for
the point robot. For the Dubins Car Model, we collected 50 paths for random
start and goal pairs for 1000 environments using SST. All maps used for
training were of size 480$\times$480.

####  4.1.4 Training

The anchor points are classified as positive if they represent a path patch
and negative otherwise. All anchor points near the given path up to 0.7 meters
are considered positive, and the rest are considered background. For the
Dubins Car Model, we used the robot’s position alone to pick the positive
anchor points. Similar to Ren et al. [2015], RPN was trained with each mini-
batch containing a single planning problem, with a 1:2 ratio of positive and
negative samples. The negative samples are randomly sampled from the set of
negative anchor points.

We trained the network by minimizing the cross-entropy loss for the select few
anchor points using the Adam optimizer [Kingma and Ba, 2015] with
$\beta_{1}=0.9$, $\beta_{2}=0.98$, and $\epsilon=1e^{-9}$. We varied the
learning rate as proposed in Vaswani et al. [2017] with warm-up steps of 3200.
Each model was trained for 70k epochs with a batch size of 400. The models
were trained on one machine with 4 NVIDIA 2080GTX graphics card and took 30hrs
to train.

####  4.1.5 Metrics

To evaluate our planning algorithm, we tested it across different environments
for random start and goal pairs. Since all the planners have some form of
optimality criteria, it is unclear when to stop searching. A common practice
in the motion planning community is to set a cost threshold for each path. The
planners stopped searching when the solution path length is below this
threshold or a time limit is reached. The cost threshold was the path length
from an RRT* planner that searched the space till a solution was found and was
collected as part of the validation set. The metrics we are interested in are
the planning accuracy, the time it takes to find a solution, and the number of
vertices on the planning tree. All the statistics we report, apart from
accuracy, consider successful solutions which terminated before the time
limit.

###  4.2 Point Robot Model

Figure 3: Planned paths for two different planning problems. Top Row: Random
Forest environment. Bottom Row: Maze environment. Left Column: The ground
truth path from the validation data. Center Column: The prediction probability
from the MPT. Right Column: The masked map and the planned trajectory using
RRT*.

To evaluate the performance of our method, we compared it with RRT* and IRRT*
without the aid of MPT. The aided versions of these planners are called MPT-
RRT* and MPT-IRRT*. The first set of experiments examined the network’s
ability to generalize to environments of the same size as the training
dataset. Hence we evaluated it across 2500 different environments with a
single start and goal pair for both the Maze and Random Forest environment.
The summary statistics of the experiment are reported in Table 1. The MPT
planners reduce the planning time and vertex count of the planning tree
substantially. In Fig. 3, we show two examples of a planned path from this
experiment.

As the start and goal points are sampled randomly, the statistic’s median
values only provide a partial picture of the benefits of MPT. To better
understand the advantages of MPT, we visualize the distribution of the
planning time and vertices on the planning tree in Fig. 4 for the Random
Forest and Maze environment using Letter-value plots [Hofmann et al., 2011].
These plots help to observe the tail of the distribution of the metrics. Naive
RRT* has a heavier tail distribution than MPT-RRT*, which is as expected
because, for start and goal pairs further away, the planner needs to generate
a denser graph, requiring more time. The MPT-RRT* only searches on the region
proposed by the MPT, and as a result, needs fewer tree vertices to reach the
goal.

We also observe that IRRT* also performs considerably better than RRT* in the
Random Forest environment and achieves similar planning time and planning tree
vertices compared to the aided planners. This is because IRRT*, like MPT’s,
reduces the planning search space once an initial solution is found by
bounding the initial path with an Ellipse. Although this works for the Random
Forest environments, for more complex maps such as the Maze, we see that MPT-
RRT* and MPT-IRRT* perform considerably better than IRRT*.

Figure 4: Planning statistics for the Point Robot Model. Left: The planning time for different SMP and MPT aided SMPs. Right: Number of vertices in the planning tree for SMP and MPT aided SMPs. MPT-RRT* and MPT-IRRT* consistently reduce the planning time and the vertices in the planning tree, resulting in a lower variance of these statistics for these planners. Table 1: Comparing planning accuracy, median planning time, and median number of vertices in the planning tree for the Point Robot Model on environments of the same size as the training data. Environment | Random Forest | Maze
---|---|---
Accuracy | Time (sec) | Vertices | Accuracy | Time (sec) | Vertices
RRT* | 99.88% | 5.44 | 3227.5 | 100% | 5.36 | 2042
IRRT* | 99.88% | 0.42 | 267 | 100% | 3.13 | 1393.5
MPT-RRT* | 97.68% | 0.20 | 251 | 98.96% | 0.83 | 615
MPT-IRRT* | 97.68% | 0.07 | 133 | 98.96% | 0.74 | 557

The next experiment we conducted evaluated the planner’s ability to generalize
to map sizes that were not part of the training data. We tested the model on
four different map sizes, 240$\times$240, 360$\times$240, 480$\times$240, and
560$\times$560, of the Random Forest environment on 1000 randomly generated
maps with an increasing number of obstacles. The same MPT model as before was
used without any retraining or fine-tuning. Three successful planned paths are
shown in Fig. 5. The planning statistics are summarized in Table 2.

We notice that the MPT-RRT* and MPT-IRRT* can achieve lower planning time and
tree vertices count compared to RRT* without any additional training or fine-
tuning. Thus MPT generalizes well to maps of different sizes. Similar to the
previous experiment with fixed map sizes, IRRT* achieves similar performance
to MPT-RRT*, although MPT-IRRT* achieves the best performance of the bunch.
Better performance of MPT-IRRT* could be attributed to the fact that unlike
IRRT*, which has to search the entire space to find an initial solution, MPT-
IRRT* samples only within the proposed region, thus arriving at an initial
solution faster with fewer vertices on the tree.

For the map size of 560x560, we do see a significant drop in accuracy. This
drop in performance might be because the current depth of MSA cannot capture
the global information for tokens that are far away. With larger map sizes,
the number of input tokens increases quadratically, and as noted in Vaswani et
al. [2017], a deeper MSA model performs better for longer input sequences. The
problem requires further investigation and is an interesting future question
to look at.

Figure 5: Plot of paths for Random Forest environments of different size. The architecture of the MPT Model allows flexibility in planning for environments of different sizes. Left: 240$\times$240\. Center: 360$\times$240\. Right: 480$\times$240. Table 2: Comparing planning accuracy, planning time, and number of vertices in the tree for Point Robot on maps of the different sizes of the Random Forest environment. Map Sizes |  | Number of
---
Obstacles
| RRT* | IRRT* | MPT-RRT* | MPT-IRRT*
240$\times$240 | 20 | Accuracy | 100% | 100% | 95.4% | 95.4%
Time (sec) | 4.862 | 0.25 | 0.27 | 0.015
Vertices | 3199 | 229 | 402.5 | 52
360$\times$240 | 35 | Accuracy | 98.9% | 98.9% | 97.4% | 97.4%
Time (sec) | 5.93 | 0.29 | 0.27 | 0.06
Vertices | 3660 | 257 | 377 | 118
480$\times$240 | 50 | Accuracy | 99.4% | 99.4% | 96.3% | 96.3%
Time (sec) | 6.31 | 0.39 | 0.27 | 0.07
Vertices | 3480 | 291 | 348 | 129.5
560$\times$560 | 100 | Accuracy | 100% | 100% | 75.6% | 75.6%
Time (sec) | 6.71 | 0.28 | 0.22 | 0.051
Vertices | 3854 | 203 | 262 | 101
Table 3: Planning accuracy and time for Dubins Car Model for the Random Forest Environment. Planner | SST | MPT-SST
---|---|---
Accuracy | 92.2% | 84.9%
Time (sec) | 28 | 15.8

###  4.3 Dubins Car Model

To examine if the MPT can aid in the planning of non-holonomic systems, we
trained a new MPT model to plan for a Dubins Car Model. We tested the trained
model on 1000 Random Forest environments and compared the metrics with those
from the naive SST planner. The SST aided by the MPT is called MPT-SST. Two
trajectories from this experiment is plotted in Fig. 6 (Left and Center). We
can observe that the MPT highlights regions so that SST can generate a
kinematically feasible path.

Planners such as SST, which sample in the control space, do not scale well to
larger planning spaces, and MPT can help with reducing the search space. The
timing results reported in Table 3 support this claim. MPT-SST can plan nearly
twice as fast as naive SST. The distribution of the planning time for MPT-SST
and SST is given in Fig. 6 (Right).

Figure 6: MPT can also be trained to aid SMP planners for non-holonomic
robots. Left and Center: Planned paths on Random Forest environment using MPT-
SST. Right: The Letter-value plot comparing the planning time between SST and
MPT-SST. MPT is able to reduce the planning time for SST.

Figure 7: Path generated by MPT-RRT* on two random start and goal pairs on a
map of building 079 at University of Freiburg. The model used to generate the
patches were trained only on synthetic environments.

###  4.4 Real World Maps

The final experiment we performed was to test the model trained on Point Robot
Model on a real-world map. We obtained a map of building 079 University of
Freiburg from publicly available 2D lidar scan data. We used the MPT model,
trained on the synthetic Random Forest and Maze environment for the Point
Robot, to plan on 20 randomly sampled start and goal pairs. MPT-RRT* was able
to solve 17 of the 20 planning problems, and we plot two of these trajectories
in Fig. 7. This result is encouraging since it shows MPT’s can solve planning
problems in real-world maps by training on synthetic maps. Robots that have
such generalization capabilities will be able to adapt to dynamic environments
easily.

##  5 Conclusion

In this work, we have shown an application of the transformer model for motion
planning tasks. Unlike prior methods that need to retrain models for maps of
different sizes, we leverage the ability of transformers to handle sequences
of different lengths and parallelize long-term dependencies without recursions
for planning problems with different map sizes. By combing MPT with SMP, we
could generate paths faster, with fewer tree nodes for different environments
and robot models.

While the initial results are impressive, there are loads of places to
improve. One being, the performance of the planner with maps that are larger
than training data size. In NLP and vision tasks, larger transformer models
seem to improve performance. Similarly, larger MPT models may result in better
accuracy for these maps. Extending the planner to higher dimensional robots
such as robotic arms or drones is also an interesting problem with many real-
world applications. Because of the higher dimensions of these spaces, sparse
Transformer models[Roy et al., 2021, Beltagy et al., 2020] would be an ideal
choice for these problems.

## References

  * Arslan et al. [2017] Oktay Arslan, Karl Berntorp, and Panagiotis Tsiotras.  Sampling-based algorithms for optimal motion planning using closed-loop prediction.  In _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 4991–4996, 2017.  doi: 10.1109/ICRA.2017.7989581.
  * Banino et al. [2018] Andrea Banino, Caswell Barry, Benigno Uria, Charles Blundell, Timothy Lillicrap, Piotr Mirowski, Alexander Pritzel, Martin Chadwick, Thomas Degris, Joseph Modayil, Greg Wayne, Hubert Soyer, Fabio Viola, Brian Zhang, Ross Goroshin, Neil Rabinowitz, Razvan Pascanu, Charlie Beattie, Stig Petersen, and Dharshan Kumaran.  Vector-based navigation using grid-like representations in artificial agents.  _Nature_ , 557, 05 2018.  doi: 10.1038/s41586-018-0102-6.
  * Beltagy et al. [2020] Iz Beltagy, Matthew E. Peters, and Arman Cohan.  Longformer: The long-document transformer, 2020.
  * Chaplot et al. [2021] Devendra Singh Chaplot, Deepak Pathak, and Jitendra Malik.  Differentiable spatial planning using transformers, 2021.  URL https://openreview.net/forum?id=n4IMHNb8_f.
  * Chen et al. [2016] Tianqi Chen, Bing Xu, Chiyuan Zhang, and Carlos Guestrin.  Training deep nets with sublinear memory cost, 2016.
  * Dosovitskiy et al. [2021] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=YicbFdNTTy.
  * Emmons et al. [2020] Scott Emmons, Ajay Jain, Misha Laskin, Thanard Kurutach, Pieter Abbeel, and Deepak Pathak.  Sparse graphical memory for robust planning.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 5251–5262. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/385822e359afa26d52b5b286226f2cea-Paper.pdf.
  * Gammell et al. [2014] Jonathan D. Gammell, Siddhartha S. Srinivasa, and Timothy D. Barfoot.  Informed rrt*: Optimal sampling-based path planning focused via direct sampling of an admissible ellipsoidal heuristic.  In _2014 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pages 2997–3004, 2014.  doi: 10.1109/IROS.2014.6942976.
  * Gammell et al. [2015] Jonathan D Gammell, Siddhartha S Srinivasa, and Timothy D Barfoot.  Batch informed trees (bit*): Sampling-based optimal planning via the heuristically guided search of implicit random geometric graphs.  In _2015 IEEE international conference on robotics and automation (ICRA)_ , pages 3067–3074. IEEE, 2015.
  * Hart et al. [1968] Peter E. Hart, Nils J. Nilsson, and Bertram Raphael.  A formal basis for the heuristic determination of minimum cost paths.  _IEEE Transactions on Systems Science and Cybernetics_ , 4(2):100–107, 1968.  doi: 10.1109/TSSC.1968.300136.
  * Hofmann et al. [2011] Heike Hofmann, Karen Kafadar, and Hadley Wickham.  Letter-value plots: Boxplots for large data.  Technical report, had.co.nz, 2011.
  * Ichter et al. [2018] Brian Ichter, James Harrison, and Marco Pavone.  Learning sampling distributions for robot motion planning.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 7087–7094. IEEE, 2018.
  * Johnson et al. [2020] Jacob J. Johnson, Linjun Li, Fei Liu, Ahmed H. Qureshi, and Michael C. Yip.  Dynamically constrained motion planning networks for non-holonomic robots.  In _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pages 6937–6943, 2020.  doi: 10.1109/IROS45743.2020.9341283.
  * Karaman and Frazzoli [2011] Sertac Karaman and Emilio Frazzoli.  Sampling-based algorithms for optimal motion planning.  _The international journal of robotics research_ , 30(7):846–894, 2011.
  * Kavraki et al. [1996] L.E. Kavraki, P. Svestka, J.-C. Latombe, and M.H. Overmars.  Probabilistic roadmaps for path planning in high-dimensional configuration spaces.  _IEEE Transactions on Robotics and Automation_ , 12(4):566–580, 1996.
  * Kingma and Ba [2015] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.  URL http://arxiv.org/abs/1412.6980.
  * LaValle and James J. Kuffner [2001] Steven M. LaValle and Jr. James J. Kuffner.  Randomized kinodynamic planning.  _The International Journal of Robotics Research_ , 20(5):378–400, 2001.  doi: 10.1177/02783640122067453.  URL https://doi.org/10.1177/02783640122067453.
  * Li et al. [2021] Linjun Li, Yinglong Miao, Ahmed H Qureshi, and Michael C Yip.  Mpc-mpnet: Model-predictive motion planning networks for fast, near-optimal planning under kinodynamic constraints.  _IEEE Robotics and Automation Letters_ , 6(3):4496–4503, 2021.
  * Li et al. [2016] Yanbo Li, Zakary Littlefield, and Kostas E. Bekris.  Asymptotically optimal sampling-based kinodynamic planning.  _The International Journal of Robotics Research_ , 35(5):528–564, 2016.  doi: 10.1177/0278364915614386.  URL https://doi.org/10.1177/0278364915614386.
  * Liu et al. [2021] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows, 2021.
  * Qureshi et al. [2020] A. H. Qureshi, J. Dong, A. Choe, and M. C. Yip.  Neural manipulation planning on constraint manifolds.  _IEEE Robotics and Automation Letters_ , 5(4):6089–6096, 2020.
  * Qureshi et al. [2019] Ahmed H. Qureshi, Anthony Simeonov, Mayur J. Bency, and Michael C. Yip.  Motion planning networks.  In _2019 International Conference on Robotics and Automation (ICRA)_ , pages 2118–2124, 2019.  doi: 10.1109/ICRA.2019.8793889.
  * Qureshi and Ayaz [2015] Ahmed Hussain Qureshi and Yasar Ayaz.  Intelligent bidirectional rapidly-exploring random trees for optimal motion planning in complex cluttered environments.  _Robotics and Autonomous Systems_ , 68:1–11, 2015.
  * Qureshi and Ayaz [2016] Ahmed Hussain Qureshi and Yasar Ayaz.  Potential functions based sampling heuristic for optimal path planning.  _Autonomous Robots_ , 40(6):1079–1093, 2016.
  * Qureshi et al. [2020] Ahmed Hussain Qureshi, Yinglong Miao, Anthony Simeonov, and Michael C Yip.  Motion planning networks: Bridging the gap between learning-based and classical motion planners.  _IEEE Transactions on Robotics_ , 2020.
  * Ren et al. [2015] Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun.  Faster r-cnn: Towards real-time object detection with region proposal networks.  In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.  URL https://proceedings.neurips.cc/paper/2015/file/14bfa6bb14875e45bba028a21ed38046-Paper.pdf.
  * Roy et al. [2021] Aurko Roy, Mohammad Saffar, Ashish Vaswani, and David Grangier.  Efficient content-based sparse attention with routing transformers.  _Transactions of the Association for Computational Linguistics_ , 9:53–68, 02 2021.  doi: 10.1162/tacl_a_00353.
  * Srinivas et al. [2018] Aravind Srinivas, Allan Jabri, Pieter Abbeel, Sergey Levine, and Chelsea Finn.  Universal planning networks: Learning generalizable representations for visuomotor control.  In _International Conference on Machine Learning_ , pages 4732–4741. PMLR, 2018.
  * Şucan et al. [2012] Ioan A. Şucan, Mark Moll, and Lydia E. Kavraki.  The Open Motion Planning Library.  _IEEE Robotics & Automation Magazine_, 19(4):72–82, December 2012.  doi: 10.1109/MRA.2012.2205651.  https://ompl.kavrakilab.org.
  * Tahir et al. [2018] Zaid Tahir, Ahmed H Qureshi, Yasar Ayaz, and Raheel Nawaz.  Potentially guided bidirectionalized rrt* for fast optimal path planning in cluttered environments.  _Robotics and Autonomous Systems_ , 108:13–27, 2018.
  * Tamar et al. [2016] Aviv Tamar, Yi Wu, Garrett Thomas, Sergey Levine, and Pieter Abbeel.  Value iteration networks.  In _Advances in Neural Information Processing Systems_ , pages 2154–2162, 2016.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.
  * Williams [1992] Ronald J Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine learning_ , 8(3-4):229–256, 1992.
  * Zucker et al. [2008] Matt Zucker, James Kuffner, and J Andrew Bagnell.  Adaptive workspace biasing for sampling-based planners.  In _2008 IEEE International Conference on Robotics and Automation_ , pages 3757–3762. IEEE, 2008.

## Appendix

##  Appendix A Network Architecture

Table 4: Network architecture of the Feature Extractor Layer | Dimension
---|---
2D Convolution | [2, 6, 5, 0]
2D Maxpool | [2]
ReLU |
2D Convolution | [6, 16, 5, 0]
2D Maxpool | [2]
ReLU |
Convolution | [16, 512, 5, 5]

In this section, we detail the network architecture used in our experiments.
The Transformer architecture is similar to the ones proposed in Vaswani et al.
[2017]. We used 6 layers of encoder block, each consisting of 3 heads. The
dimension of the keys and queries was set at 512, and the dimension of the
value was set at 256. The architecture of our feature extractor is given in
Table 4. For the convolution layer, the dimensions in the brackets represent
[Input Channel Size, Output Channel Size, Kernel Size, Stride], and for the
Maxpool layer, it represents the Kernel Size.

##  Appendix B Successful Trajectories

In this section we plot few more successful trajectories of fixed map sizes
for the Random Forest and Maze environments in Fig. 8 and 9 respectively.

Figure 8: Three different trajectories planned successfully using MPT on the
Random Forest environment. Left Column: The ground truth path in the
validation dataset. Middle Column: The probability estimate made by the MPT.
Right Column: The planned path using RRT* on the region proposed by the MPT.

Figure 9: Three different trajectories planned successfully using MPT on the
Maze environment. Left Column: The ground truth path in the validation
dataset. Middle Column: The probability estimate made by the MPT. Right
Column: The planned path using RRT* on the region proposed by the MPT.

##  Appendix C Distribution of Metrics for Maps of Different Sizes

The distribution of metrics reported in Table 2 is given in Fig. 10.

Figure 10: The distribution of metrics for maps of different sizes. Top: The
distribution of planning time for different algorithms for maps of different
sizes. Bottom: The distribution of the number of vertices in the planning tree
for different algorithms for maps of different sizes. MPT aided planners can
achieve faster planning times and lower variance in the metrics across maps of
different sizes, which were not part of the training data.

##  Appendix D Failed Trajectories

In Fig. 11 and 12, we visualize some of the trajectories that MPT could not
solve for the Random Forest and Maze environment, respectively. These failures
were mainly caused by a single patch that was not classified as part of the
path, resulting in a map with no valid paths.

Figure 11: Three different trajectories planned unsuccessfully using MPT on
the Random Forest environment. Left Column: The ground truth path in the
validation dataset. Middle Column: The probability estimate made by the MPT.
Right Column: The planned path using RRT* on the region proposed by the MPT.
We observe that in all three cases the failure to predict a few patches
resulted in a scenario with no valid paths.

Figure 12: Three different trajectories planned unsuccessfully using MPT on
the Maze environment. Left Column: The ground truth path in the validation
dataset. Middle Column: The probability estimate made by the MPT. Right
Column: The planned path using RRT* on the region proposed by the MPT. We
observe that in all three cases the failure to predict a few patches resulted
in a scenario with no valid paths.


# Scintillation Timing Characteristics of Common Plastics for Radiation
Detection Excited With 120 GeV Protons

Burak BILKI1,2, Nilay BOSTAN2, Ohannes Kamer KÖSEYAN2, Emrah TIRAŞ2,3, James
WETZEL2,4
1Department of Mathematics<EMAIL_ADDRESS>Beykent University  Istanbul
Turkey
2Department of Physics and Astronomy  The University of Iowa  Iowa City  Iowa
USA
3Department of Physics and Astronomy  Iowa State University  Ames  Iowa  USA
4Department of Physics  Augustana College  Rock Island  Illinois  USA

###### Abstract

The timing characteristics of scintillators must be understood in order to
determine which applications they are appropriate for. Polyethylene
naphthalate (PEN) and polyethylene teraphthalate (PET) are common plastics
with uncommon scintillation properties. Here, we report the timing
characteristics of PEN and PET, determined by exciting them with 120 GeV
protons. The test beam was provided by Fermi National Accelerator Laboratory,
and the scintillators were tested at the Fermilab Test Beam Facility. PEN and
PET are found to have dominant decay constants of 34.91 ns and 6.78 ns,
respectively.

###### keywords:

Scintillators, scintillation and light emission processes (solid, gas and
liquid scintillators), Scintillators and scintillating fibres and light guides

##  1 Introduction

After observing that common plastics polyethylene naphthalate (PEN) and
polyethylene teraphthalate (PET) have very good scintillation properties [1,
2], and reasonable tolerance to gamma radiation exposure [3], interest in
using PEN and PET in scintillation applications has grown. Of note is their
potential application in collider physics experiments, where huge amounts of
active media are required in very large detectors. These plastics’ low cost
and ease of manufacture are attractive properties. However, for these
applications, scintillators need to be relatively fast and bright.

In order to study their properties for future applications, similar samples of
PEN and PET were prepared for testing at Fermilab Test Beam Facility (FTBF) by
exciting them with 120 GeV protons. The proton beam at FTBF produces $5\times
10^{5}$ protons per spill with a spot size of 6 mm. The experiment was
designed to measure the decay constants of PEN and PET. The decay constant is
the time it takes for the scintillator to diminish to 1/e of its maximum
signal amplitude. As a particle traverses a scintillator, the particle excites
the molecules within the material, which then radiate photons during de-
excitation.

Timing is important for experiments like those at CERN, where the LHC [4] has
a design collision frequency of 40 MHz, or one collision every 25 ns.
Considering the potential for unmanageable pileup, or the excessive number of
proton-proton collisions per beam bunch crossing, a scintillator must dispatch
its entire signal ideally within, in the case of LHC experiments, 25 ns.

In addition, PEN was recently identified as a relatively efficient wavelength-
shifter for vacuum ultraviolet (VUV) photons and found an immediate
implementation at the CERN ProtoDUNE-DP [5] experiment. 30 out of 36
photomultiplier tubes (PMTs) of ProtoDUNE-DP have PEN sheets installed on top
of their windows, and the PMTs have been operating in liquid argon
successfully since Summer 2019 [6]. Studying the scintillation properties of
such common plastics under various test conditions will have critical
importance for a wide range of implementations in the future.

##  2 Experimental Setup

Three samples each of PEN and PET were prepared for testing. The PEN sample
was cut to a square 100 mm x 100 mm x 1 mm from an 8.5” x 11” sheet of Teijin
Plastics Scinterex using an end mill. The edges were polished with a polishing
wheel. The PET sample was delivered in 100 mm x 100 mm x 1 mm squares from
Goodfellow.com, and the edges were polished with a polishing wheel.

Each tile had a ‘keyhole’ groove machined into its surface around the
perimeter to hold an optical fiber. A ‘Y11’ Kuraray wavelength shifting (WLS)
fiber was inserted into the groove, Fig 1. An injection molded connector was
glued to the end of the fiber, and the fiber edge was polished. Each tile was
placed into a 3D printed housing, Fig 2, and the WLS fiber was coupled to a
clear fiber. The clear fiber was connected to a clear acrylic ‘cookie’ for
coupling to an R7600 PMT for readout.

Figure 1: Scintillator tile with keyhole groove and Y11 fiber.

(a)

(b)

Figure 2: 3D printed frame for holding the scintillators to be tested (a) and
experimental setup in the beam area showing the sample holder, the clear
fiber, and the PMT (b).

The signals recorded from the R7600 PMT were fed into a Tektronix TDS 5034
digital oscilloscope, and events were recorded in a Tektronix proprietary
‘.wfm’ format , using the ‘Fast Frame’ mode of the oscilloscope. This allowed
the recording of twenty-five thousand events to be saved in a single file,
which we could later analyze. The events were triggered in the central 1 cm x
1 cm area of the tiles by the coincidence of two external scintillator paddle
signals.

##  3 Experimental Results

The twenty-five thousand events were extracted from the .wfm format and
ntupled using the ROOT framework [7], collating all events into a single ROOT
file. The pedestal was subtracted from each event, and the amplitudes of the
event waveforms were averaged using ROOT’s TProfile function in order to
create an average waveform profile. The TProfile function averages the
contents of each bin over the twenty-five thousand events, thereby providing a
very nice picture of the scintillator’s photoluminescence and decay profile.

Generally, scintillators have a _fast_ and _slow_ component to their time-
dependent luminescence profile, corresponding respectively to its
recombination centers and electron traps.

The waveform profile of a scintillator can therefore be modeled as a
combination of slow and fast components [8]:

$\frac{A}{\tau_{fast}}e^{-t/\tau_{fast}}+\frac{B}{\tau_{slow}}e^{-t/\tau_{slow}}$
(3.1)

The waveform profile of PET was fit with equation 3.1, though the fast
component of PET was an order of magnitude greater in intensity than the slow
component.

PEN was found to be well described by a single exponential function of the
form:

$\frac{A}{\tau_{fast}}e^{-t/\tau_{fast}}$ (3.2)

implying dominance of a single photoluminescence mechanism in PEN.

Figure 3 shows the waveform profiles and fits of PEN (a) and PET (b). PEN and
PET were found to have dominant decay constants of 34.91 ns and 6.78 ns,
respectively.

(a)  (b)

Figure 3: Average waveform profiles with decay constants of PEN (a) and PET
(b) using 120 GeV protons.

##  4 Conclusions

Waveform profiles of scintillating plastic sheets polyethylene naphthalate
(PEN) and polyethylene teraphthalate (PET) were recorded by placing the
scintillators in and perpendicular to an accelerated beam of 120 GeV protons.
The waveform profiles were measured using a PMT read out with an oscilloscope.
The exponential fits to the average waveforms indicate that PEN and PET have
dominant decay constants of 34.91 ns and 6.78 ns, respectively.

##  5 Acknowledgements

The authors would like to thank the staff at Fermi National Accelerator
Laboratory, in particular Mandy Rominsky, manager of the Fermilab Test Beam
Facility.

## References

  * [1] Nakamura, H., Shirakawa, Y., Takahashi, S., & Shimizu, H. (2011). Evidence of deep-blue photon emission at high efficiency by common plastic. EPL (Europhysics Letters), 95(2), 22001.
  * [2] Nakamura, H., Shirakawa, Y., Kitamura, H., Yamada, T., Shidara, Z., Yokozuka, T., … & Takahashi, S. (2013). Blended polyethylene terephthalate and polyethylene naphthalate polymers for scintillation base substrates. Radiation Measurements, 59, 172-175.
  * [3] Wetzel, J., Tiras, E., Bilki, B., Onel, Y., & Winn, D. (2016). Radiation damage and recovery properties of common plastics PEN (Polyethylene Naphthalate) and PET (Polyethylene Terephthalate) using a 137Cs gamma ray source up to 1.4 Mrad and 14 Mrad. Journal of Instrumentation, 11(08), P08023.
  * [4] Evans, L., & Bryant, P. (2008). LHC machine. Journal of instrumentation, 3(08), S08001.
  * [5] Scarpelli, A. (2018) ProtoDUNE and a Dual-Phase LArTPC, Neutrino Oscillation Workshop (NOW2018), 9 - 16 September, 2018, Rosa Marina (Ostuni, Brindisi, Italy), https://arxiv.org/pdf/1902.04780.pdf
  * [6] Cuesta, C. (2019) Status of ProtoDUNE Dual Phase, European Physical Society Conference on High Energy Physics - EPS-HEP2019 -10-17 July, 2019, Ghent, Belgium, https://arxiv.org/pdf/1910.10115.pdf
  * [7] Antcheva, I., Ballintijn, M., Bellenot, B., Biskup, M., Brun, et al. (2009). ROOT-A C++ framework for petabyte data storage, statistical analysis and visualization. Computer Physics Communications, 180(12), 2499-2512.
  * [8] Swiderski, L., Moszynski, M., Syntfeld-Kazuch, A., Szawlowski, M., & Szczesniak, T. (2014). Measuring the scintillation decay time for different energy depositions in NaI: Tl, LSO: Ce and CeBr3 scintillators. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 749, 68-73.