# Deep learning in physics: a study of dielectric quasi-cubic particles in a
uniform electric field

Zhe Wang
Energy Research Institute
Nanyang Technological University
637141 Singapore
<EMAIL_ADDRESS>
&Claude Guet
Energy Research Institute
Nanyang Technological University
637141 Singapore
School of Materials Science and Engineering
Nanyang Technological University
639798 Singapore
<EMAIL_ADDRESS>

###### Abstract

Solving physics problems for which we know the equations, boundary conditions
and symmetries can be done by deep learning. The constraints can be either
imposed as terms in a loss function or used to formulate a neural ansatz. In
the present case study, we calculate the induced field inside and outside a
dielectric cube placed in a uniform electric field, wherein the dielectric
mismatch at edges and corners of the cube makes accurate calculations
numerically challenging. The electric potential is expressed as an ansatz
incorporating neural networks with known leading order behaviors and
symmetries and the Laplace’s equation is then solved with boundary conditions
at the dielectric interface by minimizing a loss function. The loss function
ensures that both Laplace’s equation and boundary conditions are satisfied
everywhere inside a large solution domain. We study how the electric potential
inside and outside a quasi-cubic particle evolves through a sequence of shapes
from a sphere to a cube. The neural network being differentiable, it is
straightforward to calculate the electric field over the whole domain, the
induced surface charge distribution and the polarizability. The neural network
being retentive, one can efficiently follow how the field changes upon
particle’s shape or dielectric constant by iterating from any previously
converged solution. The present work’s objective is two-fold, first to show
how an a priori knowledge can be incorporated into neural networks to achieve
efficient learning and second to apply the method and study how the induced
field and polarizability change when a dielectric particle progressively
changes its shape from a sphere to a cube.

##  1 Introduction

Solving physics problems for which we know the equations, boundary conditions
and symmetries by deep learning follows from the universal approximation
theorem [1, 2, 3], which states a sufficiently deep artificial neural network
(ANN) can approximate any well-behaved function with a finite number of
parameters. In 1994, Meade and Fernandez approximated the solution of linear
[4] and nonlinear [5] ordinary differential equations using a single-layer
perceptron. This approach was soon generalized by Lagaris et. al. [6] and
applied to two-dimensional Poisson equations with various source terms in a
rectangular domain. Their trial functions consisted of two terms: the first
one which contained no trainable variable satisfied the boundary conditions,
whereas the second one involved neural networks, whose parameters, e.g.
weights and biases, were adjusted to minimize a loss function. A neural
network solution of Poisson’s equation in a 3D domain with irregular
boundaries was achieved by including the boundary conditions into the loss
function [7]. Alternatively, McFall and Mahan constructed a proper trial
function so that the boundary conditions were automatically satisfied [8].
Massive growth of available scientific data has introduced a new flavor into
the ANN approach to differential equations. By incorporating data and
governing equations into the loss functions, the physics-informed neural
networks enable inferring hidden physics from measured data [9, 10, 11] with
successful applications on the visualization of turbulent flows [12, 13, 14],
where multiple scales interact, and the design of metamaterials in nano-optics
[15], where the finite size effect dominates. However, as pointed out by Wong
et al. [16], high-dimensional, non-convex loss functions require significant
optimization efforts, calling for sophisticated hyper-parameter optimization
[17, 18, 19] and transfer neuroevolution [16] algorithms.

In this work, we combine physics knowledge with recent advances in deep
learning to calculate the induced electric field of semiconductor colloidal
nanocrystals in an external photon field. Semiconductor colloidal nanocrystals
show numerous advantages for powerful opto-electronic materials, as their
optical properties change with shape and size due to quantum effects [20, 21].
In order to characterize correctly their absorption/emission properties, one
needs to know the electric field induced by the external photon field [22,
23]. As the particle’s size is much smaller than the photon wavelength, the
inner field results from a homogeneous applied electric field of amplitude
given by the laser intensity. Whereas there are analytical solutions of
Laplace’s equation for dielectric particles with shapes that allow to define a
system of curvilinear coordinates such as sphere, ellipsoids, torus, etc. [20,
24, 25, 26], numerical solvers are required for other shapes. The case of a
cube is challenging because the edges and corners lead to sharp variations of
the induced fields. To our knowledge there is only one published calculation
based on finite element methods applied to a cube with relative dielectric
constant ranging from to $2$ to $10$ [22]. The authors claimed a $1\%$
accuracy at the far field of their finite element simulation domain. At low
values of the relative dielectric constant considered, they found that the
electric field at the center of the particle is lower for a cube than for a
sphere. This would imply a lower polarizability for a cube than for a sphere,
at variance with accurate calculations [27, 28, 29, 30, 31, 32, 33]. Note that
these herein cited accurate calculations all focused on solving a surface
integral equation to estimate the dielectric polarizability and consequently
did not provide the induced electric field inside and outside the cube.

In this paper, we introduce an alternative method to calculate the electric
field inside and outside a dielectric nanoparticle embedded in another
dielectric medium. From the field inside, the polarizability is extracted.
More specifically, we approximate the solution of Laplace’s equation by a
function combining a known analytical term and an ANN. Then, instead of
solving numerically Laplace’s equation with boundary conditions, we express
the problem as an optimization problem and construct a loss function which can
be minimized by machine learning algorithms to yield the full electric
potential inside and outside the particle. One clear advantage of the ANN
approach over finite element methods is that it provides the solution of
Laplace’s equation as a differentiable function that can be used as such.
Additionally, the retentive nature of neural networks allows for a systematic
tracking of the evolution of induced fields as dielectric particles deviate
from canonical spherical and cubic shapes. Note that physical nanoparticles
usually have rounded edges and corners, see e.g. Fig. 2s in [34] and in [35]
for perovskite nanocrystals.

Most previous works considered partial differential equations in a homogeneous
domain with explicit boundary conditions. The electrostatic problem of a
colloidal particle requires to treat accurately the discontinuity of the
displacement field at the interface. To the best of our knowledge, the present
work is the first ANN-based attempt to address this problem, namely solving
Laplace’s equation in a three-dimensional piece-wise homogeneous domain with
implicit boundary conditions on an irregular interface. In addition to
formulating the best loss functions, one can substantially reduce the
optimization effort by incorporating known symmetries and specific features,
e.g. leading order behaviors, into the ANN ansatz.

The rest of the paper is organized as follows: in Sec. 2, we present the
physics model. Starting from a general solution in the form of a linear
combination of spherical harmonics, we replace the higher order terms by a
function containing neural networks and construct a loss function that
includes all constraints. The architecture of the neural network is then
discussed. In Sec. 3.1 , we benchmark the proposed ANN method by studying
dielectric spheroids for which analytical solutions are known [24]. In Sec.
3.2, we study the evolution of polarizabilities and induced fields inside and
outside a dielectric particle through a sequence of shapes from a sphere to a
cube. Finally, conclusions are drawn in Sec. 4, with a highlight on future
works.

##  2  Neural Network solution of Laplace’s equation

###  2.1 Governing equations

Consider a neutral homogeneous dielectric 3D particle, with surface $S$ and
permittivity $\epsilon_{1}$, embedded in a homogeneous medium with
permittivity $\epsilon_{0}$. In a spherical coordinates $(r,\theta,\varphi)$
system whose origin coincides with the center of the particle, we take the
direction of an uniform electric field $\boldsymbol{E}_{\text{ext}}$ to be the
axis from which the polar angle $\theta$ is measured. In the absence of
external charges, the electric potential obeys Laplace’s equation inside and
outside the particle,

$\displaystyle\nabla^{2}\phi=0.$ (1)

The continuity of the tangential components of the electric field
($\boldsymbol{E}=-\nabla\phi$) and the normal component of the displacement
field ($\boldsymbol{D}=\epsilon\boldsymbol{E}$) at the interface $S$ leads to
the following boundary conditions for the potential:

$\displaystyle\phi_{0}$ $\displaystyle=\phi_{1},$ (2a)
$\displaystyle\nabla\phi_{0}\cdot\hat{\boldsymbol{n}}$
$\displaystyle=\epsilon_{r}\nabla\phi_{1}\cdot\hat{\boldsymbol{n}}.$ (2b)

Here, $\hat{\boldsymbol{n}}=\nabla S/|\nabla S|$, is the unit normal vector to
$S$, $\epsilon_{r}=\epsilon_{1}/\epsilon_{0}$ is the relative dielectric
constant, and subscripts 0 and 1 denote quantities outside and inside the
particle, respectively. Furthermore, the electric field tends asymptotically
to the applied field and the potential is zero at the origin, yielding

$\displaystyle\phi_{0}$
$\displaystyle=-E_{\text{ext}}r\cos(\theta),\quad\mbox{as}\quad r\to\infty,$
(3a) $\displaystyle\phi_{1}$ $\displaystyle=0,\quad\mbox{at}\quad r=0,$ (3b)

where $E_{\text{ext}}=|\boldsymbol{E}_{\text{ext}}|$.

Figure 1:  $(a)$ Sequences of shapes from a sphere to a quasi-cube with
increasing values of $N$, cf. Eq. (4) with $a=b=c$. $(b)$ Plot of normalized
volume and surface area of quasi-cubes by the same of the sphere as a function
of $N$.

Within the perspective of approaching a real cube and benchmarking the present
numerical method with analytical solutions, we consider a super-ellipsoidal
inclusion whose surface is described by the equation

$\displaystyle
S(x,y,z)\equiv\left|\frac{x}{a}\right|^{2N}+\left|\frac{y}{b}\right|^{2N}+\left|\frac{z}{c}\right|^{2N}=R^{2N},$
(4)

