# Fast Fourier Transform periodic interpolation method for superposition sums
in a periodic unit cell

Fangzhou Ai  Vitaliy Lomakin Corresponding Author. Email<EMAIL_ADDRESS>
###### Abstract

We propose a Fast Fourier Transform based Periodic Interpolation Method (FFT-
PIM), a flexible and computationally efficient approach for computing the
scalar potential given by a superposition sum in a unit cell of an infinitely
periodic array. Under the same umbrella, FFT-PIM allows computing the
potential for 1D, 2D, and 3D periodicities for dynamic and static problems,
including problems with and without a periodic phase shift. The computational
complexity of the FFT-PIM is of $O(N\log N)$ for $N$ spatially coinciding
sources and observer points. The FFT-PIM uses rapidly converging series
representations of the Green’s function serving as a kernel in the
superposition sum. Based on these representations, the FFT-PIM splits the
potential into its near-zone component, which includes a small number of
images surrounding the unit cell of interest, and far-zone component, which
includes the rest of an infinite number of images. The far-zone component is
evaluated by projecting the non-uniform sources onto a sparse uniform grid,
performing superposition sums on this sparse grid, and interpolating the
potential from the uniform grid to the non-uniform observation points. The
near-zone component is evaluated using an FFT-based method, which is adapted
to efficiently handle non-uniform source-observer distributions within the
periodic unit cell. The FFT-PIM can be used for a broad range of applications,
such as periodic problems involving integral equations in computational
electromagnetic and acoustic, micromagnetic solvers, and density functional
theory solvers.

Department of Electrical and Computer Engineering, University of California
San Diego,

San Diego 92093, U.S.A

Center for Memory and Recording Research, University of California San Diego,

San Diego 92093, U.S.A

##  1 Introduction

A common category of problem in computational physics is calculating the
potential, referred to as periodic scalar potential (PSP), that is generates
by a non-uniform distribution of sources arranged in a unit cell that is
extended to an infinite periodic array (Fig. 1). Examples of such problems are
phased-array antennas, crystals, periodic gratings, periodically modulated
waveguides, to name a few [1, 2, 3, 4]. These problems may present challenges
for obtaining accurate and rapid solutions. For a standalone unit cell without
periodicity, this task can be accomplished via a standard superposition sum
involving free-space Green’s function as its kernel. Such sums can be
evaluated rapidly in $O(N)$ or $O(N\log N)$ computational cost for $N$ sources
and observers using several fast methods, such as fast multipole like method
[5, 6], H-matrix method [7], interpolation-based methods [8, 9, 10], and FFT-
based methods [11, 12, 13, 14]. However, for the periodic case, the
superposition sums involve infinite number of images, which poses additional
challenges to the slow converge or even divergence as well as a potentially
high computational cost.

Figure 1: Illustration for a 2D periodic problem consisting of an infinite 2D
array of cubes. The central red cube is the zeroth unit cell. Surrounding it
along the $x$ and $y$ axes are the $1^{st}$ order image cubes in green and the
$2^{nd}$ order image cubes in blue.

Several existing methods can address such periodic problems with certain
assumptions and limitations. Periodic Green’s function (PGF) can be defined in
terms of an infinite sum over periodic images. In dynamic problems, e.g., in
computational electromagnetics, PGF is assumed to be defined with a phase
shift or non-vanishing wavenumber [15, 16, 17]. However, for static problems
in which no phase shift and wave propagation are present, even in the 1D case,
the PGF diverges. In molecular dynamics problems [18], a special form of
static PGF assumes uniformly arranged sources [19], such as a cubic lattice.
Although various methods have been proposed to address static and dynamic
periodic unit cell problems with non-uniform sources, they often focus on
specific problem types, such as addressing only 1D or 2D periodicity,
addressing only static problems [20], requiring lattice structure or relying
on finite difference methods [21]. There is a need in a general approach
capable of effectively handling all types of these problems regardless of the
periodicity type, distribution of sources, with or without phase shift, and
for static and dynamic cases.

In this paper, we proposed a general approach, referred to as Fast Fourier
Transform Periodic Interpolation Method (FFT-PIM) that is applicable to all
these types of problems under a unified framework, including 1D, 2D, and 3D
periodicities for static and dynamic problems with and without periodic phase
shifts. The FFT-PIM separates the PSP into a near-zone component involving a
small number of images near the periodic unit cell of interest, and a
spatially slowly varying far-zone component involving the rest of an infinite
number of images. The evaluation of the near-zone PSP component is based on a
modified version of the box adaptive integral method [11]. The evaluation of
the far-zone PSP component is based on the sparse periodic interpolation
method that is adapted to allow handling dynamic and static problems. The
computational cost of the FFT-PIM is of $O(N\log N)$ and the memory
consumption is of $O(N)$.

The paper is organized as follows. Sec. 2 presents the problem formulation.
Sec. 3 presents the algorithmic foundations of the FFT-PIM, including near-
and far-zone representations of PGF and PSP as well as efficient ways of their
evaluation. Sec. 4 shows results for the use of the FFT-PIM, its error
analysis, and computational performance. Finally, Sec. 5 presents summary and
conclusions.

##  2 Formulation

Consider an infinite periodic array of unit cells residing in free space (see
Fig. 1 for a 2D periodicty example). Within each unit cell there are $N$
source points located at $\mathbf{r}_{n}$. In the zeroth unit cell, the source
amplitudes are $q_{n}$, where $n=1...N$. In each unit cell, there is the same
number of coinciding observer points at which the periodic scalar potentials
(PSPs) $u(\mathbf{r}_{i})$ are to be found. The periodic array can be 1D, 2D
or 3D with the periodicity of $L_{x}$, $L_{y}$, and $L_{z}$ in three possible
directions $x$, $y$, and $z$, respectively. The amplitudes of the sources may
be periodically phase shifted with a linear, possibly complex, phase shift
determined by the wavenumbers $k_{x0}$, $k_{y0}$, and $k_{z0}$ along the $x$,
$y$, and $z$ directions. A free space wavenumber $k_{0}$ describes the
propagation of the waves in the free space for the dynamic case. For
$k_{0}=0$, the problem is static, and for $k_{x0}=k_{y0}=k_{z0}=0$, the
sources are not phase shifted.

The PSP in the zeroth unit cell can be written as

$u(\mathbf{r}_{m})=\sum_{n=1}^{N}G^{p}(r_{m}-r_{n})q(\mathbf{r}_{n}),$ (1)

where $m=1...N$ and $G^{p}$ represents the scalar periodic Green’s function
(PGF), given for the 1D, 2D, and 3D cases as:

$\displaystyle G^{p}_{1D}(\mathbf{r})=$
$\displaystyle\sum_{i_{x}=-\infty}^{\infty}e^{-jk_{x0}(i_{x}L_{x})}G_{0}(\mathbf{r}-i_{x}L_{x}\hat{x}),$
(2a) $\displaystyle G^{p}_{2D}(\mathbf{r})=$
$\displaystyle\sum_{i_{x}=-\infty}^{\infty}\sum_{i_{y}=-\infty}^{\infty}e^{-j[k_{x0}(i_{x}L_{x})+k_{y0}(i_{y}L_{y})]}G_{0}(\mathbf{r}-i_{x}L_{x}\hat{x}-i_{y}L_{y}\hat{y}),$
(2b) $\displaystyle G^{p}_{3D}(\mathbf{r})=$
$\displaystyle\sum_{i_{x}=-\infty}^{\infty}\sum_{i_{y}=-\infty}^{\infty}\sum_{i_{z}=-\infty}^{\infty}e^{-j[k_{x0}(i_{x}L_{x})+k_{y0}(i_{y}L_{y})+k_{z0}(i_{z}L_{z})]}G_{0}(\mathbf{r}-i_{x}L_{x}\hat{x}-i_{y}L_{y}\hat{y}-i_{z}L_{z}\hat{z}),$
(2c)

respectively. Here, $G_{0}$ is the free space scalar Green’s function given by
$G_{0}(\mathbf{r})=\exp(-jk_{0}|\mathbf{r}|)/4\pi|\mathbf{r}|$.

The series in Eq. (1) for the PSP is important in multiple computational
physics problems. For example, the evaluation of the superposition integrals
appearing in electromagnetic integral equation or micromagnetic solvers [22]
can be represented as a product of sparse matrices describing the mesh of the
geometry and the result of the superposition sum of Eq. (1). The sum of Eq.
(1) can also be directly apply to a discrete set of sources and observers,
e.g., when considering a set of point sources.

The series of Eq. (2) for the PGF is slowly convergent, which is because of a
slow spatial decay of $G_{0}$. An alternative spectral representation of the
PGF of Eq. (2) can be based on the Floquet mode expansion [17]:

$\displaystyle G^{p}_{1D}(\mathbf{r})=$
$\displaystyle\sum_{m=-\infty}^{\infty}\frac{1}{4jL_{x}}e^{-jk_{xm}x}H_{0}^{(2)}(k_{\rho
m}\sqrt{y^{2}+z^{2}}),$ (3a) $\displaystyle G^{p}_{2D}(\mathbf{r})=$
$\displaystyle\sum_{m,n=-\infty}^{\infty}\frac{e^{-jk_{xm}x-jk_{yn}y-jk_{zmn}|z|}}{2jk_{zmn}L_{x}L_{y}},$
(3b) $\displaystyle G^{p}_{3D}(\mathbf{r})=$
$\displaystyle\sum_{m,n=-\infty}^{\infty}\frac{e^{-jk_{xm}x-jk_{yn}y}}{2jk_{zmn}L_{x}L_{y}}$
$\displaystyle\times\left(\frac{e^{-jk_{zmn}|z|}}{1-e^{-j(k_{zmn}-k_{z0})L_{z}}}-\frac{jk_{zmn}|z|}{1-e^{-j(k_{zmn}+k_{z0})L_{z}}}\right),$
(3c)

where $H_{0}^{(2)}$ is the zeroth order Hankel function of the second kind,
$k_{xm}=k_{x0}+2\pi m/L_{x}$, $k_{ym}=k_{y0}+2\pi n/L_{z}$, $k_{\rho
m}=(k_{0}^{2}-k_{xm}^{2})^{1/2}$ and
$k_{zmn}=(k_{0}^{2}-k_{xm}^{2}-k_{yn}^{2})^{1/2}$, and the square roots are
chosen such that their imaginary parts are non-positive. The spectral series
of Eq. (3) converges exponentially fast when $(y^{2}+z^{2})^{1/2}$ is non-
vanishing and $k_{0},k_{x0},k_{y0},k_{z0}$ are not all vanishing.

A special attention needs to be paid to the case of
$k_{0}=k_{x0}=k_{y0}=k_{z0}=0$, which we refer to as static case. The static
case is important as it corresponds to several practical problem, such as
using periodic extensions to represent infinite 1D, 2D, and 3D domains in
various static problems, e.g., to characterize micromagnetic [23] or in
molecular dynamics problems. In this static periodic case, the PGF calculated
via series of Eq. (2) or Eq. (3) diverges and the PSP for a general source
distribution diverges as well. However, a finite PSP can be obtained for the
special case of a neutral source, i.e., under the condition

$\sum_{n=1}^{N}q_{n}=0.$ (4)

This neutrality condition corresponds to various physical problems, e.g., when
computing the magnetic scalar potential generated by a magnetized objects or
electrostatic potential in polarized objects [24]. In this case, the PSP
converges and we can neglect any non-relevant components (e.g., constant
terms) in the PGF, so that it is convergent as well. The result is that the
PGF for the neutral static case can be calculated via:

$\displaystyle G^{p}_{1D}(\mathbf{r})=$ $\displaystyle-\frac{1}{2\pi
L_{x}}\ln\left(\sqrt{y^{2}+z^{2}}\right)+\frac{1}{\pi
L_{x}}\sum_{m=1}^{\infty}K_{0}\left(\frac{2m\pi}{L_{x}}\sqrt{y^{2}+z^{2}}\right)\cos\left(\frac{2m\pi
x}{L_{x}}\right)$ (5a) $\displaystyle G^{p}_{2D}(\mathbf{r})=$
$\displaystyle-\frac{|z|}{2L_{x}L_{y}}-\frac{1}{4\pi
L_{x}}\ln\left(1-2e^{-\frac{2\pi|z|}{L_{y}}}\cos\left(\frac{2\pi
y}{L_{y}}\right)+e^{-\frac{4\pi|z|}{L_{y}}}\right)$
$\displaystyle+\frac{1}{\pi
L_{x}}\sum_{m=1}^{\infty}\sum_{n=-\infty}^{\infty}K_{0}\left(\frac{2m\pi}{L_{x}}\sqrt{(nL_{y}+y)^{2}+z^{2}}\right)\cos\left(\frac{2m\pi
x}{L_{x}}\right)$ (5b) $\displaystyle G^{p}_{3D}(\mathbf{r})=$
$\displaystyle\frac{1}{2L_{x}L_{y}L_{z}}(z^{2}-|z|L_{z})$
$\displaystyle-\frac{1}{4\pi
L_{x}}\sum_{k=-\infty}^{\infty}\ln\left(1-2e^{-\frac{2\pi|kL_{z}+z|}{L_{y}}}\cos\left(\frac{2\pi
y}{L_{y}}\right)+e^{-\frac{4\pi|kL_{z}+z|}{L_{y}}}\right)$
$\displaystyle+\frac{1}{\pi
L_{x}}\sum_{m=1}^{\infty}\sum_{n,k=-\infty}^{\infty}K_{0}\left(\frac{2m\pi}{L_{x}}\sqrt{(nL_{y}+y)^{2}+(kL_{z}+z)^{2}}\right)\cos\left(\frac{2m\pi
x}{L_{x}}\right)$ (5c)