in the Cartesian coordinates. Here, $a,b,c,R\in\mathbb{R}^{+}$ and the
exponent $N\geq 1$. For $a=b=c$ and $N=1$, Eq. (4) defines a sphere with
radius $R$. A continuous deformation of a sphere with radius $R$ to a cube of
edge length $2R$ as $N$ increases towards $N\in\infty$, is shown in Fig. 1.
For cases where $a$, $b$, $c$ are not equal to each other, one obtains an
ellipsoid for $N=1$ and a rectangular parallelopiped as $N\to\infty$. Note by
passing that “real” nanoparticles have rounded corners, resembling quasi-cubes
with $N\in[1,8]$ in Fig. 1. For instance, by fitting $2D$ images of $100$
perovskite nanocrystals to a superellipse, cf. Eq. (4) with $z=0$, Tang et al.
[35] found that $N\approx 2.65$ for freshly synthesized and $N\approx 1.8$ for
aged samples, respectively.

A general solution to Eq. (1) can be expressed as a linear combination of
spherical harmonics $Y_{l}^{m}(\theta,\varphi)$ weighted by appropriate
scaling factors $r^{l}$ inside and $r^{-(l+1)}$ outside of the dielectric
particle

$\phi=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left[A_{l}^{m}r^{-(l+1)}+B_{l}^{m}r^{l}\right]Y_{l}^{m}(\theta,\varphi),$
(5)

where coefficients $A_{l}^{m}$ and $B_{l}^{m}$ are determined by the boundary
conditions. The homogeneity and the direction of the external electric field
$\boldsymbol{E}_{\text{ext}}$ along the $z$-axis impose an anti-symmetry on
the electric potential with respect to the mid-plane $\theta=\pi/2$, which is
orthogonal to $\boldsymbol{E}_{\text{ext}}$ and crosses the center of the
particle. Thus,

$\displaystyle\phi(r,\theta,\varphi)=-\phi(r,\pi-\theta,\varphi),$ (6)

such that only odd terms remain in the spherical harmonics expansions (5). For
sake of simplicity, we select one of the principal axis of the dielectric
particle to be aligned with $\boldsymbol{E}_{\text{ext}}$. As such, there is a
further mirror symmetry,

$\displaystyle\phi(r,\theta,\varphi)$ $\displaystyle=\phi(r,\theta,-\varphi),$
(7a) $\displaystyle\phi(r,\theta,\varphi)$
$\displaystyle=\phi(r,\theta,\pi+\varphi),$ (7b)
$\displaystyle\phi(r,\theta,\varphi)$
$\displaystyle=\phi(r,\theta,\pi-\varphi).$ (7c)

For $N=1$, only the dipole moment contributes to $\phi_{1}$, leading to a
constant electric field inside the particle [24]. Deviating from a sphere,
$N\geq 2$, however, brings in contributions from higher-order moments with no
clear optimal order of truncation, rendering an analytical solution
unfeasible.

Figure 2:  $(a)$ Sketch for a dielectric particle immersed in a uniform
electric field along the $z$-axis, $\boldsymbol{E}_{\text{ext}}$. $(b)$ A
superposition of collocation points sampled from a uniform distribution and a
Gaussian distribution with mean centered at the dielectric interface $r_{S}$.
$(c)$ Proposed ANN ansatz for electric potentials, see Sec. 2.2. $(d)$ Loss
function which measures the deviation from the Laplace equation (1) and the
boundary conditions (2) on the dielectric interface, see Sec. 2.3. $(e)$
Visualization of radial, polar, and azimuthal components of $\nabla^{2}\phi$
and its deviation from Eq. (1). $(f)$ Recovery of electric fields
$\boldsymbol{E}$ and charge density $\sigma$ from $\phi$. Here and thereafter,
results are plotted in a circular domain of radius $3$, wherein the blue curve
indicates the interface.

###  2.2 ANN ansatz

Our ansatz for solving Laplace’s equation with the boundary conditions given
by Eqs. (2,3) is to explicitly retain the dipole contributions and group the
infinite series of higher order terms into two ANN functions:

$\displaystyle\phi_{0}$
$\displaystyle=-\left[E_{\text{ext}}r-d_{0}r^{-2}\right]\cos(\theta)+\frac{r-R_{\text{max}}}{r_{S}-R_{\text{max}}}H_{0},$
(8a) $\displaystyle\phi_{1}$
$\displaystyle=d_{1}r\cos(\theta)+\frac{r-R_{\text{min}}}{r_{S}-R_{\text{min}}}H_{1},$
(8b)

where the dipole coefficients $d_{0}$ and $d_{1}$ are unknown. Here,
$r_{S}=r_{S}(\theta,\varphi)$ is an implicit solution to Eq. (4) measuring the
distance from the origin to points on the interface, $R_{\text{min}}$ and
$R_{\text{max}}$ are the inner and outer boundaries of the solution domain and
finally, $H_{0}$ and $H_{1}$ denote ANN based functions modeling a deviation
from the leading order dipolar behavior of the induced electric field. A
breakdown of the proposed ansatz (8) is detailed below.

To get around the singularity at the origin of the Laplacian operator in
spherical coordinates [36], we assume a practical origin of the radial
coordinate: $R_{\text{min}}\ll 1$. Similarly, to enable a numerical
calculation, the infinity $r\to\infty$ is replaced by a practical infinity:
$R_{\text{max}}\gg 1$. To mitigate the effect of such a truncation, we modify
the boundary conditions (3) by including the dipole contributions at the
practical infinity and practical origin, leading to the first terms on the
right hand side of the ANN ansatz (8). Then, the approximate solution domain
is a spherical shell with inner and outer radii $R_{\text{min}}$ and
$R_{\text{max}}$, wherein the modified boundary conditions are satisfied by
construction. In the limit, $R_{\text{min}}=0$ and $R_{\text{max}}\to\infty$,
boundary conditions (3) are recovered.

To inform ANN ansatz (8) with the leading octupolar radial trends, namely
$r^{+3}$ inside and $r^{-4}$ outside the particle, as well as the symmetry
constraints of Eqs. (6,7), we select

$\displaystyle H_{0}$
$\displaystyle=\cos(\theta)r^{-4}\text{atanh}\left[\eta\text{NN}_{0}(r^{*},\theta^{*},\varphi^{*};\boldsymbol{\xi}_{0})\right],$
(9a) $\displaystyle H_{1}$
$\displaystyle=\cos(\theta)r^{+3}\text{atanh}\left[\eta\text{NN}_{1}(r^{*},\theta^{*},\varphi^{*};\boldsymbol{\xi}_{1})\right],$
(9b)

so that the parameters $\boldsymbol{\xi}_{i}$ of neural networks
$\text{NN}_{i}$, with $i=0,1$, are adapted to learn deviations from octupolar
radial trends. The symmetry constraints are imposed on $\text{NN}_{i}$ through
a reparameterization of spatial variables $(r^{*},\theta^{*},\varphi^{*})$ and
the inclusion of $\cos(\theta)$. The neural network architectures, the
constant $\eta$, and the reparameterized coordinate variables
$(r^{*},\theta^{*},\varphi^{*})$, are discussed below.

We describe both $\text{NN}_{i}$ as multilayer perceptrons [37] each
consisting of four densely connected hidden layers with $16$ neurons per layer

$\displaystyle\boldsymbol{x}^{[k]}=\tanh(\boldsymbol{W}^{[k]}\cdot\boldsymbol{x}^{[k-1]}+\boldsymbol{b}^{[k]}),$
(10)

where $k$ is the index of the current layer, and the weight
$\boldsymbol{W}^{[k]}$ and bias $\boldsymbol{b}^{[k]}$ operate a linear
transformation of the input vector $\boldsymbol{x}^{[k-1]}$. As the bias of
both input and output layer are selected to be zero vectors, the neural
network eventually consists of $880$ trainable variables. Instead of the
canonical rectified linear unit (ReLU) which vanishes upon second-order
differentiation, we assign the activation function to be the hyperbolic
tangent function. Since the output of $\text{NN}_{i}$ measures a deviation
from the octupolar trend, it must remain bounded, with an amplitude that
depends on the strength of the external field, the geometry of the dielectric
inclusion and the mismatch at the dielectric interface. Therefore, an
additional operation $\text{atanh}(\eta\cdot)$ is included in Eqs. (9) to
transform the output of $\text{NN}_{i}$, ranging from $[-1,1]$, into an
interval $[-\text{atanh}(\eta),\text{atanh}(\eta)]$. A comparison of the
proposed activation function with the canonical linear and $\tanh$ activation
is discussed in Appendix A. In this work, with $E_{\text{ext}}=1$, we select
$\eta=0.99$. The model parameters
$\boldsymbol{\xi}_{i}=[\boldsymbol{W}_{i},\boldsymbol{b}_{i}]$, as well as
$d_{i}$ defined in Eqs. (8a, 8b), are determined by minimizing the loss
function discussed in §2.3.

Noting that the Laplacian operator is of second order, we consider the
following continuous and differentiable transformation of angular variables

$\displaystyle\theta^{*}=-\cos(2\theta)\quad\mbox{and}\quad\varphi^{*}=-\cos(2\varphi),$
(11)

which simultaneously maps $\theta\in[0,\pi]$ and $\varphi\in[-\pi,\pi]$ to the
first quadrant and rescales them to the range $[-1,1]$. Similarly, the radial
variable $r\in[R_{\text{min}},R_{\text{max}}]$ is normalized to the range
$[-1,1]$ using the min-max normalization