where $K_{0}$ is the zeroth order modified Bessel function of the second kind.

It can be shown that Eq. (5a) is obtained from Eq. (3a) under the assumption
of Eq. (4). In this case, by taking the limit of ($k_{0}\to 0,k_{x0}\to 0$)
for all $n\neq 0$, the Hankel function becomes the modified Bessel function
and the symmetric $\pm n$ terms can be combined into the positive $n$ terms.
For the $n=0$ term, the Hankel function exhibits a logarithmic behavior but it
is summed up to $0$ due to the source neutrality of Eq. (4). The PGF for the
2D and 3D cases can be then obtained by a spatial sum of 1D PGFs. When the
sources are arranged on a lattice, the PGF for the periodic static case
sometimes is also referred to as the lattice Green’s function [19].

The PGFs converge exponentially fast provided the $(y^{2}+z^{2})^{1/2}$ is not
too small. For a desired error of $\epsilon$, the number of terms required Eq.
(3a) and Eq. (5a) scales as $L_{x}\log(\epsilon^{-1})/z$ for the 1D case and
$L_{x}L_{y}\log^{2}(\epsilon^{-1})/z^{2}$ for the 2D and 3D cases. To
demonstrate the convergence, Fig. 2 shows the PGFs and their convergence when
using Eq. (3) and Eq. (5) for $L_{x},L_{y},L_{z}=1$, $k_{0}=-1-j$ and
$k_{x0}=1-j,k_{y0}=1+j,k_{z0}=-1+j$ in Fig. 2 and the static case in Fig. 2.

Figure 2: Convergence of PGFs, $L_{x},L_{y},L_{z}=1$, $x=y=z=0.5$. (a) The
relative error of the sum of the first $m$ terms from Eq. (3), with
$k_{0}=-1-j$ and $k_{x0}=1-j,k_{y0}=1+j,k_{z0}=-1+j$, (b) The relative error
of the sum of the first $m$ terms from Eq. (5). Exponential convergence is
schieved in both cases.

##  3 Algorithm

###  3.1 Representation of PGF and PSP

We now describe the FFT-PIM that allows rapidly computing the PSP in Eq. (1).
We start by representing the PGF in terms of its near-zone component
$G_{near}^{p}$ and far-zone component $G^{p}_{far}$:

$G^{p}=G_{near}^{p}+G^{p}_{far}.$ (6)

Here, the near-zone component $G_{near}^{p}$ is given by Eq. (2) but with the
sums over $i_{x}$, $i_{y}$, and $i_{z}$ in the range of $[-i_{xd},i_{xd}]$,
$[-i_{yd},i_{yd}]$, and $[-i_{zd},i_{zd}]$, respectively, where
$i_{xd},i_{yd},i_{zd}$ are small integer numbers. This components represent
the contribution from a small number of images around the zeroth unit cell of
interest. Assuming that the periodicities are close to each other, we can use
the same value of $i_{d}$ for $i_{xd},i_{yd},i_{yd}$. For the sake of clarity
we assume this case in the following discussions. $G_{near}^{p}$ includes the
self-term and a few terms from the surrounding images (see Fig. 1 showing the
near-zone images for the case of $i_{xd}=i_{yd}=1$ as green boxes). Therefore,
$G_{near}^{p}$ may be spatially singular and rapidly varying. The evaluation
of $G_{near}^{p}$ involves only a small number of the sum terms.

The far-zone component $G_{far}^{p}$ is the rest of the images from far-away
unit cells and can be evaluated as $G^{p}_{far}=G^{p}-G_{near}^{p}$ (blue
boxes in Fig. 1). The far-zone PGF $G_{far}^{p}$ is smoothly varying in space
because it corresponds to large source-observer spatial separations. The rate
of the variations is reduced by using a greater $i_{d}$.

Following the near- and far-zone decomposition of the PGF, the potential is
also expressed in terms of its near- and far-zone components:

$\begin{split}u(\mathbf{r}_{m})&=u_{near}(\mathbf{r}_{m})+u_{far}(\mathbf{r}_{m}),\\\
u_{near}(\mathbf{r}_{m})&=\sum_{n=1,n\neq
m}^{N}G_{near}^{p}(\mathbf{r}_{m}-\mathbf{r}_{n})q_{n},\\\
u_{far}(\mathbf{r}_{m})&=\sum_{n=1,n\neq
m}^{N}G^{p}_{far}(\mathbf{r}_{m}-\mathbf{r}_{n})q_{n},\end{split}$ (7)

where, again, the far-zone potential $u_{far}$ is a spatially slowly varying
function. These potential and PGF decompositions lead to a fast approach for
the numerical evaluation of the PSP $u(\mathbf{r}_{m})$.

###  3.2 Evaluation of the far-zone PSP

The far-zone PGF and PSP components vary slowly in space, which enables us to
calculate the far-zone PGF and PSP on a uniform grid within the observation
domain and subsequently interpolate to the required observer points. We
construct a sparse uniform grid of observation points
{$\mathbf{r}_{n(l,p,q)}^{o}={x_{l}^{o},y_{p}^{o},z_{q}^{o}}$} (red grid in
Fig. 3) and a spatially shifted uniform grid of source points
{$\mathbf{r}_{n(l,p,q)}^{s}={x_{l}^{s},y_{p}^{s},z_{q}^{s}}$} (black grid in
Fig. 3):

$\displaystyle
x_{n(l,p,q)}^{s}=\frac{(l-1)(D_{x}-\Delta_{x})}{N_{gx}-1},y_{n(l,p,q)}^{s}=\frac{(p-1)(D_{y}-\Delta_{y})}{N_{gy}-1},z_{n(l,p,q)}^{s}=\frac{(q-1)(D_{z}-\Delta_{z})}{N_{gz}-1},$
(8a) $\displaystyle x_{n(l,p,q)}^{o}=x_{n(l,p,q)}^{s}+\frac{\Delta_{x}}{2},\ \
\ \ \ y_{n(l,p,q)}^{o}=y_{n(l,p,q)}^{s}+\frac{\Delta_{y}}{2},\ \ \ \ \ \ \
z_{n(l,p,q)}^{o}=z_{n(l,p,q)}^{s}+\frac{\Delta_{z}}{2},$ (8b)
$\displaystyle\Delta_{x}=\frac{D_{x}}{N_{gx}},\
\Delta_{y}=\frac{D_{y}}{N_{gy}},\ \Delta_{z}=\frac{D_{z}}{N_{gz}};\
n=lN_{gz}N_{gy}+pN_{gz}+q.$ (8c)