$\displaystyle
r^{*}=\left\\{\begin{array}[]{c}-1+2\dfrac{r-\text{min}(r_{S})}{R_{\text{max}}-\text{min}(r_{S})},\quad\mbox{for}\quad
r>r_{S},\\\\[10.0pt]
-1+2\dfrac{r-R_{\text{min}}}{\text{max}(r_{S})-R_{\text{min}}},\quad\mbox{for}\quad
r<r_{S}.\end{array}\right.$ (14)

With $(r^{*},\theta^{*},\varphi^{*})$, the outputs of neural networks are
symmetric with respect to the mid-plane $\theta=0$, and the anti-symmetry is
enforced by the multipliers $\cos(\theta)$ in Eqs. (9). A flowchart for the
proposed ANN ansatz is sketched in Fig. 2$(c)$

###  2.3 Loss function

Given governing equation (1) and boundary conditions (2), we write the loss
function as:

$\displaystyle L=L_{ge}+L_{bc},$ (15)

where $L_{ge}$ and $L_{bc}$ measure the mean squared deviations of the ansatz
functions (8a, 8b) from the exact solutions of the governing equation (1) and
the dielectric interface boundary conditions (2), respectively. Although an
exact minimization of the loss function (15) ensures the uniqueness of the
solution to the boundary value problem, an approximation of that solution by
neural networks yields a small nonzero loss. The first term writes as:

$\displaystyle L_{ge}$
$\displaystyle=\frac{w_{g_{0}}}{N_{0}+N_{b}}\sum_{j=1}^{N_{0}+N_{b}}\left[\tilde{r}_{j}^{\beta_{0}}\sin^{2}(\theta_{j})\nabla^{2}\phi_{0}(\boldsymbol{r}_{j})\right]^{2}+\frac{w_{g_{1}}}{N_{1}+N_{b}}\sum_{j=1}^{N_{1}+N_{b}}\left[\tilde{r}_{j}^{\beta_{1}}\sin^{2}(\theta_{j})\nabla^{2}\phi_{1}(\boldsymbol{r}_{j})\right]^{2},$
(16)

where the $\boldsymbol{r}_{j}=(r_{j},\theta_{j},\varphi_{j})$ denote the
$j$-th collocation points sampled from a superposition of a uniform
distribution and a Gaussian distribution centered at $r_{S}$ as shown in Fig.
2$(b)$; $N_{0}$, $N_{1}$ and $N_{b}$ denote the number of collocation points
outside, inside, and on the interface of the dielectric particle,
respectively. The multiplier $\sin^{2}(\theta_{j})$ is included to compensate
the singularity of the Laplacian operator at $\theta=0,\pi$. The Laplacian has
vanishing magnitude at large $r$ values, whereas it is diverging near
$R_{\text{min}}$. Therefore, the solutions $\phi_{i}$ are not equally
optimized throughout the solution domain. Inspired by van der Meer et al.
[38], we introduce the scaling factors $\tilde{r}_{j}^{\beta_{i}}$, with
$\tilde{r}_{j}=r_{j}/r_{S}(\theta_{j},\varphi_{j})$ defined at each
collocation point. We select the exponents $\beta_{0}=4$ and $\beta_{1}=1$, in
order to ensure that the radial, polar, and azimuthal components of $L_{ge}$
are of the same order throughout. A breakdown of the Laplacian
$\nabla^{2}\phi$ componentwise and its visualization are presented in Fig. 2
$(d)$ and $(e)$. With increasing $N$, larger and larger derivatives associated
with sharper and sharper edges and corners cause a strong mismatch among
components of $\nabla^{2}\phi_{i}$ near the interface. Consequently, large
oscillations of the loss function emerge during the training process, leading
to a slow convergence. In order to balance such a mismatch, weights
$w_{g_{0}}$ and $w_{g_{1}}$ are introduced in Eq. (16) in addition to a
normalization by $r_{S}(\theta_{i},\varphi_{i})$.

The dielectric interface boundary conditions (2) lead to the following loss
function:

$\displaystyle L_{bc}$
$\displaystyle=\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}\left\\{w_{b_{t}}\left[\phi_{0}(\boldsymbol{r}_{j})-\phi_{1}(\boldsymbol{r}_{j})\right]^{2}+\frac{w_{b_{n}}}{|\nabla
S|^{2}}\left[{\epsilon_{r}}^{-1}\nabla\phi_{0}(\boldsymbol{r}_{j})\cdot\nabla
S-\nabla\phi_{1}(\boldsymbol{r}_{j})\cdot\nabla S\right]^{2}\right\\}.$ (17)

The weights $w_{b_{t}}$ and $w_{b_{n}}$ help balance the losses of tangential
and normal boundary conditions during the training. Convergence and accuracy
are significantly improved when the normalization factor $|\nabla
S(r,\theta,\varphi)|$ is included.

##  3 Computational results

Throughout the experiments, the loss function (15) was evaluated over a
sampling of $2^{n}$ collocation points, with $n=13$ on the boundary, $n=14$
inside, and $n=15$ outside of the particle, respectively. The numbers of
collocation points were selected to accommodate the GPU memory. The gradients
of the loss function with respect to model parameters ($d_{i}$ and
$\boldsymbol{\xi}_{i}$) were computed using automatic differentiation [39];
they were subsequently applied to update the model parameters by using the
ADAM optimizer [40] with a starting learning rate of $10^{-3}$. At each
$2,000$ iterations, the learning rates were adjusted to be the twice and half
of the current loss magnitude for $d_{i}$ and $NN_{i}$, respectively. The
collocation points were re-sampled after each $10,000$ iterations, suggesting
an unsupervised multi-task learning using mini-batch gradient descent with an
infinite set of collocation points. Our numerical models were implemented
using Python language and TensorFlow backend [41]. During the training, we
keep $w_{g_{1}}=w_{b_{n}}=w_{b_{t}}=1$ constant, but with the value of
$w_{g_{0}}\in(0,1]$ varies for different cases. The selection of $w_{g_{0}}$
is detailed in Appendix B. As a reference, on a workstation equipped with two
Nvdia GeForce GTX 1080 Ti graphics cards, each iteration takes around $0.1$
second.

Figure 3:  Dielectric $(a)$ oblate $(a=2/3,b=c=1)$ and $(b)$ prolate
$(a=3/2,b=c=1)$ particles with $\epsilon_{r}=6$ placed in a homogeneous field
$\boldsymbol{E}_{\text{ext}}=\hat{\boldsymbol{z}}$. From top to bottom: the
electric potential $\phi$, the electric fields $E_{x}/E_{\text{ext}}$ and
$E_{z}/E_{\text{ext}}$, and a distribution of $E_{z}/E_{\text{ext}}$ along
principal axes on $xz$-plane spanned by $r$ and $\theta$. As a reminder, the
exact electric potential outside a dielectric spheroid ($a,b=c=1$) is:
$\phi_{0}=-E_{0}z\left[1+(\epsilon_{r}-1)n_{z}^{\xi}\right]/\left[1+(\epsilon_{r}-1)n_{z}^{\infty}\right]$,
where
$n_{z}^{a}=(a/2)\int_{0}^{a}\mathrm{d}s/\left[(s+1)^{2}\sqrt{s+a^{2}}\right]$
and $\xi=z^{2}-1$ along $z$-axis, cf. [26].

The present model allows us to calculate the induced potential over the whole
domain. By differentiation, we obtain the total electric field,
$\boldsymbol{E}$, the polarization field
$\boldsymbol{P}=\frac{\epsilon_{r}-1}{4\pi}\boldsymbol{E}$ inside the
particle, and the source of it, which is the induced surface charge,
$\sigma=-\boldsymbol{P}\cdot\hat{\boldsymbol{n}}$, as illustrated in Fig.
2$(f)$. The polarizability,
$\boldsymbol{\alpha}=\boldsymbol{p}/E_{\text{ext}}$, where $\boldsymbol{p}$ is
the induced moment given by

$\boldsymbol{p}=\frac{\epsilon_{r}-1}{4\pi}\int_{V}\boldsymbol{E}\mathrm{d}V.$
(18)

Considering our problem symmetry, all terms but the dipole one cancel upon
integration over volume, leading to:

$p_{x}=p_{y}=0\quad\mbox{and}\quad
p_{z}=\frac{\epsilon_{r}-1}{4\pi}E_{z}|_{r=0}V,$ (19)

where $E_{z}|_{r=0}$ is the amplitude of the electric field at the center of
the particle. In order to compare with previous works, we introduce the volume
normalized polarizability defined by Sihvola et al. [30] as:
$\overline{\alpha}_{j}=4\pi\alpha_{j}/V$, with $j=x,y,z$. Thus,

$\displaystyle\overline{\alpha}_{z}\approx(\epsilon_{r}-1)\frac{E_{z}|_{r=R_{\text{min}}}}{E_{\text{ext}}},$
(20)

because $r=0$ is excluded from the solution domain.

In the following, we fix the parameters $E_{\text{ext}}=R=b=c=1$, and let
$R_{\text{min}}=0.01$. In order to keep the dipole contributions on the
boundaries of the solution domain of the same order, we take
$R_{\text{max}}=10$. The remaining three parameters: $N$, $a$, and
$\epsilon_{r}$, enable one to assess the emergence of edges and corners, the
stretching/squeeze of the geometry, and the dielectric mismatch at the
interface.

Let us first assess the present method with dielectric spheroids for which
there are analytical solutions and then apply it to quasi-cubes with
increasing values of $N$.

###  3.1 Spheroid