Here, $N_{gx}$, $N_{gy}$, $N_{gz}$ are the number of the grid points in the
$x$, $y$, $z$ dimensions, $N_{g}=N_{gx}N_{gy}N_{gz}$ is the total number of
the grid points, which is of $O(1)$, with $n\in[1,N_{g}]$, and $\Delta_{x}$,
$\Delta_{y}$, $\Delta_{z}$ are grid spacing on in the $x$, $y$, $z$
dimensions. The source grid is shifted from the observer grid by
$\Delta_{x}/2$, $\Delta_{y}/2$, $\Delta_{z}/2$ to make the grids non-
overlapping, which results in a fast convergence of the sum in Eq. (3) and Eq.
(5) when calculating PGF.

Based on this grid, the PGF at any observation point $\mathbf{r}_{m}$ from any
source point $\mathbf{r}_{n}$ in the zeroth unit cell can be calculated by
interpolation:

$G^{p}_{far}(\mathbf{r}_{m}-\mathbf{r}_{n})=\sum_{m^{\prime}=1}^{N_{go}(\mathbf{r}_{m})}\sum_{n^{\prime}=1}^{N_{gs}(\mathbf{r}_{n})}\omega^{o}(\mathbf{r}_{m},\mathbf{r}_{m^{\prime}}^{o})G^{p}_{far}(\mathbf{r}_{m^{\prime}}^{o}-\mathbf{r}_{n^{\prime}}^{s})\omega^{s}(\mathbf{r}_{n^{\prime}}^{s},\mathbf{r}_{n}).$
(9)

Here, $\omega^{s}(\mathbf{r}_{n^{\prime}}^{s},\mathbf{r}_{n})$ represents the
interpolation coefficients from the source grid points
$\mathbf{r}_{n^{\prime}}^{s}$ to original non-uniform source points
$\mathbf{r}_{n}$ and $\omega^{o}(\mathbf{r}_{m},\mathbf{r}_{m^{\prime}}^{o})$
represents the interpolation coefficients from the observer grid points
$\mathbf{r}_{m^{\prime}}^{o}$ to the original non-uniform observer points
$\mathbf{r}_{m}$. The number of the source interpolation points
$N_{gs}(\mathbf{r}_{n})$ and observer interpolation points
$N_{go}(\mathbf{r}_{n})$ represents the number of the grid points that need to
be used for the interpolation. For example, we can choose the $q$-th order
Lagrange interpolation such that all grid points participate in the
interpolation for any of the non-uniform points. For this choice,
$N_{q}=(q+1)^{3}$ , $N_{go}=N_{gs}=N_{g}$, and the coefficients $\omega^{o}$
and $\omega^{s}$ are the same for the same points. For a more general case of
a grid with $N_{q}$ greater than $N_{go}$ and $N_{gs}$ , the interpolation for
each non-uniform point involve only a subset of the grid points, and the
coefficients $\omega^{o}$ and $\omega^{s}$ can be represented as sparse
matrices. The corresponding sparse matrices are generally transpose versions
of each other.

We can substitute $G^{p}_{far}(\mathbf{r}_{m}-\mathbf{r}_{n})$ from Eq. (9)
into Eq. (1) and re-arrange the sums to obtain an alternative representation:

$u(\mathbf{r}_{m})=\sum_{m^{\prime}=1}^{N_{go}(\mathbf{r}_{m})}\omega^{o}(\mathbf{r}_{m},\mathbf{r}_{m^{\prime}}^{o})\sum_{n^{\prime}=1}^{N_{g}}G^{p}_{far}(\mathbf{r}_{m^{\prime}}^{o}-\mathbf{r}_{n^{\prime}}^{s})\sum_{n=1}^{N_{s}(\mathbf{r}_{n^{\prime}}^{s})}\omega^{s}(\mathbf{r}_{n^{\prime}}^{s},\mathbf{r}_{n})q(\mathbf{r}_{n}).$
(10)

Here, $N_{s}$ is the number of the non-uniform source points that participate
in the process of interpolations with a source grid point
$\mathbf{r}_{n^{\prime}}^{s}$. These grid points can be found from the
interpolation procedure of Eq. (9). For example, when choosing the $q$-th
order Lagrange interpolation with $N_{q}=(q+1)^{3}$, $N_{s}=N$. For a more
general choice of a greater $N_{q}$, $N_{s}$ is found from the sparse matrix
representation of the interpolation coefficients $\omega^{o}$ and $\omega^{s}$
.

Based on this grid construction in Eq. (8) and PSP representation in Eq. (10),
we can first pre-compute the coefficients $\omega^{o}$ and $\omega^{s}$, e.g.,
as a sparse matrix, and pre-compute
$G^{p}_{far}(\mathbf{r}_{m^{\prime}}^{o}-\mathbf{r}_{n^{\prime}}^{s})$ table
at the grid points. The computation of $\omega^{o}$ and $\omega^{s}$ involves
$O(N)$ operations since the interpolations are spatially local. The
computation of
$G^{p}_{far}(\mathbf{r}_{m^{\prime}}^{o}-\mathbf{r}_{n^{\prime}}^{s})$
involves $O(N_{g})$ operations since the grid is uniform and $G^{p}_{far}$ is
shift invariant in term of the differences
$\mathbf{r}_{m^{\prime}}^{o}-\mathbf{r}_{n^{\prime}}^{s}$. The computation of
$G^{p}_{far}(\mathbf{r}_{m^{\prime}}^{o}-\mathbf{r}_{n^{\prime}}^{s})$ can be
done via the rapidly convergent sums of Eq. (3) and Eq. (5) due to the choice
of the source and observer grid in Eq. (8) with a spatial separation. The far-
zone PSP $u_{far}$ can, then, be evaluated via the following three steps,
which are illustrated in Fig. 3.

Step 1, Projection from the source points to the source grid: This step
calculates effective charges at the uniform source grid. Substituting the
representation of Eq. (10) for $G^{p}_{far}$ into Eq. (1) for the far-zone PGF
and potential components allows obtaining an expression for the effective
charge at the source grid points as projections from the original non-uniform
sources (black arrows in Fig. 3 represent a contribution of a single non-
uniform source to its relevant source grid points):

$q_{g}(\mathbf{r}_{n^{\prime}}^{s})=\sum_{n=1}^{N_{s}(\mathbf{r}_{n^{\prime}}^{s})}\omega^{s}(\mathbf{r}_{n^{\prime}}^{s},\mathbf{r}_{n})q(\mathbf{r}_{n}),$
(11)

where $n^{\prime}\in[1,N_{s}]$.

Step 2, Evaluation of the PSP at the observer grid: This step calculates the
far-zone PSP at the sparse observer grid points from the sparse source
effective charges via the following convolution (green arrows in Fig. 3):