For $N=1$, Eq. (4) defines a unit sphere for $a=1$, and it deforms into a
spheroid as $a$ deviates from $a=1$. In both cases, the induced electric field
inside the dielectric particle is uniform and given by an analytical
expression [20, 24, 25, 26], thereby providing a benchmark for the accuracy of
the proposed ANN approach. Here, we consider a sphere ($a=1$), an oblate shape
with $a=2/3$ and a prolate one with $a=3/2$.

The training process is stopped when the loss function (15) drops below
$10^{-5}$. From the electric potential $\phi(r,\theta,\varphi)$, we calculate
the Cartesian components of the electric field $E_{x}$, $E_{y}$ and $E_{z}$ by
means of the chain rule. Their distributions inside and outside an oblate
spheroid with $a=2/3$ and a prolate spheroid with $a=3/2$ with relative
dielectric constant $\epsilon_{r}=6$ are plotted in Fig. 3. Indeed, the
calculated induced electric field is quite uniform inside the dielectric
particle and decays towards $E_{\text{ext}}$ outside, with values in good
agreement with theory [20, 24, 26].

In Fig. 4 we show the distribution of charges on the dielectric interface for
spheroidal particles. The accumulation of positive charges on the north pole
and negative charges on the south pole leads to an induced electric field
which counteracts $\boldsymbol{E}_{\text{ext}}$. Our calculated surface charge
distribution agrees nicely with the exact solution except for a small
deviation at the tips of the oblate spheroid in Fig. 4$(d)$, where the high
curvature in the $xz$-plane leads to large variation of $L_{ge}$ components of
$\phi_{0}$, degrading the numerical accuracy.

Figure 4:  Visualization of surface charge distribution for: $(a)$ an oblate
$(a=2/3,b=c=1)$, $b$ a spherical $(a=1,b=c=1)$, and $(c)$ a prolate
$(a=3/2,b=c=1)$ dielectric particles with $\epsilon_{r}=6$. Comparison of the
obtained surface charge distributions with their exact solutions on a cross
section with $(d)$ $xz$-plane $(\varphi=0)$ and $(e)$ $yz$-plane
$(\varphi=\pi/2)$. Table 1: Calculated normalized polarizabilities of
dielectric spheroidal particles compared with the corresponding exact values
and relative errors.

| $\epsilon_{r}=2$ | $\epsilon_{r}=6$ | $\epsilon_{r}=10$
---|---|---|---
This work | Exact | Error(%) | This work | Exact | Error(%) | This work | Exact | Error(%)
$L\leq 10^{-4}$ | Oblate | 0.7845 | 0.78306 | 0.185 | 2.1028 | 2.09623 | 0.313 | 2.5888 | 2.57627 | 0.486
Sphere | 0.7492 | 0.75000 | 0.107 | 1.8711 | 1.87500 | 0.208 | 2.2521 | 2.25000 | 0.093
Prolate | 0.7245 | 0.72280 | 0.235 | 1.7179 | 1.71377 | 0.241 | 2.0391 | 2.02175 | 0.858
$L\leq 10^{-5}$ | Oblate | 0.78355 | 0.78306 | 0.063 | 2.09661 | 2.09623 | 0.017 | 2.57710 | 2.57627 | 0.032
Sphere | 0.75001 | 0.75000 | 0.001 | 1.87507 | 1.87500 | 0.004 | 2.25008 | 2.25000 | 0.004
Prolate | 0.72314 | 0.72280 | 0.047 | 1.71379 | 1.71377 | 0.001 | 2.02157 | 2.02175 | 0.009