$u_{far}^{g}({\mathbf{r}_{m^{\prime}}^{o}})=\sum_{n^{\prime}=1}^{N_{g}}G^{p}_{far}(\mathbf{r}_{m^{\prime}}^{o}-\mathbf{r}_{n^{\prime}}^{s})q_{g}(\mathbf{r}_{n^{\prime}}^{s}),$
(12)

which can be evaluated rapidly directly since $N_{g}=O(1)$.

Step 3, Interpolation from the observer grid to the observer points: The far-
zone component of the PSP at the original non-uniform points can be finally
calculated by interpolation, which can be viewed as an inverse procedure of
step 1 (red arrows in Fig. 3):

$u_{far}(\mathbf{r}_{m})=\sum_{m^{\prime}=1}^{N_{go}}\omega^{o}(\mathbf{r}_{m},\mathbf{r}_{m^{\prime}}^{o})u_{far}^{g}(\mathbf{r}_{m^{\prime}}^{o}),$
(13)

where $m\in[1,N]$.

Figure 3: Illustration of the source (black) and observer (red) grids for the
far-zone PSP component. Black and red circles are source/observer grid points
on the Cartesian lattices. Black and red dots are source/observer points with
arbitrary distribution. Black arrows denote projections from source points to
source grid points. Green arrow is the grid-to-grid interaction, performed by
convolution. Red arrows represent interpolations from observer grid points to
observer points.

###  3.3 Evaluation of the near-zone PSP

To evaluate the near-zone PSP component, $u_{near}$, we note that for
$i_{d}=0$, we only consider the zeroth order unit cell, which is equivalent to
a non-periodic problem. The evaluation of $u_{near}$ for a free-space
standalone problem can be efficiently handled using various fast computational
approaches, such as the fast multipole like methods [5, 6] or interpolation-
based approaches [8, 9, 10]. For $i_{d}>0$, neighbor images in the near-zone
evaluation are included and they require modifying the fast evaluation
algorithms. Here, we use the box-adaptive integral method of [11, 25], which
is also related to the pre-corrected FFT method [14] and adaptive integral
method [12, 13]. We modify this approach to efficiently account for the
inclusion of the periodic images. The idea is similar to the far-zone
procedure in Sec. 3.2 in its first three steps, but it also adds another
modified correction step.

Similar to the grid construction in Sec. 3.2, we construct uniform grids.
However, there are two notable differences. The first difference is that the
source and observer uniform grids coincide, i.e., there is a single uniform
grid {$\mathbf{r}_{m}={x_{m},y_{m},z_{m}}$} =
{$\mathbf{r}_{m}^{o}={x_{m}^{o},y_{m}^{o},z_{m}^{o}}$} =
{$\mathbf{r}_{n}^{s}={x_{n}^{s},y_{n}^{s},z_{n}^{s}}$} by removing
$\Delta_{x}$, $\Delta_{y}$, $\Delta_{z}$ in Eq. (8a) and Eq. (8b). The reason
for using a single grid is that the near-zone PGF only involves a small number
of terms in its evaluation and does not need to have a spatial separation
between source and observer points. The second difference is that the number
of the grid points $N_{g}$ is significantly larger than what is used in the
far-zone evaluation. Typically, $N_{g}$ is comparable to the number of the
non-uniform points $N$ for the optimal performance. This large $N_{g}$ is
required to have $N_{s}$ (as defined in Eq. (10)) to be of $O(1)$ to result in
a reduced number of computations in the correction step 4. Based on this grid
definition, the near-zone evaluation proceeds in the following four steps.

Step 1, Projection from the source points to the source grid: This step is
similar to step 1 in Sec. 3.2, where $G^{p}_{far}$ is replaced with
$G^{p}_{near}$. Here, $N_{g}$ is typically much larger than $N_{s}$, and the
matrices of the interpolation coefficients are sparse. The resulting
computational cost is of $O(N_{g})=O(N)$.

Step 2, Evaluation of the PSP at the observer grid: This step is similar to
step 2 in Sec. 3.2, where $G^{p}_{far}$ is replaced with $G^{p}_{near}$.
However, because $N_{g}$ is large, the convolution sum of Eq. (12) is computed
via FFT, which is allowed because the grids are uniform. Since $G^{p}_{near}$
is not periodic, using FFT requires making the equivalent matrix circular,
which involves doubling its size in each periodicity dimension. For performing
FFTs in the computation stage, $G_{near}^{p}$ is tabulated in the pre-
processing stage. The tabulation computational cost is of $O(N_{g})=O(N)$ and
the FFT computational cost is of $O(N\log N)$

Step 3, Interpolation from the observer grid to the observer points: The
interpolation step for $u_{near}$ is identical to step 3 in Sec. 3.2 for
$u_{far}$ and it can be viewed as an inverse of step 1 using the same sparse
matrix for the interpolation coefficients. The resulting computational cost is
of $O(N_{g})=O(N)$.

Step 4: Error correction: The ability to perform the projection/interpolation
in steps 1 and 3 rely on the assumption of slow variations of the Green’s
function. For the far-zone evaluation, $G^{p}_{far}$ indeed varies slowly in
the domain of the zeroth unit cell because it excludes the direct and near-
zone image interactions between the sources and observers. The near-zone
component with $G^{p}_{near}$, however, involves direct and near-image
interactions, which result in rapid spatial variations, e.g., the rate of such
spatial variations is unbounded when the source and observer coincide.
Therefore, the interpolations for the nearby source-observer pair
contributions are highly inaccurate, and they need to be corrected. To
eliminate these errors, we define an error correction range $D^{ER}$ for each
box, determined by the uniform grid. This range is of the same order as
$max\\{\Delta_{x},\Delta_{y},\Delta_{z}\\}$. For all observer points in a box,
the corrections are performed in the error-correction region,
$\Omega^{ER}(\mathbf{r}_{m})$, for all the non-uniform sources in the same box
and a certain number of surrounding boxes that are withing the $D^{ER}$
distance from the box. The correction involves subtracting the contributions
due to steps 1-3 and adding the exact point-to-point superposition sums via
appropriate contributions in Eq. (1). The error correction procedure can be
expressed as

$\displaystyle
u_{near}(\mathbf{r}_{m})=u_{near}(\mathbf{r}_{m})+\sum_{n}^{\mathbf{r}_{n}\in\Omega^{ER}(\mathbf{r}_{m})}[G^{p}_{near}(\mathbf{r}_{m}-\mathbf{r}_{n})-G_{grid}(\mathbf{r}_{m}-\mathbf{r}_{n})]q(\mathbf{r}_{n}),$
(14a) $\displaystyle
G_{grid}(\mathbf{r}_{m},\mathbf{r}_{n})=\sum_{m^{\prime}}^{\omega^{o}(\mathbf{r}_{m},\mathbf{r}_{m^{\prime}})\neq
0}\>\sum_{n^{\prime}}^{\omega^{s}(\mathbf{r}_{n^{\prime}},\mathbf{r}_{n})\neq
0}\omega^{o}_{nufft}(\mathbf{r}_{m},\mathbf{r}_{m^{\prime}})G^{p}_{near}(\mathbf{r}_{m^{\prime}}-\mathbf{r}_{n^{\prime}})\omega^{s}_{nufft}(\mathbf{r}_{n^{\prime}},\mathbf{r}_{n}),$
(14b)

where the error is corrected in the region $\Omega^{ER}(\mathbf{r}_{m})$,
i.e., for any point $\mathbf{r}$ within this region,
$|\mathbf{r}-\mathbf{r}_{m}|\leq D^{ER}$. Similar to Eq. (11) and Eq. (13),
the $\omega^{s}_{nufft}$ and $\omega^{o}_{nufft}$ correspond to the
interpolation matrices containing the weights for the above step 1 and 2. The
function $G_{grid}$ represents the near-zone Green’s function between the
source and observer points obtained via the grid interactions with
interpolation, which needs to be subtracted in the error correction process.
While the procedure in Eq. (14) can be efficient, it still can significantly
increase the computational cost due to the need to account for all the $i_{d}$
near-zone images via $G^{p}_{near}$. For example, in the 3D periodicity case,
$i_{d}=125$ when including 2 surrounding boxes per each dimension, which leads
to a high additional computational cost.

This computational cost can be much reduced by including only a single extra
periodic image per periodicity dimension. This is possible due to the fact
that $\\{\Delta_{x},\Delta_{y},\Delta_{z}\\}$ are typically much smaller than
$\\{D_{x},D_{y},D_{z}\\}$, so that most of the images correspond to a large
spatial separation. For example, in Fig. 4 showing an example for a 1D
$x$-dimension periodicity case, for the sources in the internal boxes of the
zeroth (middle) unit cell, only the $G_{0}(\mathbf{r}_{m}-\mathbf{r}_{n})$
component of the $G^{p}_{near}(\mathbf{r}_{m}-\mathbf{r}_{n})$ can be used for
the near-zone correction. On the other hand, for the left box, which is next
to the left edge of the zeroth unit cell (the left green box with red dot and
arrows), in addition to the $G_{0}(\mathbf{r}_{m}-\mathbf{r}_{n})$ component,
the component of $G_{0}(\mathbf{r}_{m}-\mathbf{r}_{n}-L_{x}\hat{x})$ of
$G^{p}_{near}(\mathbf{r}_{m}-\mathbf{r}_{n})$ needs to be accounted for. The
latter component corresponds to the sources from the right box (right green
box with a black dot and arrows) contributing to the left box observers via
the $i_{x}=1$ image. These sources are equivalent to the sources in the right
box of the left image cell (the dashed green box with a grey dot) contributing
to the left box observers via the image Green’s function
$G_{0}(\mathbf{r}_{m}-\mathbf{r}_{n})$, i.e., they are equivalent to the
geometrically close boxes. This procedure is generalized to 2D and 3D
periodicities by including corresponding required images. The computational
cost of this procedure is only slightly higher than that for a non-periodic
case, in which only the $G_{0}(\mathbf{r}_{m}-\mathbf{r}_{n})$ contribution is
required. The cost increase is associated with the need to include additional
images at the sources associated with the boxes located at the corners, edges,
and sides of the computational domain. There is a small number of such boxes
as compared to the total number of near-zone images, so that the computational
cost increase is insignificant. The overall computational cost of this step is
of $O(N)$.

Figure 4: Illustration for the error correction step 4 for the near-zone PSP
component. The black dot is a source point and red dot is an observer point.
The green region is the error-correction region $\Omega^{ER}$ . The black
circles are grid points. For the observers in $\Omega^{ER}$, PSPs are
calculated through direct calculations, this is done by subtracting the grid
interaction inside the green region (dotted arrows) and adding the direct
interactions (dashed purple arrow). The shaded green region and grey dot is
the image of the right boundary in the $i_{d}=-1$ image cell, which is close
to the left boundary of the unit cell within $\Omega^{ER}$.

##  4 Numerical results

We implemented FFT-PIM both in a Central Processing Unit (CPU) and Graphics
Processing Unit (GPU) based code. This section presents numerical test results
demonstrating the accuracy and performance of the FFT-PIM. Sec. 4.1 presents
error analysis, comparing the results from the FFT-PIM with results obtained
via the direct sum of Eq. (1) with Eq. (5a) for different problem sizes,
uniform grid sizes, and interpolation orders. Sec. 4.2 examines the
computational speed of the FFT-PIM for different parameters, such as problem
size, uniform grid size, and the periodicity dimension.

We start by showing results for the calculated potential due to a source
distribution (Fig. 5). The structure is 1D periodic in the $x$-direction with
the periodicity of $L_{x}=1$ and the unit cell comprises a section of coaxial
cable along the $x$-axis that extends through the entire period. The inner and
outer radius is $r_{1}=1$ and $r_{2}=2$, respectively. Negative and positive
sources (e.g., charges) are distributed uniformly on the inner surface with
charge density $\rho_{1}=-1$ and $\rho_{2}=1/2$, respectively, resulting in a
neutral unit cell. We first calculate the potential on the central $x$ axis
for a non-periodic case, as shown in Fig. 5 and we observe significant edge
effects in terms of the potential spatial dependence. We then use the FFT-PIM
with a 1D periodicity, making it into an infinitely long coaxial source
distribution. With the static condition and no phase shift, we find that the
potential becomes near zero everywhere (blue curve in Fig. 5), the edge effect
is eliminated, which agrees with the analytical result. We also calculate the
potential for the dynamic case with the wavenumber $k_{0}=1$ and a complex
phase shift $k_{x0}=1-j$, in which case the potential becomes complex with
non-zero magnitude and corresponding real/imaginary parts (black, red, and
green curves in Fig. 5).

Figure 5: Potential on the $x$ axis of a coaxial structure . (a) Illustration
of coaxial structure unit cell ($L=1$) and potential of the standalone unit
cell. Inner radius $r_{1}=1$ with negative line charge density $\rho_{1}=-1$,
outer radius $r_{2}=2$ with positive line charge density $\rho_{2}=1/2$. (b)
PSP with a 1D periodicity along the $x$ axis with $L_{x}$ = 1. The blue line
is the magnitude of PSP for the static case with $k_{0}=k_{x0}=0$. The red,
green and black lines are real part, imaginary part, and magnitude of PSP for
the dynamic case with $k_{0}=1$, $k_{x0}=1-j$.

We, then, proceed with showing results for the error analysis and
computational performance.

###  4.1 Error analysis