To obtain the polarizability, we compute the integral in Eq. (18) using a
Monte Carlo method [42]. The comparison is made for two values reached by the
loss function $L$ during optimization. As a reminder, the normalized
polarizability (columns “Exact") of a spheroid with axis ($a,b=c=1$) is:
$\overline{\alpha}_{z}=(\epsilon_{r}-1)/\left[1+(\epsilon_{r}-1)n_{z}^{\infty}\right]$,
with $n_{z}^{\infty}$ defined in caption of Fig. 3, cf. [24]. Calculated
normalized polarizabilities are within less than $1\%$ their exact values
depending upon the loss function accuracy, see Table. 1. More precisely, a
$99\%$ accuracy is found for $L\leq 10^{-4}$ and it reaches $99.9\%$ when the
loss function is minimized below $L\leq 10^{-5}$. A further minimization to
$L\leq 10^{-6}$ would require a significant increase in computation time for a
marginal increase in accuracy. Therefore, we shall limit ourselves to a loss
level $L\leq 10^{-5}$ for applications to quasi-cubic inclusions in the next
section.

Figure 5:  Contours for the normalized electric field $E_{z}/E_{\text{ext}}$
associated with dielectric quasi-cubes ($N=6$) with $(a)$ $\epsilon_{r}=2$ and
$(b)$ $\epsilon_{r}=8$ on the $xz$-plane. The insets reveal a saddle-shaped
field inside particles. Distribution of $E_{z}/E_{\text{ext}}$ along $(c)$
$x$\- and $(d)$ $z$-axes for dielectric particles with various values of
$\epsilon_{r}=2,4,6,8$.

###  3.2 Quasi-cube

For quasi-cubic particles with $N>1$, the emergence of edges and corners
strongly modifies the electric potential and its derivatives. Fig. 5 shows the
normalized induced electric field $E_{z}/E_{\text{ext}}$ in a quasi-cubic
particle with $N=6$ for values of $\epsilon_{r}=2,4,6,8$. The strong rise of
the electric field along the $z$-axis as one approaches the particle from
outside, its discontinuous drop at the interface, and the continuous decrease
of $E_{z}$ along the $x$-axis from infinity to the origin, all together lead
to a saddle-shaped field inside the particle, as shown in Fig. 5$(a,b)$. With
varying values of $\epsilon_{r}$, despite an apparent difference of amplitude,
the electric field remains qualitatively unchanged.

Figure 6:  Shape dependent polarizability for quasi-cubic particles with
$\epsilon_{r}=2,4,6,8$ where, as shown in the inset, points indicate values
computed using Eq. (21) and curves are fitted polynomial functions with an
asymptote $\tilde{\alpha}_{z}|_{N\to\infty}=1$. The shaded regions indicate
deviation of the computed data points to the fitted curve.

Figure 7:  $(a)$ Visualization of surface charge on dielectric quasi-cubic
particles with $\epsilon_{r}=6$ and $N=2,4,6,8$. $(b)$ Comparison of the
obtained surface charge distribution on cross sections $\varphi=0$ and
$\varphi=\pi/4$ of quasi-cubes with $N\in[1,10]$.

In Fig. 5$(c,d)$, it is observed that the induced electric field
$E_{z}/E_{\text{ext}}$ at the center of the particle is higher for quasi-cubes
than for a sphere, which implies a higher value of the normalized
polarizability for quasi-cubes than that for a sphere, cf. Eq. (20). To
quantify the relation between the polarizability and the shape of the
particle, we plot

$\displaystyle\tilde{\alpha}_{z}(N)=\frac{\overline{\alpha}_{z}(N)-\overline{\alpha}_{z}^{1}}{\overline{\alpha}_{z}^{\infty}-\overline{\alpha}_{z}^{1}},$
(21)

as a function of $N$ in Fig. 6. In most cases, $\overline{\alpha}_{z}$
obtained using the integral (18) is systematically smaller than a direct
evaluation using Eq. (20) by $0.2\%$. This small deviation can arise from
multiple origins, e.g. numerical accuracy, finite size effect, nonzero dipolar
component in neural networks, complicating the task of sourcing out the prime
contribution. In Eq. (21), the polarizability of quasi-cubic particles is re-
scaled by that of a sphere $\overline{\alpha}^{1}_{z}$ and that of a cube
$\overline{\alpha}_{z}^{\infty}$ to enable a comparison with different values
of $\epsilon_{r}$; hence, $\tilde{\alpha}_{z}=0$ for a sphere and
$\tilde{\alpha}_{z}=1$ for a cube, independent of $\epsilon_{r}$. For the
dependence of $\overline{\alpha}_{z}^{\infty}$ upon $\epsilon_{r}$ we use the
approximation formula given in Ref. [30]. Then, the obtained polarizabilities
are fitted to a polynomial function of $N$ which displays an asymptotic
behavior $\tilde{\alpha}_{z}(N)\to 1$ as $N\to\infty$. It is observed that
$\tilde{\alpha}_{z}(N)$ exhibits a rapid transition to almost $80\%$ and
$90\%$ of the asymptotic value for $N\leq 10$ and $N\leq 20$, respectively.
This numerical result is in accordance with the fact that the geometry of Eq.
(4) converges virtually to a cube for $N>8$, as seen in Fig. 1. Note that the
transition of $\tilde{\alpha}_{z}(N)$ to its asymptotic value is slower as
$\epsilon_{r}$ increases.

In order to investigate the behavior of induced surface charges with the
emergence of ever sharper edges and corners, we visualize $\sigma$ for quasi-
cubic particles with fixed $\epsilon_{r}=6$ and increasing values of
$N\in[1,10]$ in Fig. 7. Compared with the sinusoidal charge distribution on
the surface of a sphere, charges accumulate towards edges and are peaked at
corners as $N$ increases. Such a re-distribution of surface charges leads to
an enhanced dipole moment, which in turns leads to the higher value of
polarizability observed in Fig. 6, thereby a higher electric field at the
center of dielectric quasi-cubic particles observed in Fig. 5.

As a baseline comparison, we compare our ansatz with a vanilla ansatz. For the
latter, we replaced $H_{0}$ and $H_{1}$, cf. Eqs. (9), by multilayer
perceptrons with the same network architecture as $\text{NN}_{i}$ but a linear
activation at the output layer. The input variables $r$, $\theta$ and
$\varphi$ were rescaled to $[-1,1]$ using the min-max normalization. As
explicitly shown in Fig. 8 $(a)$ and $(b)$, the inclusion of physical
constraints in ansatz (9) comparatively reduces the optimization effort by a
large degree. It is worth noting that in order to capture the transition
regime from sphere to cube, classical numerical methods, e.g. finite element,
would require a re-mesh and a re-simulation for each different values of $N$.
Alternatively, the proposed ANN approach is able to tackle more efficiently
this progressive transition. Since each change in $N$ will only affect the
boundary conditions and few collocation points near the dielectric interface,
while Laplace’s equation remains satisfied in the bulk of the solution domain,
a minimization initiated from a previous converged solution for, e.g. $N-1$,
leads to a drastic reduction in computational time as shown in Fig. 8 $(c)$
and $(d)$. For piece-wise homogeneous dielectric media considered in this
work, $\epsilon_{r}$ appears only through the normal boundary condition (2b).
Thus, for simulations with different values of $\epsilon_{r}$ and a fixed
value of $N$, neural networks are committed to learn the change in mismatch,
which is induced by a continuous variation of $\epsilon_{r}$, at the
dielectric interface, cf. Fig. 8 $(e)$ and $(f)$. Once the loss function is
minimized below the target value, only the converged solutions, which consist
of the model variables and the structure of the ANN ansatz, are stored
independently of collocation points. Unlike finite element methods whose
solutions are defined on meshes, our ANN approach provides a continuous
mapping: $(r,\theta,\varphi)\to\phi$ over the entire solution domain
$r\in[R_{\text{min}},R_{\text{max}}]$, resembling analytical solutions.

Figure 8:  Decay of loss function $L$ for the case $N=3$ and $\epsilon_{r}=4$
initiated from $(a)$ random conditions using the proposed ansatz (8); $(b)$
random conditions using the vanilla ansatz; $(c)$ converged solution of $N=2$
and $\epsilon_{r}=4$; $(d)$ converged solution of $N=4$ and $\epsilon_{r}=4$;
$(e)$ converged solution of $N=3$ and $\epsilon_{r}=2$; $(f)$ converged
solution of $N=3$ and $\epsilon_{r}=6$. Iteration stops at $L\leq 10^{-4}$ or
after $150,000$ iterations. In $(a)$ $L$ decays below $10^{-4}$ within
$60,000$ iterations; whereas in $(b)$ $L$ remains of order of $10^{-3}$ after
$150,000$ iterations. In $(a)$ and $(b)$, the first $100$ iterations are
removed for a better visualization. Due to the slow convergence in $(b)$, the
learning rate is decreased by a factor of $1.2$ at every $2000$ iterations
until a minimal learning rate $5\times 10^{-5}$ is reached.

##  4 Conclusions

Recent advances in machine learning enable a revisit of longstanding physics
problems from a new perspective. We have presented a neural networks based
calculation for the response of dielectric particles with shapes varying from
sphere to cube, placed in an external uniform electric field. Our ansatz
intertwined boundary conditions at the borders of the solution domain, as well
as symmetry constraints resulting from the geometry of the particle and the
external field, with neural networks. Then, solving Laplace’s equation with
dielectric boundary conditions was translated to a minimization of a loss
function defined in Sec. 2.3. To evaluate the accuracy, we applied the
proposed ANN approach to spheroids with various aspect ratios and relative
dielectric constants. An overall $99.9\%$ accuracy for the normalized
polarizability was achieved, with however, a slight deviation of surface
charge from the exact solution at the tips of an oblate spheroid with large
curvature, as shown in Fig. 4$(d)$.

As a sphere progressively deforms into a cube, the accumulation of surface
charges towards the ever sharper edges and corners leads to a rapid transition
of polarizability to its asymptotic value at lower values of $N\leq 10$. This
implies that the shape effect has a significant impact on determining the
induced polarizability by dielectric nano-particles. The enhanced
polarizability with increasing values of $N$ leads to an amplified dipole
moment, which, in turn, results in a higher electric field at the center of a
quasi-cubic particle than that of a sphere, independent of the relative
dielectric constant of the particle.

Since neural networks are infinitely differentiable, instead of interpolating
among discrete values, auto-differentiation enables a conversion of the
obtained ANN solution to higher order derivatives, which are again continuous
functions. This feature can be advantageous in a broad range of physical
applications where higher-order derivatives are of interest. By contrast to
the finite element method, where the entire mesh is required for computation,
the ANN methods during training use only a fraction of collation points to
compute the loss function and update the model parameters by mini-batch
stochastic gradient descent techniques. A successive re-sampling enables a
complete covering of the solution domain with sufficient number of iterations.
Therefore, the mesh-free ANN approach overcomes two major deficiencies of the
finite element method stemming from the mesh dependence: (i) the finite
differentiability; and (ii) the high memory usage associated with fine mesh.

Finally, it is worth emphasizing some limitations which need to be overcome in
future works. Loss functions stemming from physics problems often consist of
several components, e.g. governing equation, initial and boundary conditions,
symmetries and conservation laws. Since the loss function is not strictly
zero, an imbalance in its components can deteriorate the accuracy of the
solution. In this work, we introduced weights to balance the amplitudes of
each loss component during the training process. However, it is unclear, how
does the selection of these weights affect the numerical accuracy. Therefore,
it would be beneficial to establish a relation between the accuracy and the
relative amplitudes of each components, and devise an adaptive algorithm which
automatically adjusts the weights during the training process.

## Appendix

##  Appendix A Output activation

The linear $x$, $\tanh(x)$, and $g(x)=\text{atanh}(\eta\tanh(x))$ activation
functions and their derivatives are shown in Fig. 9. By varying
$\eta\in(0,1)$, we adjust the output of ANN to an arbitrary bounded interval.
Alternatively, one can achieve the same purpose by scaling the $\tanh$ output.
However, as observed from Fig. 9$(b)$ that a scaling up of $\tanh$ output
leads to a sharp variation of gradients. Therefore, we select the activation
function at the output layer to be $g(x)$.

Figure 9:  Comparison of $(a)$ linear $x$, $\tanh$, and proposed $g(x)$
activation functions; and $(b)$ their derivatives with respect to $x$.

##  Appendix B Weight tuning

The weight $w_{g_{0}}$ is selected to balance the components of $L_{ge}$ on
the dielectric interface. It is observed from Fig. 10 that, due to the
emergence of ever sharper corners with increasing $N$, the components of
$L_{ge}$ inside and outside of the particle differ in magnitudes
substantially. To enable an optimization of $\phi_{0}$ and $\phi_{1}$ on the
same footing, we take the spherical case $N=1$ as a reference and select
$w_{g_{0}}=0.3,0.2,0.15$ for $N=3,6,9$, respectively. However, since the
relative magnitudes of each loss component is not known a priori, an adaptive
algorithm which adjusts $w_{g_{0}}$ during the training can be beneficial.

Figure 10:  Visualization of the radial (blue), the polar (orange) and the
azimuthal (green) components of $L_{ge}$ against $\theta$ on the interface
with $\varphi=\pi/4$ and $\epsilon_{r}=4$. The solid and the dashed lines are
associated with $\phi_{0}$ and $\phi_{1}$, respectively. The gray vertical
lines mark the locations of corners for a canonical cube.

## Acknowledgment

The authors would like to thank Steven Blundell and Tan Nguyen for discussions
that initiated this work, Liang Mong Koh, Sean Ooi and Kavitha Srinivasan for
providing computational resources and Subodh Mhaisalkar for support. The
computational work for this article was partially performed on resources of
the National Supercomputing Centre, Singapore (https://www.nscc.sg).

## References

  * [1] A. N. Gorban and D. C. Wunsch, “The general approximation theorem,” in _Proceedings of the International Joint Conference on Neural Networks_ , 1998\.
  * [2] D. A. Winkler and T. C. Le, “Performance of deep and shallow neural networks, the universal approximation theorem, activity cliffs, and qsar,” _Molecular Informatics_ , vol. 36, no. 1-2, p. 1600118, 2017.
  * [3] H. Lin and S. Jegelka, “Resnet with one-neuron hidden layers is a universal approximator,” in _Advances in Neural Information Processing Systems_ , 2018\.
  * [4] A. J. Meade Jr. and A. A. Fernandez, “The numerical solution of linear ordinary differential equations by feedforward neural networks,” _Mathematical and Computer Modelling_ , vol. 19, no. 12, pp. 1–25, 1994.
  * [5] ——, “Solution of nonlinear ordinary differential equations by feedforward neural networks,” _Mathematical and Computer Modelling_ , vol. 20, no. 9, pp. 19–44, 1994.
  * [6] I. E. Lagaris, A. Likas, and D. I. Fotiadis, “Artificial neural networks for solving ordinary and partial differential equations,” _IEEE Transactions on Neural Networks_ , vol. 9, no. 5, pp. 987 – 1000, 1998.
  * [7] I. E. Lagaris, A. C. Likas, and D. G. Papageorgiou, “Neural-network methods for boundary value problems with irregular boundaries,” _IEEE Transactions on Neural Networks_ , vol. 11, no. 5, pp. 1041 – 1049, 2000.
  * [8] K. S. McFall and J. R. Mahan, “Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions,” _IEEE Transactions on Neural Networks_ , vol. 20, no. 8, pp. 1221 – 1233, 2009.
  * [9] M. Raissi, “Deep hidden physics models: Deep learning of nonlinear partial differential equations,” _Journal of Machine Learning Research_ , vol. 19, no. 25, pp. 1–24, 2018.
  * [10] Y. Yang and P. Perdikaris, “Adversarial uncertainty quantification in physics-informed neural networks,” _Journal of Computational Physics_ , vol. 394, pp. 136–152, 2019.
  * [11] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, “Universal differential equations for scientific machine learning,” arXiv:2001.04385.
  * [12] M. Raissi, Z. Wang, M. S. Triantafyllou, and G. E. Karniadakis, “Deep learning of vortex-induced vibrations,” _Journal of Fluid Mechanics_ , vol. 861, pp. 119–137, 2018.
  * [13] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 378, pp. 686–707, 2019.
  * [14] M. Raissi, A. Yazdani, and G. E. Karniadakis, “Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations,” _Science_ , vol. 367, no. 6481, pp. 1026–1030, 2020.
  * [15] Y. Chen, L. Lu, G. E. Karniadakis, and L. D. Negro, “Physics-informed neural networks for inverse problems in nano-optics and metamaterials,” _Optics Express_ , vol. 28, no. 8, p. 11618, 2020.
  * [16] J. C. Wong, A. Gupta, and Y. S. Ong, “Can transfer neuroevolution tractably solve your differential equations?” arXiv:2101.01998.
  * [17] S. Wang, X. Yu, and P. Perdikaris, “When and why pinns fail to train: A neural tangent kernel perspective,” arXiv:2007.14527.
  * [18] S. Wang, Y. Teng, and P. Perdikaris, “Understanding and mitigating gradient pathologies in physics-informed neural networks,” arXiv:2001.04536.
  * [19] M. Elhamod, J. Bu, C. Singh, M. Redell, A. Ghosh, W. C. L. V. Podolskiy, and A. Karpatne, “Cophy-pgnn: Learning physics-guided neural networks withcompeting loss functions for solving eigenvalue problems,” arXiv:2007.01420.
  * [20] V. Klimov, _Nanoplasmonics_. CRC Press, Taylor and Francis Group, 2013.
  * [21] C. Delerue and M. Lannoo, _Nanostructures: theory and modeling_. Springer Science & Business Media, 2013.
  * [22] M. A. Becker et al., “Bright triplet excitons in caesium lead halide perovskites,” _Nature_ , vol. 553, no. 1, pp. 189–193, 2018.
  * [23] T. P. T. Nguyen, S. A. Blundell, and C. Guet, “One-photon absorption by inorganic perovskite nanocrystals: a theoretical study,” _Physical Review B_ , vol. 101, p. 195414, 2020.
  * [24] L. D. Landau, E. M. Lifshitz, and L. P. Pitaevshiǐ, _Electrodynamics of continuous media_. Pergamon Press, 1984\.
  * [25] J. D. Jackson, _Classical Electrodynamics_. John Wiley & Sons, 2007.
  * [26] J. A. Stratton, _Electromagnetic Theory_. John Wiley & Sons, 2007.
  * [27] T. W. Edwards and J. V. Bladel, “Electrostatic dipole moment of a dielectric cube,” _Applied Scientific Research, Section B_ , vol. 9, no. 9, pp. 151–154, 1961.
  * [28] D. F. Herrick and T. A. Senior, “The dipole moment of a dielectric cube,” _IEEE Transactions on Antennas and Propagation_ , vol. 25, pp. 590–592, 1977\.
  * [29] L. Eyges and P. Gianino, “Polarizabilities of rectangular dielectric cylinders and of a cube,” _IEEE Transactions on Antennas and Propagation_ , vol. 27, pp. 557–560, 1979.
  * [30] J. Avelin, H. H. R. Sharma, and A. H. Sihvola, “Polarizability analysis of cubical and square-shaped dielectric scatterers,” _IEEE Transactions on Antennas and Propagation_ , vol. 49, no. 3, pp. 451–457, 2001.
  * [31] A. Sihvola, P. Ylä-Oijala, S. Järvenpää, and J. Avelin, “Polarizabilities of platonic solids,” _IEEE Transactions on Antennas and Propagation_ , vol. 52, no. 3, pp. 451–457, 2004.
  * [32] A. Sihvola, “Dielectric polarization and particle shape effects,” _Journal of Nanomaterials_ , vol. 2007, pp. 1–9, 2007.
  * [33] J. Helsing and K.-M. Perfekt, “On the polarizability and capacitance of the cube,” _Applied and Computational Harmonic Analysis_ , vol. 34, p. 445–468, 2013.
  * [34] M. V. Kovalenko, L. Protesescu, and M. I. Bodnarchuk, “Properties and potential optoelectronic applications of lead halide perovskite nanocrystals,” _Science_ , vol. 358, no. 6364, pp. 745–750, 2017.
  * [35] Y. Tang et al., “Highly stable perovskite supercrystals via oil-in-oil templating,” _Nano Letters_ , vol. 20, no. 8, p. 5997–6004, 2020.
  * [36] A. A. Khelashvili and T. P. Nadareishvili, “Singular behavior of the laplace operator in polar spherical coordinates and some of its consequences for the radial wave function at the origin of coordinates physics of particles and nuclei letters,” _Physics of Particles and Nuclei Letters_ , vol. 12, pp. 11–25, 2015.
  * [37] G. Cybenko, “Approximation by superpositions of a sigmoidal function,” _Mathematics of Control, Signals and Systems_ , vol. 2, pp. 303–314, 1989\.
  * [38] R. van der Meer, C. Oosterlee, and A. Borovykh, “Optimally weighted loss functions for solving pdes with neural networks,” _arXiv:2002.06269_ , 2020\.
  * [39] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, “Automatic differentiation in machine learning: a survey,” _Journal of Machine Learning Research_ , vol. 18, no. 1, pp. 5595–5637, 2017.
  * [40] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _International Conference on Learning Representations_ , 2015.
  * [41] M. Abadi et al., “Tensorflow: Large-scale machine learning on heterogeneous distributed systems,” _arXiv:1603.04467_ , 2016.
  * [42] R. Y. Rubinstein and D. P. Kroese, _Simulation and the Monte Carlo method_. John Wiley & Sons, 2016.


# Aperiodic subshifts of finite type on groups which are not finitely
generated

Sebastián Barbieri

###### Abstract

We provide an example of a non-finitely generated group which admits a
nonempty strongly aperiodic SFT. Furthermore, we completely characterize the
groups with this property in terms of their finitely generated subgroups and
the roots of their conjugacy classes.

_Keywords: symbolic dynamics, subshift of finite type, aperiodicity._

_MSC2020:_ 37B10.

Let $A$ be a finite set and $G$ a group with identity $1$. Consider the set of
configurations $A^{G}=\\{x\colon G\to A\\}$ endowed with the left shift action
$G\curvearrowright A^{G}$ given by

$(gx)(h)=x(g^{-1}h)\mbox{ for every }x\in A^{G}\mbox{ and }g,h\in G.$

A set $X\subset A^{G}$ is called a subshift of finite type (SFT) if there
exists a finite set $F\subset G$ and $L\subset A^{F}$ such that $x\in X$ if
and only if for every $g\in G$ we have that the restriction of $gx$ to $F$ is
in $L$. An SFT can be though of as a set of colorings of $G$ using colors from
$A$ which satisfy a finite set of local rules encoded by $L$.

An SFT $X$ is called strongly aperiodic (SA) if the shift action is free,
namely, given $x\in X$ we have that $gx=x$ can only hold for $g=1$. Clearly
$X=\varnothing$ is SA, but this case is not very interesting. This definition
begs the question of whether nonempty SA SFTs actually exist and the answer
turns out to depend on which group $G$ we are considering.

Let us illustrate the failure of the existence of nonempty SA SFT in the case
$G=\mathbb{Z}$. Given $A,F$ and $L$ which define an SFT $X$ we may assume
without loss of generality that $F$ is of the form $\\{0,\dots,n-1\\}$ for
some $n\geq 1$. Consider the finite directed graph with vertices $A^{n+1}$ and
such that there is an edge from $u_{0}\dots u_{n}$ to $v_{0}\dots v_{n}$ if
and only if $u_{1}\dots u_{n}=v_{0}\dots v_{n-1}\in L$. It is clear that the
elements of $X$ are precisely the projections of the bi-infinite walks in this
graph to the first symbol in each coordinate (for further details, see [17]).
If $X$ is nonempty, it follows that the graph must contain a cycle, say of
length $m>0$, and consequently there must exist $x\in X$ for which $mx=x$,
showing that $X$ is not SA.

The case of $G=\mathbb{Z}^{2}$ is much more interesting. In [20] Wang studied
the algorithmic problem of deciding, given as inputs $A,F$ and $L$ whether the
SFT $X$ they define is empty or not111To be rigorous, Wang studied an
equivalent problem where the SFT is described by square tiles with colored
edges and in any tiling of $\mathbb{Z}^{2}$ the colors at the edges must
match.. Wang did not solve this problem but showed that if $\mathbb{Z}^{2}$
does not admit any nonempty SA SFT, then the problem would be algorithmically
decidable. Years later Berger solved the problem [7] showing that it was
algorithmically undecidable, and in doing so constructed the first known
example of a nonempty SA SFT in $\mathbb{Z}^{2}$. After this breakthrough,
several beautiful SA SFTs on $\mathbb{Z}^{2}$ have been constructed, see [15,
16, 19].

In recent years the problem of classifying the groups which admit a nonempty
SA SFT has gained notoriety. This has led to wonderful discoveries that relate
this property with the algorithmic and geometric properties of groups. For
instance, the property of admitting a nonempty SA SFT is invariant under
commensurability [8] and under quasi-isometries of finitely presented groups
[9]. It is also known that finitely generated groups with infinitely many ends
(in particular, virtually free groups) cannot admit nonempty SA SFTs [9], and
that finitely generated and recursively presented groups which admit SA SFTs
must necessarily have decidable word problem [13]. Within the groups that
satisfy these constraints, several are known to admit SA SFTs. For instance
polycyclic-by-finite groups which are not virtually cyclic [3, 14], one-ended
word-hyperbolic groups [10, 11], some Baumslag-Solitar groups [1, 2, 12],
groups of the form $\mathbb{Z}^{d}\rtimes_{\varphi}G$ with $d\geq 2$,
$\varphi\in\operatorname{GL}_{d}(\mathbb{Z})$ and $G$ infinite, finitely
generated and with decidable word problem [5], groups which are the direct
product of three infinite, finitely generated groups with decidable word
problem [4], the Grigorchuk group, and more generally, any finitely generated
branch group with decidable word problem [4], and any self-simulable group
with decidable word problem such as Thompson’s $V$,
$\operatorname{GL}_{n}(\mathbb{Z})$ and $\operatorname{SL}_{n}(\mathbb{Z})$
for $n\geq 5$, or the direct product of any two finitely generated non-
amenable groups with decidable word problem [6].

A common point shared by all the existence results mentioned above is that
they apply to finitely generated groups. For obvious reasons this is not very
surprising: if $X\subset A^{G}$ is an SFT on a group $G$ described by
$L\subset A^{F}$ and $H=\langle F\rangle$ is the subgroup of $G$ finitely
generated by $F$, we may consider instead $Y\subset A^{H}$ as the SFT which is
defined by $L$ but on $H$. It turns out that configurations in $X$ can be seen
as independent copies of configurations of $Y$ on each coset of $H$.

######  Proposition 1.

Let $H\leqslant G$, $F\subset H$ finite, $L\subset A^{F}$ and let $X,Y$ be the
SFTs defined by $L$ on $G$ and $H$ respectively. Let $(g_{i})_{i\in I}$ be a
set of left coset representatives of $H$ in $G$. Then $x\in X$ if and only if
for every $i\in I$ we have $y_{i}=(x(g_{i}h))_{h\in H}\in Y$.

###### Proof.

Let $x\in X$ and $i\in I$. For every $u\in H$ we have
$(uy_{i})(h)=y_{i}(u^{-1}h)=x(g_{i}u^{-1}h)=(ug_{i}^{-1}x)(h)$. As $x\in X$,
it follows that $(ug_{i}^{-1}x)|_{F}\in L$ and thus $(uy_{i})|_{F}\in L$ and
$y_{i}\in Y$. Conversely, suppose $y_{i}\in Y$ for every $i\in I$ and let
$g\in G$. We may choose $i\in I$ and $v\in H$ such that $g^{-1}=g_{i}v$ and
thus we have
$(gx)(h)=(v^{-1}g_{i}^{-1}x)(h)=x(g_{i}vh)=y_{i}(vh)=(v^{-1}y_{i})(h)$ for
every $h\in H$. As $y_{i}\in Y$, it follows that $v^{-1}y_{i}|_{F}\in L$ and
thus $gx|_{F}\in L$, showing that $x\in X$. ∎

Based on the former argument, one should expect that non-finitely generated
groups would have a hard time admitting nonempty SA SFTs, as one could always
pick the same configuration on each coset and use that to create periodicity.
Let us exemplify that this intuition is wrong with a real world situation.

Let us suppose the reader is walking on the street, while suddenly they get
jumped by a thief who demands: “Quickly! give me an example of a non-finitely
generated group which admits a nonempty SA SFT or I shall take your wallet!”.
The reader, scared by the gravity of the situation, might wrongly answer that
such examples do not exist by explaining the argument above. The thief, with a
victorious grin on their face, would take out their portable blackboard and
write down the following example.

######  Example 2.

Consider $G=\mathbb{Q}^{2}$ and let $Y\subset A^{\mathbb{Z}^{2}}$ be a
nonempty SA SFT. Let $X\subset A^{\mathbb{Q}^{2}}$ be given by the condition

$x\in X\mbox{ if and only if }(sx)|_{\mathbb{Z}^{2}}\in Y,\mbox{ for every
}s\in\mathbb{Q}^{2}.$

As $Y$ is a nonempty SFT, it follows that $X$ is a nonempty SFT. For
$q\in\mathbb{Q}^{2}$, we may write
$q=\left(\frac{p_{1}}{r_{1}},\frac{p_{2}}{r_{2}}\right)$ with
$p_{1},p_{2},r_{1},r_{2}\in\mathbb{Z}$ and $r_{1}r_{2}\neq 0$. It follows that
$r_{1}r_{2}q\in\mathbb{Z}^{2}$. Now suppose we have $qx=x$, then we also have
$(r_{1}r_{2}q)x=x$ and from the fact that $x|_{\mathbb{Z}^{2}}\in Y$, we
obtain that necessarily $r_{1}r_{2}q=(0,0)$. As $r_{1}r_{2}\neq 0$, we deduce
that $q=(0,0)$ and thus $X$ is SA. $\ocircle$

Faced with this example, the reader would have no other choice but to
surrender their wallet to the thief. The argument for nonexistence sketched
above is incomplete: while it is true that one may choose independently any
configuration in every coset, it may happen that powers of the coset
representatives end up inevitably falling on a non-trivial element of the
finitely generated subgroup and thus destroy any global translational
symmetry. As we shall show, it turns out that modulo conjugacy this is
essentialy the sole way that SA SFTs can arise in non-finitely generated
groups. In fact, we shall show that global aperiodicity may arise even if the
subshift in the finitely generated subgroup is not strongly aperiodic.

For an SFT $X\subset A^{G}$, define its free part as the set

$\operatorname{Free}(X)=G\setminus\bigcup_{x\in
X}\operatorname{Stab}_{G}(x)=\\{g\in G:gx\neq x\mbox{ for every }x\in X\\}.$

In particular, a nonempty SFT $X$ is SA if and only if
$\operatorname{Free(X)}=G\setminus\\{1\\}$. For a subset $M\subset G$, let us
denote the set of its roots in $G$ by

$R_{G}(M)=\\{g\in G:\mbox{ there is }n>0\mbox{ such that }g^{n}\in M\\}.$

Finally, given $g\in G$ let us denote its conjugacy class by
$\operatorname{Cl}(g)=\\{tgt^{-1}:t\in G\\}$.

######  Proposition 3.

Let $H\leqslant G$, $F\subset H$ finite, $L\subset A^{F}$ and let $X,Y$ be
nonempty SFTs defined by $L$ on $G$ and $H$ respectively. We have that
$g\in\operatorname{Free}(X)$ if and only if

$\operatorname{Cl}(g)\cap R_{G}(\operatorname{Free}(Y))\neq\varnothing.$

###### Proof.

If $g\notin\operatorname{Free}(X)$, then there exists $x\in X$ such that
$gx=x$. Suppose there is $n>0$, $h\in H$ and $t\in G$ such that
$g^{n}=t^{-1}ht$. It follows that we would have $t^{-1}htx=g^{n}x=x$ and thus
$htx=tx$. Letting $z=tx$, it follows that $hz=z$. As $z\in X$, it follows that
if we let $y=z|_{H}\in Y$ we have $hy=y$ and thus
$h\notin\operatorname{Free}(Y)$. This shows that $\operatorname{Cl}(g)\cap
R_{G}(\operatorname{Free}(Y))=\varnothing$.

Conversely, let $g\in G$ be such that $\operatorname{Cl}(g)\cap
R_{G}(\operatorname{Free}(Y))=\varnothing$ and choose (using the axiom of
choice) a set $(\ell_{i})_{i\in I}$ of left coset representatives of $H$.
There is a unique well-defined permutation $\varphi\colon I\to I$ which
satisfies $g\ell_{i}=\ell_{\varphi(i)}h$ for some $h\in H$. In particular,
$\varphi$ induces an action of $\mathbb{Z}$ on $I$ given by $m\cdot
i=\varphi^{m}(i)$ for $m\in\mathbb{Z}$ and $i\in I$.

Choose (using the axiom of choice) $J\subset I$ such that it contains exactly
one representative of each orbit of $\varphi$, that is, for every $i\in I$
there is a unique $j\in J$ such that $i=\varphi^{m}(j)$ for some (non
necessarily unique) $m\in\mathbb{Z}$. For $n\in\mathbb{Z}$ and $j\in J$ let
$h_{n,j}=\ell_{\varphi^{n}(j)}^{-1}g^{n}\ell_{j}\in H.$

We will now define a configuration $x\in A^{G}$. Let $y^{*}\in Y$ be a fixed
configuration and let $j\in J$, there are two cases to consider:

  1. 1.

If $\\{\varphi^{m}(j)\\}_{m\in\mathbb{Z}}$ is infinite, we let $y_{j}=y^{*}$
and define $x(\ell_{\varphi^{m}(j)}s)=h_{m,j}y_{j}(s)=y_{j}(h^{-1}_{m,j}s)$
for every $s\in H$ and $m\in\mathbb{Z}$.

  2. 2.

If $\\{\varphi^{m}(j)\\}_{m\in\mathbb{Z}}$ is finite, let $n>0$ be the least
positive integer such that $\varphi^{n}(j)=j$. As
$h_{n,j}=\ell_{j}^{-1}g^{n}\ell_{j}\in\operatorname{Cl}(g^{n})$ it follows by
our assumption that $h_{n,j}\notin\operatorname{Free}(Y)$. Therefore there
exists $y_{j}\in Y$ such that $h_{n,j}y_{j}=y_{j}$. We define
$x(\ell_{\varphi^{m}(j)}s)=h_{m,j}y_{j}(s)=y_{j}(h_{m,j}^{-1}s)$ for every
$s\in H$ and $m\in\mathbb{Z}$. Notice that as $h_{n,j}y_{j}=y_{j}$, this is
well-defined.

By construction, we have that $(x(\ell_{i}h))_{h\in H}\in Y$ for every $i\in
I$ and thus by Proposition 1 we have that $x\in X$. Let us show that $gx=x$.
Indeed, for $j\in J,m\in\mathbb{Z}$ and $s\in H$ we have

$gx(\ell_{\varphi^{n}(j)}s)=x(g^{-1}\ell_{\varphi^{n}(j)}s)=x(\ell_{\varphi^{n-1}(j)}\ell_{\varphi^{n-1}(j)}^{-1}g^{-1}\ell_{\varphi^{n}(j)}s)$

As $\ell_{\varphi^{n-1}(j)}^{-1}g^{-1}\ell_{\varphi^{n}(j)}\in H$, we obtain
that
$x(\ell_{\varphi^{n-1}(j)}\ell_{\varphi^{n-1}(j)}^{-1}g^{-1}\ell_{\varphi^{n}(j)}s)=y_{j}(h_{n-1,j}^{-1}\ell_{\varphi^{n-1}(j)}^{-1}g^{-1}\ell_{\varphi^{n}(j)}s)$.
From here we obtain that

$gx(\ell_{\varphi^{n}(j)}s)=y_{j}(h_{n-1,j}^{-1}\ell_{\varphi^{n-1}(j)}^{-1}g^{-1}\ell_{\varphi^{n}(j)}s)=y_{j}(h_{n,j}^{-1}s)=x(\ell_{\varphi^{n}(j)}s)$

This shows that $gx=x$ and thus we conclude that
$g\notin\operatorname{Free}(X)$. ∎

The previous proposition tells us that we may determine whether a group admits
a nonempty SA SFT by just looking at its finitely generated subgroups.

######  Theorem 4.

A group $G$ admits a nonempty strongly aperiodic subshift of finite type if
and only if there exists a finitely generated subgroup $H\leqslant G$ and a
nonempty SFT $Y\subset A^{H}$ such that for every $g\in G\setminus\\{1\\}$ we
have

$\operatorname{Cl}(g)\cap R_{G}(\operatorname{Free}(Y))\neq\varnothing.$

###### Proof.

If $X\subset A^{G}$ is a nonempty SA SFT given by $F\subset G$ finite and
$L\subset A^{F}$, then $H=\langle F\rangle$ is a finitely generated subgroup
and by Proposition 3 for every $g\in\operatorname{Free}(X)=G\setminus\\{1\\}$
we have $\operatorname{Cl}(g)\cap
R_{G}(\operatorname{Free}(Y))\neq\varnothing$, where $Y\subset A^{H}$ is the
nonempty SFT given by $L$ on $H$. Conversely, if $H\leq G$ is a finitely
generated subgroup, $Y\subset A^{H}$ a nonempty SFT given by $F\subset H$
finite and $L\subset A^{F}$, and such that for every $g\in G\setminus\\{1\\}$
we have $\operatorname{Cl}(g)\cap
R_{G}(\operatorname{Free}(Y))\neq\varnothing$, then again by Proposition 3 we
have that $G\setminus\\{1\\}\subset\operatorname{Free}(X)$ for the nonempty
SFT $X\subset A^{G}$ induced by $L$ on $G$. Thus $X$ is SA. ∎

######  Corollary 5.

Let $G$ be a group and $H\leqslant G$ a finitely generated subgroup which
admits a nonempty strongly aperiodic SFT and such that for every $g\in
G\setminus\\{1\\}$ we have $\operatorname{Cl}(g)\cap
R_{G}(H\setminus\\{1\\})\neq\varnothing$. Then $G$ admits a nonempty strongly
aperiodic SFT.

Notice that for any nonempty subshift $Y\subset A^{H}$ we have
$\operatorname{Free}(Y)\subset H\setminus\\{1\\}$ and thus the condition
$\operatorname{Cl}(g)\cap R_{G}(\operatorname{Free}(Y))\neq\varnothing$ always
implies that $\operatorname{Cl}(g)\cap
R_{G}(H\setminus\\{1\\})\neq\varnothing$. This trivial remark suggests the
question of whether there are any examples where $X$ is SA but $Y$ is not. The
following example due to Salo (also in an article by Jeandel [13] with a
slightly different construction) shows that this case can indeed occur.

######  Example 6.

A result of Osin [18] shows that every countable torsion-free group can be
embedded in a $2$-generated group with exactly two conjugacy classes. In
particular there is a $2$-generated group $G$ with two conjugacy classes and
such that there is $t\in G$ with $\langle t\rangle\simeq\mathbb{Z}$. Let
$A=\\{\mathtt{a},\mathtt{b}\\}$, $F=\\{1,t\\}$ and
$L=\\{(1\mapsto\mathtt{a},t\mapsto\mathtt{b}),(1\mapsto\mathtt{b},t\mapsto\mathtt{a})\\}$,
that is, we let $X\subset\\{\mathtt{a},\mathtt{b}\\}^{G}$ be the SFT which
consists of all maps for which $x(gt)\neq x(g)$ for every $g\in G$.

Let $Y\subset\\{\mathtt{a},\mathtt{b}\\}^{\langle t\rangle}$ be the SFT
defined by $F$ and $L$ on $\langle F\rangle\simeq\mathbb{Z}$. On the one hand,
we have that $Y$ is nonempty and consists on the two periodic configurations
which alternate symbols, thus
$\operatorname{Free}(Y)=\\{t^{2n+1}:n\in\mathbb{Z}\\}$ consists of the odd
powers of $t$, in particular $Y$ is not SA. On the other hand, as there is
only one non-trivial conjugacy class, we have that for every $g\in
G\setminus\\{1\\}$, $t\in\operatorname{Cl}(g)\cap
R_{G}(\operatorname{Free}(Y))$ which is nonempty, and thus $X$ is SA by
Theorem 4. $\ocircle$

Acknowledgments: The author is grateful to Ville Salo for finding a subtle
error in the first version of this paper and suggesting the Osin group
example. The author also wishes to thank the referees for suggesting various
improvements. The author was supported by the FONDECYT grant 11200037.

S. Barbieri, Departamento de Matemática y ciencia de la computación,
Universidad de Santiago de Chile, Santiago. Chile.

E-mail address<EMAIL_ADDRESS>
## References

  * [1] N. Aubrun and J. Kari.  Tiling problems on Baumslag-Solitar groups.  Electronic Proceedings in Theoretical Computer Science, 128:35–46, 2013.
  * [2] N. Aubrun and J. Kari.  On the domino problem of the Baumslag-Solitar groups.  Theoretical Computer Science, 894:12–22, 2021.
  * [3] A. Ballier and M. Stein.  The domino problem on groups of polynomial growth.  Groups, Geometry, and Dynamics, 12(1):93–105, 2018.
  * [4] S. Barbieri.  A geometric simulation theorem on direct products of finitely generated groups.  Discrete Analysis, pages 1–25, Paper No. 9, 2019.
  * [5] S. Barbieri and M. Sablik.  A generalization of the simulation theorem for semidirect products.  Ergodic Theory and Dynamical Systems, 39(12):3185–3206, 2019.
  * [6] S. Barbieri, M. Sablik, and V. Salo.  Groups with self-simulable zero-dimensional dynamics.  arXiv:2104.05141, 2021.
  * [7] R. Berger.  The Undecidability of the Domino Problem.  American Mathematical Society, 1966.
  * [8] D. Carroll and A. Penland.  Periodic points on shifts of finite type and commensurability invariants of groups.  New York Journal of Mathematics, 21:811–822, 2015.
  * [9] D. B. Cohen.  The large scale geometry of strongly aperiodic subshifts of finite type.  Advances in Mathematics, 308:599–626, 2017.
  * [10] D. B. Cohen and C. Goodman-Strauss.  Strongly aperiodic subshifts on surface groups.  Groups, Geometry, and Dynamics, 11(3):1041–1059, 2017.
  * [11] D. B. Cohen, C. Goodman-Strauss, and Y. Rieck.  Strongly aperiodic subshifts of finite type on hyperbolic groups.  Ergodic Theory and Dynamical Systems, pages 1–44, 2021.
  * [12] S. J. Esnay and E. Moutot.  Aperiodic SFTs on Baumslag-Solitar groups.  Theoretical Computer Science, 2022.
  * [13] E. Jeandel.  Aperiodic subshifts of finite type on groups.  arXiv:1501.06831, 2015.
  * [14] E. Jeandel.  Aperiodic subshifts on polycyclic groups.  arXiv:1510.02360, 2015.
  * [15] E. Jeandel and M. Rao.  An aperiodic set of 11 Wang tiles.  Advances in Combinatorics, pages Paper No. 1, 37, 2021.
  * [16] J. Kari.  A small aperiodic set of Wang tiles.  Discrete Mathematics, 160:259–264, 1996.
  * [17] D. A. Lind and B. Marcus.  An Introduction to Symbolic Dynamics and Coding.  Cambridge University Press, 2 edition, 2021.
  * [18] D. Osin.  Small cancellations over relatively hyperbolic groups and embedding theorems.  Annals of Mathematics. Second Series, 172(1):1–39, 2010.
  * [19] R. M. Robinson.  Undecidability and nonperiodicity for tilings of the plane.  Inventiones Mathematicae, 12:177–209, 1971.
  * [20] H. Wang.  Proving theorems by pattern recognition, II.  In Computation, Logic, Philosophy, pages 159–192. Springer Netherlands, 1961.