To test the error of the FFT-PIM, we use an infinite long 1D bar along the $x$
axis, and the target object is a cube of the size $D_{x},D_{y},D_{z}=50$ with
1D $x$-direction periodicity with $L_{x}=D_{x}=50$. We assume a static case
without a phase shift, i.e., $k=k_{x0}=0$. We mesh the domain with a
tetrahedral mesh and the total number of the mesh points (vertices) is
$N=7189$. In order to calculate the error, we randomly set 7000 points with
finite values as source points while the rest of the sources are set to $0$,
and make sure that $\sum_{n=1}^{N}q_{n}=0$. The rest of the points are
considered as observer points, and the location of the observer points are
chosen such that they are spatially shifted as compared to the source points
such that the PGF in Eq. (5a) converges for all source-observer pairs when the
direct sum is used. We compare the results for errors between the FFT-PIM and
the direct computation via Eq. (5a), which is equivalent to Eq. (3a) as shown
in Sec. 2. We tested the errors at the observer points with different
interpolation orders, uniform grid size, and number of subtracted near-zone
images. The results are shown in Fig. 6. We find that with cubic interpolation
order and a uniform grid of size $N_{gx},N_{gy},N_{gz}=10$, by subtracting the
zeroth order and only the nearest neighbors ($i_{d}=1$), we can achieve an
error at the level of $10^{-3}$, which is sufficient for many practical cases.
The error is further reduced by increasing the number of grid points, order,
and the near-zone images.

Figure 6: Error of the far-zone PSP for a 1D periodicity along the $x$-axis.
(a) Error for 1st, 3rd, and 6th order interpolation with varying uniform grid
sizes with $i_{d}=1$. (b) Error of cubic interpolation with 1000 sparse grid
points and different numbers of subtracted near-zone unit cells.

###  4.2 Computational performance analysis

Figure 7: Preprocessing time for the no-periodicity case, 1D periodicity case,
and for the difference between the two, i.e., for the overhead due to the
periodicity.

Figure 8: Execution time of evaluating the potential for the standalone case
and for PSP. (a) Execution time on a single-core CPU with (red line) and
without (blue line) periodicity, highlighting a marginal difference (black
dashed line) between the two configurations. (b) Execution time on GPU with
(red line) and without (blue line) periodicity, where the difference (black
dashed line) is negligible.

To benchmark the performance of the FFT-PIM, we conducted a set of tests on a
workstation equipped with AMD 5950X and NVIDIA 3080 Ti, with the CPU tests
utilizing a single core. The results are shown for a 3D periodicity static
problem without phase shifts. The target object is a cube with dimensions
$D_{x},D_{y},D_{z}=100$ and $L_{x,y,z}=D_{x,y,z}+1$. This problem has
$L_{x,y,z}\approx D_{x,y,z}$ and it represents the maximal computational
complexity as compared to 1D and 2D periodicity cases. The configuration of
the far-zone calculation is set as cubic interpolation with uniform grid size
of $N_{gx},N_{gy},N_{gz}=10$, and $i_{d}=1$. The near-zone is handled using
the procedures from the FastMag micromagnetic simulator [22, 25], which is
modified to allow for an efficient handling of periodicities, as outlined in
Sec. 3.3.

Fig. 7 shows the computational time on CPU for pre-processing needed for
computing the interpolation coefficients and PGF tables. We compare the CPU
time with the case of the same source - observer distributions but without any
periodicity, which demonstrates the overhead due to the periodicity. The
predominant additional time complexity arises from substituting the original
free-space Green’s function, $G_{0}$, with a more complicated $G_{near}$.
Specifically in this 3D case with $i_{d}=1$, $G_{near}$ is $(2\times 1+1)^{3}$
times slower than $G_{0}$. However, the overall discrepancy is marginal,
resulting in only a 15% overhead. This small impact is due to the fact that
tabulating the Green’s function constitutes only a small part of the pre-
processing time, and the time complexities of other segments do not increase
as the periodicity is introduced. We also evaluated 1D and 2D periodic cases,
and obtained similar performance with a reduced overhead.

Fig. 8 shows the execution time for a single evaluation of PSP on CPU and GPU.
When compared to a non-periodic problem, the major execution time increase is
due to the evaluation of the far-zone PSP component, which can be performed
separately from the near-zone evaluation. The presence or absence of
periodicity has minimal impact on the near-zone evaluation since all Green’s
function values are tabulated during the pre-processing step. Consequently, in
the single-core CPU benchmark, only a minor difference is observed between
configurations with and without periodicity. In contrast, in the GPU
benchmark, NVIDIA CUDA multi-stream concurrency enables simultaneous
implementation of both near-zone and far-zone evaluations. This approach
effectively eliminates the overhead induced by the far-zone evaluation, which
is substantially slower, resulting in negligible observable differences
between cases with and without periodicity, even in the 3D periodicity cases.
The 1D and 2D periodicity cases exhibit a similar performance, indicating a
highly efficient implementation for all 1D, 2D, 3D periodicity types.

##  5 Summary and discussion

We introduced an efficient and flexible FFT-PIM for computing PSP. FFT-PIM can
be used for a wide variety of problem types with the same numerical
implementation. It allows computing PSP for a non-uniform source distribution,
works for arbitrary 1D, 2D, and 3D periodicities, can operate with or without
a phase shift between the periodic boundaries, and is applicable to dynamic
and static problems. A requirement is imposed for using FFT-PIM for static
problems with a vanishing phase shift, which is that the source must be
neutral within the unit cell. Such a requirement is natural to many practical
problems, such as problems dealing with electric dipoles, magnetization, and
molecular structures.

FFT-PIM is based on a superposition sum between the source distribution and
PGF. PGF, which originally is defined as an infinite sum that does not have a
rapid convergence, can be represented as an exponentially convergent sum in
terms of spectral series expansions, which are given for the dynamic case and
static case with a periodic phase shift as well as for a static case with no
phase shift. The expansions converge under the assumption of a sufficient
separation between the source and observer points. To lead to a fast
superposition sum computation, PGF and the corresponding PSP are represented
in terms of the near- and far-zone components. The near-zone component
includes a finite, typically small, number of periodic images, whereas the
far-zone component includes all the rest of an infinite number of images. The
far-zone component is recognized as a slowly varying function of the spatial
coordinates, which allows computing the far-zone PSP component by first
computing it on a sparse uniform grid, i.e., at a small number of points, and
then interpolating it to all the required non-uniform observation points by
local interpolation. This process involves defining shifted source and
observer grids to allow for a rapidly convergence of the spectral sums for
PGF. This process requires a small number of operations and low computational
cost, while allowing for a rapid converge. The near-zone PSP component can be
evaluated by any fast method, but it required accounting for the additional
near-zone images, which may substantially increase the computational cost and
memory consumption. An approach based on the box-adaptive adaptive integral
method is presented that allows evaluating the superposition sum for the far-
zone component rapidly based on including essential images only at the
required locations near the periodic boundaries. The result is an approach
that has an overall computational cost of $O(N\log N)$ and memory consumption
of $O(N)$. The presented results demonstrate the high convergence, accuracy,
and computational performance of the FFT-PIM.

##  6 Code availability

We open-sourced our far-zone component code package Periodic Unit Far Field
(PUFF) under Apache-2.0 license on GitHub (https://github.com/UCSD-
CEM/PeriodicUnitFarField).

##  7 Acknowledgments

This work was supported in part by the Quantum Materials for Energy Efficient
Neuromorphic-Computing (Q-MEEN-C), an Energy Frontier Research Center funded
by the U.S. Department of Energy, Office of Science, Basic Energy Sciences
under Award No. DESC0019273. The work was also supported in part by Binational
Science Foundation, grant #2022346. The work used Purdue Anvil cluster at
Rosen Center for Advanced Computing (RCAC) in Purdue University through
allocation ASC200042 from the Advanced Cyberinfrastructure Coordination
Ecosystem: Services & Support (ACCESS) program [26], which is supported by
National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and
#2138296.

## References

  * [1] R. J. Mailloux, Phased array antenna handbook.  Artech house, 2017.
  * [2] D. Kalkstein and P. Soven, “A green’s function theory of surface states,” Surface Science, vol. 26, no. 1, pp. 85–99, 1971.
  * [3] Y. K. Sirenko and S. Strom, “Modern theory of gratings,” Resonant Scattering: Analysis Techniques and Phenomena, 2010.
  * [4] S. Peng, T. Tamir, and H. L. Bertoni, “Theory of periodic dielect waveguides,” IEEE transactions on microwave theory and techniques, vol. 23, no. 1, pp. 123–133, 1975.
  * [5] L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” Journal of Computational Physics, vol. 73, no. 2, pp. 325–348, 1987.
  * [6] A. D. Baczewski, D. L. Dault, and B. Shanker, “Accelerated cartesian expansions for the rapid solution of periodic multiscale problems,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 9, pp. 4281–4290, 2012.
  * [7] W. Hackbusch and B. N. Khoromskij, “A sparse h-matrix arithmetic. part ii: application to multi-dimensional problems,” Computing, vol. 64, pp. 21–47, 2000.
  * [8] A. Boag, E. Michielssen, and A. Brandt, “Nonuniform polar grid algorithm for fast field evaluation,” IEEE Antennas and Wireless Propagation Letters, vol. 1, pp. 142–145, 2002.
  * [9] S. Li, B. Livshitz, and V. Lomakin, “Fast evaluation of helmholtz potential on graphics processing units (gpus),” Journal of Computational Physics, vol. 229, no. 22, pp. 8463–8483, 2010.
  * [10] J. Meng, A. Boag, V. Lomakin, and E. Michielssen, “A multilevel cartesian non-uniform grid time domain algorithm,” Journal of Computational Physics, vol. 229, no. 22, pp. 8430–8444, 2010.
  * [11] B. Livshitz, A. Boag, H. N. Bertram, and V. Lomakin, “Nonuniform grid algorithm for fast calculation of magnetostatic interactions in micromagnetics,” Journal of Applied Physics, vol. 105, p. 07D541, 03 2009\.
  * [12] C. Liu, K. Aygün, and A. E. Yılmaz, “A parallel fft-accelerated layered-medium integral-equation solver for electronic packages,” International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, vol. 33, no. 2, p. e2684, 2020.  e2684 JNM-19-0067.R1.
  * [13] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, “Aim: Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems,” Radio Science, vol. 31, no. 5, pp. 1225–1251, 1996.
  * [14] J. Phillips and J. White, “A precorrected-fft method for electrostatic analysis of complicated 3-d structures,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 16, no. 10, pp. 1059–1072, 1997.
  * [15] F. Capolino, D. Wilton, and W. Johnson, “Efficient computation of the 3d green’s function for the helmholtz operator for a linear array of point sources using the ewald method,” Journal of Computational Physics, vol. 223, no. 1, pp. 250–261, 2007.
  * [16] S. Li, D. A. Van Orden, and V. Lomakin, “Fast periodic interpolation method for periodic unit cell problems,” IEEE Transactions on Antennas and Propagation, vol. 58, no. 12, pp. 4005–4014, 2010.
  * [17] D. Van Orden and V. Lomakin, “Rapidly convergent representations for 2d and 3d green’s functions for a linear periodic array of dipole sources,” IEEE Transactions on Antennas and Propagation, vol. 57, no. 7, pp. 1973–1984, 2009\.
  * [18] Z. Hu, “Infinite boundary terms of ewald sums and pairwise interactions for electrostatics in bulk and at interfaces,” Journal of Chemical Theory and Computation, vol. 10, no. 12, pp. 5254–5264, 2014.  PMID: 26583209.
  * [19] S. L. Marshall, “A periodic green function for calculation of coloumbic lattice potentials,” Journal of Physics: Condensed Matter, vol. 12, p. 4575, may 2000.
  * [20] K. M. Lebecki, M. J. Donahue, and M. W. Gutowski, “Periodic boundary conditions for demagnetization interactions in micromagnetic simulations,” Journal of Physics D: Applied Physics, vol. 41, p. 175005, aug 2008.
  * [21] F. Bruckner, A. Ducevic, P. Heistracher, C. Abert, and D. Suess, “Strayfield calculation for micromagnetic simulations using true periodic boundary conditions,” Scientific Reports, vol. 11, p. 9202, Apr 2021.
  * [22] R. Chang, S. Li, M. V. Lubarda, B. Livshitz, and V. Lomakin, “FastMag: Fast micromagnetic simulator for complex magnetic structures (invited),” Journal of Applied Physics, vol. 109, p. 07D358, 04 2011.
  * [23] A. Bagnérés‐Viallix and P. Baras, “On a method to calculate the demagnetizing field in a micromagnetic structure,” Journal of Applied Physics, vol. 69, pp. 4599–4601, 04 1991.
  * [24] J. Stöhr and H. C. Siegmann, “Magnetism,” Solid-State Sciences. Springer, Berlin, Heidelberg, vol. 5, p. 236, 2006.
  * [25] S. Li, R. Chang, A. Boag, and V. Lomakin, “Fast electromagnetic integral-equation solvers on graphics processing units,” IEEE Antennas and Propagation Magazine, vol. 54, no. 5, pp. 71–87, 2012.
  * [26] T. J. Boerner, S. Deems, T. R. Furlani, S. L. Knuth, and J. Towns, “Access: Advancing innovation: Nsf’s advanced cyberinfrastructure coordination ecosystem: Services & support,” in Practice and Experience in Advanced Research Computing, PEARC ’23, (New York, NY, USA), p. 173–176, Association for Computing Machinery, 2023.