qpms/notes/ewald.lyx

1418 lines
33 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#LyX 2.1 created this file. For more info see http://www.lyx.org/
\lyxformat 474
\begin_document
\begin_header
\textclass article
\begin_preamble
\usepackage{unicode-math}
% Toto je trik, jimž se z fontspec získá familyname pro následující
\ExplSyntaxOn
\DeclareExpandableDocumentCommand{\getfamilyname}{m}
{
\use:c { g__fontspec_ \cs_to_str:N #1 _family }
}
\ExplSyntaxOff
% definujeme novou rodinu, jež se volá pomocí \MyCyr pro běžné použití, avšak pro účely \DeclareSymbolFont je nutno získat název pomocí getfamilyname definovaného výše
\newfontfamily\MyCyr{CMU Serif}
\DeclareSymbolFont{cyritletters}{EU1}{\getfamilyname\MyCyr}{m}{it}
\newcommand{\makecyrmathletter}[1]{%
\begingroup\lccode`a=#1\lowercase{\endgroup
\Umathcode`a}="0 \csname symcyritletters\endcsname\space #1
}
\count255="409
\loop\ifnum\count255<"44F
\advance\count255 by 1
\makecyrmathletter{\count255}
\repeat
\renewcommand{\lyxmathsym}[1]{#1}
\end_preamble
\use_default_options true
\maintain_unincluded_children false
\language english
\language_package default
\inputencoding auto
\fontencoding global
\font_roman TeX Gyre Pagella
\font_sans default
\font_typewriter default
\font_math default
\font_default_family default
\use_non_tex_fonts true
\font_sc false
\font_osf true
\font_sf_scale 100
\font_tt_scale 100
\graphics default
\default_output_format pdf4
\output_sync 0
\bibtex_command default
\index_command default
\paperfontsize 10
\spacing single
\use_hyperref true
\pdf_title "Accelerating lattice mode calculations with T-matrix method"
\pdf_author "Marek Nečada"
\pdf_bookmarks true
\pdf_bookmarksnumbered false
\pdf_bookmarksopen false
\pdf_bookmarksopenlevel 1
\pdf_breaklinks false
\pdf_pdfborder false
\pdf_colorlinks false
\pdf_backref false
\pdf_pdfusetitle true
\papersize a5paper
\use_geometry true
\use_package amsmath 1
\use_package amssymb 1
\use_package cancel 1
\use_package esint 1
\use_package mathdots 1
\use_package mathtools 1
\use_package mhchem 1
\use_package stackrel 1
\use_package stmaryrd 1
\use_package undertilde 1
\cite_engine basic
\cite_engine_type default
\biblio_style plain
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
\justification true
\use_refstyle 1
\index Index
\shortcut idx
\color #008000
\end_index
\leftmargin 2cm
\topmargin 2cm
\rightmargin 2cm
\bottommargin 2cm
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\quotes_language english
\papercolumns 1
\papersides 1
\paperpagestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header
\begin_body
\begin_layout Standard
\begin_inset FormulaMacro
\newcommand{\uoft}[1]{\mathfrak{F}#1}
\end_inset
\begin_inset FormulaMacro
\newcommand{\uaft}[1]{\mathfrak{\mathbb{F}}#1}
\end_inset
\begin_inset FormulaMacro
\newcommand{\usht}[2]{\mathbb{S}_{#1}#2}
\end_inset
\begin_inset FormulaMacro
\newcommand{\bsht}[2]{\mathrm{S}_{#1}#2}
\end_inset
\begin_inset FormulaMacro
\newcommand{\pht}[2]{\mathfrak{\mathbb{H}}_{#1}#2}
\end_inset
\begin_inset FormulaMacro
\newcommand{\vect}[1]{\mathbf{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ud}{\mathrm{d}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\basis}[1]{\mathfrak{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\dc}[1]{Ш_{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\rec}[1]{#1^{-1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\recb}[1]{#1^{\widehat{-1}}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ints}{\mathbb{Z}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\nats}{\mathbb{N}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\reals}{\mathbb{R}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ush}[2]{Y_{#1,#2}}
\end_inset
\end_layout
\begin_layout Title
Accelerating lattice mode calculations with
\begin_inset Formula $T$
\end_inset
-matrix method
\end_layout
\begin_layout Author
Marek Nečada
\end_layout
\begin_layout Abstract
The
\begin_inset Formula $T$
\end_inset
-matrix approach is the method of choice for simulating optical response
of a reasonably small system of compact linear scatterers on isotropic
background.
However, its direct utilisation for problems with infinite lattices is
problematic due to slowly converging sums over the lattice.
Here I develop a way to compute the problematic sums in the reciprocal
space, making the
\begin_inset Formula $T$
\end_inset
-matrix method very suitable for infinite periodic systems as well.
\end_layout
\begin_layout Section
Formulation of the problem
\end_layout
\begin_layout Standard
Assume a system of compact EM scatterers in otherwise homogeneous and isotropic
medium, and assume that the system, i.e.
both the medium and the scatterers, have linear response.
A scattering problem in such system can be written as
\begin_inset Formula
\[
A_{α}=T_{α}P_{α}=T_{α}(\sum_{β}S_{α\leftarrowβ}A_{β}+P_{0α})
\]
\end_inset
where
\begin_inset Formula $T_{α}$
\end_inset
is the
\begin_inset Formula $T$
\end_inset
-matrix for scatterer α,
\begin_inset Formula $A_{α}$
\end_inset
is its vector of the scattered wave expansion coefficient (the multipole
indices are not explicitely indicated here) and
\begin_inset Formula $P_{α}$
\end_inset
is the local expansion of the incoming sources.
\begin_inset Formula $S_{α\leftarrowβ}$
\end_inset
is ...
and ...
is ...
\end_layout
\begin_layout Standard
...
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\sum_{β}(\delta_{αβ}-T_{α}S_{α\leftarrowβ})A_{β}=T_{α}P_{0α}.
\]
\end_inset
\end_layout
\begin_layout Standard
Now suppose that the scatterers constitute an infinite lattice
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\sum_{\vect bβ}(\delta_{\vect{ab}}\delta_{αβ}-T_{\vect aα}S_{\vect aα\leftarrow\vect bβ})A_{\vect bβ}=T_{\vect aα}P_{0\vect aα}.
\]
\end_inset
Due to the periodicity, we can write
\begin_inset Formula $S_{\vect aα\leftarrow\vect bβ}=S_{α\leftarrowβ}(\vect b-\vect a)$
\end_inset
and
\begin_inset Formula $T_{\vect aα}=T_{\alpha}$
\end_inset
.
In order to find lattice modes, we search for solutions with zero RHS
\begin_inset Formula
\[
\sum_{\vect bβ}(\delta_{\vect{ab}}\delta_{αβ}-T_{α}S_{\vect aα\leftarrow\vect bβ})A_{\vect bβ}=0
\]
\end_inset
and we assume periodic solution
\begin_inset Formula $A_{\vect b\beta}(\vect k)=A_{\vect a\beta}e^{i\vect k\cdot\vect r_{\vect b-\vect a}}$
\end_inset
, yielding
\begin_inset Formula
\begin{eqnarray*}
\sum_{\vect bβ}(\delta_{\vect{ab}}\delta_{αβ}-T_{α}S_{\vect aα\leftarrow\vect bβ})A_{\vect a\beta}\left(\vect k\right)e^{i\vect k\cdot\vect r_{\vect b-\vect a}} & = & 0,\\
\sum_{\vect bβ}(\delta_{\vect{0b}}\delta_{αβ}-T_{α}S_{\vect 0α\leftarrow\vect bβ})A_{\vect 0\beta}\left(\vect k\right)e^{i\vect k\cdot\vect r_{\vect b}} & = & 0,\\
\sum_{β}(\delta_{αβ}-T_{α}\underbrace{\sum_{\vect b}S_{\vect 0α\leftarrow\vect bβ}e^{i\vect k\cdot\vect r_{\vect b}}}_{W_{\alpha\beta}(\vect k)})A_{\vect 0\beta}\left(\vect k\right) & = & 0,\\
A_{\vect 0\alpha}\left(\vect k\right)-T_{α}\sum_{\beta}W_{\alpha\beta}\left(\vect k\right)A_{\vect 0\beta}\left(\vect k\right) & = & 0.
\end{eqnarray*}
\end_inset
Therefore, in order to solve the modes, we need to compute the
\begin_inset Quotes eld
\end_inset
lattice Fourier transform
\begin_inset Quotes erd
\end_inset
of the translation operator,
\begin_inset Formula
\begin{equation}
W_{\alpha\beta}(\vect k)\equiv\sum_{\vect b}S_{\vect 0α\leftarrow\vect bβ}e^{i\vect k\cdot\vect r_{\vect b}}.\label{eq:W definition}
\end{equation}
\end_inset
\end_layout
\begin_layout Section
Computing the Fourier sum of the translation operator
\end_layout
\begin_layout Standard
The problem evaluating
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:W definition"
\end_inset
is the asymptotic behaviour of the translation operator,
\begin_inset Formula $S_{\vect 0α\leftarrow\vect bβ}\sim\left|\vect r_{\vect b}\right|^{-1}e^{ik_{0}\left|\vect r_{\vect b}\right|}$
\end_inset
that makes the convergence of the sum quite problematic for any
\begin_inset Formula $d>1$
\end_inset
-dimensional lattice.
\begin_inset Foot
status open
\begin_layout Plain Layout
Note that
\begin_inset Formula $d$
\end_inset
here is dimensionality of the lattice, not the space it lies in, which
I for certain reasons assume to be three.
(TODO few notes on integration and reciprocal lattices in some appendix)
\end_layout
\end_inset
In electrostatics, one can solve this problem with Ewald summation.
Its basic idea is that if what asymptoticaly decays poorly in the direct
space, will perhaps decay fast in the Fourier space.
I use the same idea here, but everything will be somehow harder than in
electrostatics.
\end_layout
\begin_layout Standard
Let us re-express the sum in
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:W definition"
\end_inset
in terms of integral with a delta comb
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
W_{\alpha\beta}(\vect k)=\int\ud^{d}\vect r\dc{\basis u}(\vect r)S(\vect r_{\alpha}\leftarrow\vect r+\vect r_{\beta})e^{i\vect k\cdot\vect r}.\label{eq:W integral}
\end{equation}
\end_inset
The translation operator
\begin_inset Formula $S$
\end_inset
is now a function defined in the whole 3d space;
\begin_inset Formula $\vect r_{\alpha},\vect r_{\beta}$
\end_inset
are the displacements of scatterers
\begin_inset Formula $\alpha$
\end_inset
and
\begin_inset Formula $\beta$
\end_inset
in a unit cell.
The arrow notation
\begin_inset Formula $S(\vect r_{\alpha}\leftarrow\vect r+\vect r_{\beta})$
\end_inset
means
\begin_inset Quotes eld
\end_inset
translation operator for spherical waves originating in
\begin_inset Formula $\vect r+\vect r_{\beta}$
\end_inset
evaluated in
\begin_inset Formula $\vect r_{\alpha}$
\end_inset
\begin_inset Quotes erd
\end_inset
and obviously
\begin_inset Formula $S$
\end_inset
is in fact a function of a single 3d argument,
\begin_inset Formula $S(\vect r_{\alpha}\leftarrow\vect r+\vect r_{\beta})=S(\vect 0\leftarrow\vect r+\vect r_{\beta}-\vect r_{\alpha})=S(-\vect r-\vect r_{\beta}+\vect r_{\alpha}\leftarrow\vect 0)=S(-\vect r-\vect r_{\beta}+\vect r_{\alpha})$
\end_inset
.
Expression
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:W integral"
\end_inset
can be rewritten as
\begin_inset Formula
\[
W_{\alpha\beta}(\vect k)=\left(2\pi\right)^{\frac{d}{2}}\uaft{(\dc{\basis u}S(\vect{\bullet}-\vect r_{\beta}+\vect r_{\alpha}\leftarrow\vect 0))\left(\vect k\right)}
\]
\end_inset
where changed the sign of
\begin_inset Formula $\vect r/\vect{\bullet}$
\end_inset
has been swapped under integration, utilising evenness of
\begin_inset Formula $\dc{\basis u}$
\end_inset
.
Fourier transform of product is convolution of Fourier transforms, so (using
formula
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:Dirac comb uaFt"
\end_inset
for the Fourier transform of Dirac comb)
\begin_inset Formula
\begin{eqnarray}
W_{\alpha\beta}(\vect k) & = & \left(\left(\uaft{\dc{\basis u}}\right)\ast\left(\uaft{S(\vect{\bullet}-\vect r_{\beta}+\vect r_{\alpha}\leftarrow\vect 0)}\right)\right)(\vect k)\nonumber \\
& = & \frac{\left|\det\recb{\basis u}\right|}{\left(2\pi\right)^{\frac{d}{2}}}\left(\dc{\recb{\basis u}}^{(d)}\ast\left(\uaft{S(\vect{\bullet}-\vect r_{\beta}+\vect r_{\alpha}\leftarrow\vect 0)}\right)\right)\left(\vect k\right)\nonumber \\
& = & \frac{\left|\det\rec{\basis u}\right|}{\left(2\pi\right)^{\frac{d}{2}}}\sum_{\vect K\in\recb{\basis u}\ints^{d}}\left(\uaft{S(\vect{\bullet}-\vect r_{\beta}+\vect r_{\alpha}\leftarrow\vect 0)}\right)\left(\vect k-\vect K\right).\label{eq:W sum in reciprocal space}
\end{eqnarray}
\end_inset
\begin_inset Note Note
status open
\begin_layout Plain Layout
Factor
\begin_inset Formula $\left(2\pi\right)^{\frac{d}{2}}$
\end_inset
cancels out with the
\begin_inset Formula $\left(2\pi\right)^{-\frac{d}{2}}$
\end_inset
factor appearing in the convolution/product formula in the unitary angular
momentum convention.
\end_layout
\end_inset
As such, this is not extremely helpful because the the
\emph on
whole
\emph default
translation operator
\begin_inset Formula $S$
\end_inset
has singularities in origin, hence its Fourier transform
\begin_inset Formula $\uaft S$
\end_inset
will decay poorly.
\end_layout
\begin_layout Standard
However, Fourier transform is linear, so we can in principle separate
\begin_inset Formula $S$
\end_inset
in two parts,
\begin_inset Formula $S=S^{\textup{L}}+S^{\textup{S}}$
\end_inset
.
\begin_inset Formula $S^{\textup{S}}$
\end_inset
is a short-range part that decays sufficiently fast with distance so that
its direct-space lattice sum converges well;
\begin_inset Formula $S^{\textup{S}}$
\end_inset
must as well contain all the singularities of
\begin_inset Formula $S$
\end_inset
in the origin.
The other part,
\begin_inset Formula $S^{\textup{L}}$
\end_inset
, will retain all the slowly decaying terms of
\begin_inset Formula $S$
\end_inset
but it also has to be smooth enough in the origin, so that its Fourier
transform
\begin_inset Formula $\uaft{S^{\textup{L}}}$
\end_inset
decays fast enough.
(The same idea lies behind the Ewald summation in electrostatics.) Using
the linearity of Fourier transform and formulae
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:W definition"
\end_inset
and
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:W sum in reciprocal space"
\end_inset
, the operator
\begin_inset Formula $W_{\alpha\beta}$
\end_inset
can then be re-expressed as
\begin_inset Formula
\begin{eqnarray}
W_{\alpha\beta}\left(\vect k\right) & = & W_{\alpha\beta}^{\textup{S}}\left(\vect k\right)+W_{\alpha\beta}^{\textup{L}}\left(\vect k\right)\nonumber \\
W_{\alpha\beta}^{\textup{S}}\left(\vect k\right) & = & \sum_{\vect R\in\basis u\ints^{d}}S^{\textup{S}}(\vect 0\leftarrow\vect R+\vect r_{\beta}-\vect r_{\alpha})e^{i\vect k\cdot\vect R}\label{eq:W Short definition}\\
W_{\alpha\beta}^{\textup{L}}\left(\vect k\right) & = & \frac{\left|\det\rec{\basis u}\right|}{\left(2\pi\right)^{\frac{d}{2}}}\sum_{\vect K\in\recb{\basis u}\ints^{d}}\left(\uaft{S^{\textup{L}}(\vect{\bullet}-\vect r_{\beta}+\vect r_{\alpha}\leftarrow\vect 0)}\right)\left(\vect k-\vect K\right)\label{eq:W Long definition}
\end{eqnarray}
\end_inset
where both sums should converge nicely.
\end_layout
\begin_layout Section
Finding a good decomposition
\end_layout
\begin_layout Standard
The remaining challenge is therefore finding a suitable decomposition
\begin_inset Formula $S^{\textup{L}}+S^{\textup{S}}$
\end_inset
such that both
\begin_inset Formula $S^{\textup{S}}$
\end_inset
and
\begin_inset Formula $\uaft{S^{\textup{L}}}$
\end_inset
decay fast enough with distance and are expressable analytically.
With these requirements, I do not expect to find gaussian asymptotics as
in the electrostatic Ewald formula—having
\begin_inset Formula $\sim x^{-t}$
\end_inset
,
\begin_inset Formula $t>d$
\end_inset
asymptotics would be nice, making the sums in
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:W Short definition"
\end_inset
,
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:W Long definition"
\end_inset
absolutely convergent.
\end_layout
\begin_layout Standard
The translation operator
\begin_inset Formula $S$
\end_inset
for compact scatterers in 3d can be expressed as
\begin_inset Formula
\[
S_{l',m',t'\leftarrow l,m,t}\left(\vect r\leftarrow\vect 0\right)=\sum_{p}c_{p}^{l',m',t'\leftarrow l,m,t}\ush p{m'-m}\left(\theta_{\vect r},\phi_{\vect r}\right)z_{p}^{(J)}\left(k_{0}\left|\vect r\right|\right)
\]
\end_inset
where
\begin_inset Formula $Y_{l,m}\left(\theta,\phi\right)$
\end_inset
are the spherical harmonics,
\begin_inset Formula $z_{p}^{(J)}\left(r\right)$
\end_inset
some of the Bessel or Hankel functions (probably
\begin_inset Formula $h_{p}^{(1)}$
\end_inset
in the meaningful cases; TODO) and
\begin_inset Formula $c_{p}^{l,m,t\leftarrow l',m',t'}$
\end_inset
are some ugly but known coefficients (REF Xu 1996, eqs.
76,77).
\end_layout
\begin_layout Standard
The spherical Hankel functions can be expressed analytically as (REF DLMF
10.49.6, 10.49.1)
\begin_inset Formula
\begin{equation}
h_{n}^{(1)}(r)=e^{ir}\sum_{k=0}^{n}\frac{i^{k-n-1}}{r^{k+1}}\frac{\left(n+k\right)!}{2^{k}k!\left(n-k\right)!},\label{eq:spherical Hankel function series}
\end{equation}
\end_inset
so if we find a way to deal with the radial functions
\begin_inset Formula $s_{q}(r)=e^{ik_{0}r}\left(k_{0}r\right)^{-q}$
\end_inset
,
\begin_inset Formula $q=1,2$
\end_inset
in 2d case or
\begin_inset Formula $q=1,2,3$
\end_inset
in 3d case, we get absolutely convergent summations in the direct space.
\end_layout
\begin_layout Subsection
2d
\end_layout
\begin_layout Standard
Assume that all scatterers are placed in the plane
\begin_inset Formula $\vect z=0$
\end_inset
, so that the 2d Fourier transform of the long-range part of the translation
operator in terms of Hankel transforms, according to
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:Fourier v. Hankel tf 2d"
\end_inset
, reads
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{multline*}
\uaft{S_{l',m',t'\leftarrow l,m,t}^{\textup{L}}\left(\vect{\bullet}\leftarrow\vect 0\right)}(\vect k)=\\
\sum_{p}c_{p}^{l',m',t'\leftarrow l,m,t}\ush p{m'-m}\left(\frac{\pi}{2},0\right)e^{i(m'-m)\phi}i^{m'-m}\pht{m'-m}{h_{p}^{(1)\textup{L}}\left(k_{0}\vect{\bullet}\right)}\left(\left|\vect k\right|\right)
\end{multline*}
\end_inset
Here
\begin_inset Formula $h_{p}^{(1)\textup{L}}=h_{p}^{(1)}-h_{p}^{(1)\textup{S}}$
\end_inset
is a long range part of a given spherical Hankel function which has to
be found and which contains all the terms with far-field (
\begin_inset Formula $r\to\infty$
\end_inset
) asymptotics proportional to
\begin_inset Formula $\sim e^{ik_{0}r}\left(k_{0}r\right)^{-q}$
\end_inset
,
\begin_inset Formula $q\le Q$
\end_inset
where
\begin_inset Formula $Q$
\end_inset
is at least two in order to achieve absolute convergence of the direct-space
sum, but might be higher in order to speed the convergence up.
\end_layout
\begin_layout Standard
Obviously, all the terms
\begin_inset Formula $\propto s_{q}(r)=e^{ik_{0}r}\left(k_{0}r\right)^{-q}$
\end_inset
,
\begin_inset Formula $q>Q$
\end_inset
of the spherical Hankel function
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:spherical Hankel function series"
\end_inset
can be kept untouched as part of
\begin_inset Formula $h_{p}^{(1)\textup{S}}$
\end_inset
, as they decay fast enough.
\end_layout
\begin_layout Standard
The remaining task is therefore to find a suitable decomposition of
\begin_inset Formula $s_{q}(r)=e^{ik_{0}r}\left(k_{0}r\right)^{-q}$
\end_inset
,
\begin_inset Formula $q\le Q$
\end_inset
into short-range and long-range parts,
\begin_inset Formula $s_{q}(r)=s_{q}^{\textup{S}}(r)+s_{q}^{\textup{L}}(r)$
\end_inset
, such that
\begin_inset Formula $s_{q}^{\textup{L}}(r)$
\end_inset
contains all the slowly decaying asymptotics and its Hankel transforms
decay desirably fast as well,
\begin_inset Formula $\pht n{s_{q}^{\textup{L}}}\left(k\right)=o(z^{-Q})$
\end_inset
,
\begin_inset Formula $z\to\infty$
\end_inset
.
The latter requirement calls for suitable regularisation functions—
\begin_inset Formula $s_{q}^{\textup{L}}$
\end_inset
must be sufficiently smooth in the origin, so that
\begin_inset Formula
\begin{equation}
\pht n{s_{q}^{\textup{L}}}\left(k\right)=\int_{0}^{\infty}s_{q}^{\textup{L}}\left(r\right)rJ_{n}\left(kr\right)\ud r=\int_{0}^{\infty}s_{q}\left(r\right)\rho\left(r\right)rJ_{n}\left(kr\right)\ud r\label{eq:2d long range regularisation problem statement}
\end{equation}
\end_inset
exists and decays fast enough.
\begin_inset Formula $J_{\nu}(r)\sim\left(r/2\right)^{\nu}/\Gamma\left(\nu+1\right)$
\end_inset
(REF DLMF 10.7.3) near the origin, so the regularisation function should
be
\begin_inset Formula $\rho(r)=o(r^{q-n-1})$
\end_inset
only to make
\begin_inset Formula $\pht n{s_{q}^{\textup{L}}}$
\end_inset
converge.
The additional decay speed requirement calls for at least
\begin_inset Formula $\rho(r)=o(r^{q-n+Q-1})$
\end_inset
, I guess.
At the same time,
\begin_inset Formula $\rho(r)$
\end_inset
must converge fast enough to one for
\begin_inset Formula $r\to\infty$
\end_inset
.
\end_layout
\begin_layout Standard
The electrostatic Ewald summation uses regularisation with
\begin_inset Formula $1-e^{-cr^{2}}$
\end_inset
.
However, such choice does not seem to lead to an analytical solution for
the current problem
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:2d long range regularisation problem statement"
\end_inset
.
But it turns out that the family of functions
\begin_inset Formula
\[
\rho_{\kappa,c}(r)\equiv\left(1-e^{-cr}\right)^{\text{\kappa}},\quad c>0,\kappa\in\nats
\]
\end_inset
leads to satisfactory results, as will be shown below.
\end_layout
\begin_layout Subsubsection
Hankel transforms of the long-range parts
\end_layout
\begin_layout Subsection
3d
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{multline*}
\uaft{S_{l',m',t'\leftarrow l,m,t}\left(\vect{\bullet}\leftarrow\vect 0\right)}(\vect k)=\\
\sum_{p}c_{p}^{l',m',t'\leftarrow l,m,t}\ush p{m'-m}\left(\theta_{\vect k},\phi_{\vect k}\right)\left(-i\right)^{p}\usht p{z_{p}^{(J)}}\left(\left|\vect k\right|\right)
\end{multline*}
\end_inset
\end_layout
\begin_layout Section
(Appendix) Fourier vs.
Hankel transform
\end_layout
\begin_layout Subsection
Three dimensions
\end_layout
\begin_layout Standard
Given a nice enough function
\begin_inset Formula $f$
\end_inset
of a real 3d variable, assume its factorisation into radial and angular
parts
\begin_inset Formula
\[
f(\vect r)=\sum_{l,m}f_{l,m}(\left|\vect r\right|)\ush lm\left(\theta_{\vect r},\phi_{\vect r}\right).
\]
\end_inset
Acording to (REF Baddour 2010, eqs.
13, 16), its Fourier transform can then be expressed in terms of Hankel
transforms (CHECK normalisation of
\begin_inset Formula $j_{n}$
\end_inset
, REF Baddour (1))
\begin_inset Formula
\[
\uaft f(\vect k)=\frac{4\pi}{\left(2\pi\right)^{\frac{3}{2}}}\sum_{l,m}\left(-i\right)^{l}\left(\bsht{f_{l,m}}{}\right)\left(\left|\vect k\right|\right)\ush lm\left(\theta_{\vect k},\phi_{\vect k}\right)
\]
\end_inset
where the spherical Hankel transform
\begin_inset Formula $\bsht l{}$
\end_inset
of degree
\begin_inset Formula $l$
\end_inset
is defined as (REF Baddour eq.
2)
\begin_inset Formula
\[
\bsht lg(k)\equiv\int_{0}^{\infty}\ud r\,r^{2}g(r)j_{l}\left(kr\right).
\]
\end_inset
Using this convention, the inverse spherical Hankel transform is given by
(REF Baddour eq.
3)
\begin_inset Formula
\[
g(r)=\frac{2}{\pi}\int_{0}^{\infty}\ud k\,k^{2}\bsht lg(k)j_{l}(k),
\]
\end_inset
so it is not unitary.
\end_layout
\begin_layout Standard
An unitary convention would look like this:
\begin_inset Formula
\begin{equation}
\usht lg(k)\equiv\sqrt{\frac{2}{\pi}}\int_{0}^{\infty}\ud r\,r^{2}g(r)j_{l}\left(kr\right).\label{eq:unitary 3d Hankel tf definition}
\end{equation}
\end_inset
Then
\begin_inset Formula $\usht l{}^{-1}=\usht l{}$
\end_inset
and the unitary, angular-momentum Fourier transform reads
\begin_inset Formula
\begin{eqnarray}
\uaft f(\vect k) & = & \frac{4\pi}{\left(2\pi\right)^{\frac{3}{2}}}\sqrt{\frac{\pi}{2}}\sum_{l,m}\left(-i\right)^{l}\left(\usht l{f_{l,m}}\right)\left(\left|\vect k\right|\right)\ush lm\left(\theta_{\vect k},\phi_{\vect k}\right)\nonumber \\
& = & \sum_{l,m}\left(-i\right)^{l}\left(\usht l{f_{l,m}}\right)\left(\left|\vect k\right|\right)\ush lm\left(\theta_{\vect k},\phi_{\vect k}\right).\label{eq:Fourier v. Hankel tf 3d}
\end{eqnarray}
\end_inset
Cool.
\end_layout
\begin_layout Subsection
Two dimensions
\end_layout
\begin_layout Standard
Similarly in 2d, let the expansion of
\begin_inset Formula $f$
\end_inset
be
\begin_inset Formula
\[
f\left(\vect r\right)=\sum_{m}f_{m}\left(\left|\vect r\right|\right)e^{im\phi_{\vect r}},
\]
\end_inset
its Fourier transform is then (CHECK this, it is taken from the Wikipedia
article on Hankel transform)
\begin_inset Formula
\begin{equation}
\uaft f\left(\vect k\right)=\sum_{m}i^{m}e^{im\phi_{\vect k}}\pht mf_{m}\left(\left|\vect k\right|\right)\label{eq:Fourier v. Hankel tf 2d}
\end{equation}
\end_inset
where the Hankel transform of order
\begin_inset Formula $m$
\end_inset
is defined as
\begin_inset Formula
\begin{equation}
\pht mg\left(k\right)=\int_{0}^{\infty}\ud r\,g(r)J_{m}(kr)r\label{eq:unitary 2d Hankel tf definition}
\end{equation}
\end_inset
which is already self-inverse,
\begin_inset Formula $\pht m{}^{-1}=\pht m{}$
\end_inset
(hence also unitary).
\end_layout
\begin_layout Section
(Appendix) Multidimensional Dirac comb
\end_layout
\begin_layout Subsection
1D
\end_layout
\begin_layout Standard
This is all from Wikipedia
\end_layout
\begin_layout Subsubsection
Definitions
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{eqnarray*}
Ш(t) & \equiv & \sum_{k=-\infty}^{\infty}\delta(t-k)\\
Ш_{T}(t) & \equiv & \sum_{k=-\infty}^{\infty}\delta(t-kT)=\frac{1}{T}Ш\left(\frac{t}{T}\right)
\end{eqnarray*}
\end_inset
\end_layout
\begin_layout Subsubsection
Fourier series representation
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
Ш_{T}(t)=\sum_{n=-\infty}^{\infty}e^{2\pi int/T}\label{eq:1D Dirac comb Fourier series}
\end{equation}
\end_inset
\end_layout
\begin_layout Subsubsection
Fourier transform
\end_layout
\begin_layout Standard
With unitary ordinary frequency Ft., i.e.
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\uoft f(\vect{\xi})\equiv\int_{\mathbb{R}^{n}}f(\vect x)e^{-2\pi i\vect x\cdot\vect{\xi}}\ud^{n}\vect x
\]
\end_inset
we have
\begin_inset Formula
\begin{equation}
\uoft{Ш_{T}}(f)=\frac{1}{T}Ш_{\frac{1}{T}}(f)=\sum_{n=-\infty}^{\infty}e^{-i2\pi fnT}\label{eq:1D Dirac comb Ft ordinary freq}
\end{equation}
\end_inset
and with unitary angular frequency Ft., i.e.
\begin_inset Formula
\[
\uaft f(\vect k)\equiv\frac{1}{\left(2\pi\right)^{n/2}}\int_{\mathbb{R}^{n}}f(\vect x)e^{-i\vect x\cdot\vect k}\ud^{n}\vect x
\]
\end_inset
we have (CHECK)
\begin_inset Formula
\[
\uaft{Ш_{T}}(\omega)=\frac{\sqrt{2\pi}}{T}Ш_{\frac{2\pi}{T}}(\omega)=\frac{1}{\sqrt{2\pi}}\sum_{n=-\infty}^{\infty}e^{-i\omega nT}
\]
\end_inset
\end_layout
\begin_layout Subsection
Dirac comb for multidimensional lattices
\end_layout
\begin_layout Subsubsection
Definitions
\end_layout
\begin_layout Standard
Let
\begin_inset Formula $d$
\end_inset
be the dimensionality of the real vector space in question, and let
\begin_inset Formula $\basis u\equiv\left\{ \vect u_{i}\right\} _{i=1}^{d}$
\end_inset
denote a basis for some lattice in that space.
Let the corresponding lattice delta comb be
\begin_inset Formula
\[
\dc{\basis u}\left(\vect x\right)\equiv\sum_{n_{1}=-\infty}^{\infty}\ldots\sum_{n_{d}=-\infty}^{\infty}\delta\left(\vect x-\sum_{i=1}^{d}n_{i}\vect u_{i}\right).
\]
\end_inset
\end_layout
\begin_layout Standard
Furthemore, let
\begin_inset Formula $\rec{\basis u}\equiv\left\{ \rec{\vect u}_{i}\right\} _{i=1}^{d}$
\end_inset
be the reciprocal lattice basis, that is the basis satisfying
\begin_inset Formula $\vect u_{i}\cdot\rec{\vect u_{j}}=\delta_{ij}$
\end_inset
.
This slightly differs from the usual definition of a reciprocal basis,
here denoted
\begin_inset Formula $\recb{\basis u}\equiv\left\{ \recb{\vect u_{i}}\right\} _{i=1}^{d}$
\end_inset
, which satisfies
\begin_inset Formula $\vect u_{i}\cdot\recb{\vect u_{j}}=2\pi\delta_{ij}$
\end_inset
instead.
\end_layout
\begin_layout Subsubsection
Factorisation of a multidimensional lattice delta comb
\end_layout
\begin_layout Standard
By simple drawing, it can be seen that
\begin_inset Formula
\[
\dc{\basis u}(\vect x)=c_{\basis u}\prod_{i=1}^{d}\dc{}\left(\vect x\cdot\rec{\vect u_{i}}\right)
\]
\end_inset
where
\begin_inset Formula $c_{\basis u}$
\end_inset
is some numerical volume factor.
In order to determine
\begin_inset Formula $c_{\basis u}$
\end_inset
, let us consider only the
\begin_inset Quotes eld
\end_inset
zero tooth
\begin_inset Quotes erd
\end_inset
of the comb, leading to
\begin_inset Formula
\[
\delta^{d}(\vect x)=c_{\basis u}\prod_{i=1}^{d}\delta\left(\vect x\cdot\rec{\vect u_{i}}\right).
\]
\end_inset
From the scaling property of delta function,
\begin_inset Formula $\delta(ax)=\left|a\right|^{-1}\delta(x)$
\end_inset
, we get
\begin_inset Formula
\[
\delta^{d}(\vect x)=c_{\basis u}\prod_{i=1}^{d}\left\Vert \rec{\vect u_{i}}\right\Vert ^{-1}\delta\left(\vect x\cdot\frac{\rec{\vect u_{i}}}{\left\Vert \rec{\vect u_{i}}\right\Vert }\right).
\]
\end_inset
\end_layout
\begin_layout Standard
From the Osgood's book (p.
375):
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\dc A(\vect x)=\frac{1}{\left|\det A\right|}\dc{}^{(d)}\left(A^{-1}\vect x\right)
\]
\end_inset
\begin_inset Note Note
status open
\begin_layout Plain Layout
Applying both sides to a test function that is one at the origin, we get
\begin_inset Formula $c_{\basis u}=\prod_{i=1}^{d}\left\Vert \rec{\vect u_{i}}\right\Vert $
\end_inset
SRSLY?, and hence
\begin_inset Formula
\begin{equation}
\dc{\basis u}(\vect x)=\prod_{i=1}^{d}\left\Vert \rec{\vect u_{i}}\right\Vert \dc{}\left(\vect x\cdot\rec{\vect u_{i}}\right).\label{eq:Dirac comb factorisation}
\end{equation}
\end_inset
\end_layout
\end_inset
\end_layout
\begin_layout Subsubsection
Fourier series representation
\end_layout
\begin_layout Standard
\begin_inset Note Note
status open
\begin_layout Plain Layout
Utilising the Fourier series for 1D Dirac comb
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:1D Dirac comb Fourier series"
\end_inset
and the factorisation
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:Dirac comb factorisation"
\end_inset
, we get
\begin_inset Formula
\begin{eqnarray*}
\dc{\basis u}(\vect x) & = & \prod_{j=1}^{d}\left\Vert \rec{\vect u_{j}}\right\Vert \sum_{n_{j}=-\infty}^{\infty}e^{2\pi in_{i}\vect x\cdot\rec{\vect u_{i}}}\\
& = & \left(\prod_{j=1}^{d}\left\Vert \rec{\vect u_{j}}\right\Vert \right)\sum_{\vect n\in\mathbb{Z}^{d}}e^{2\pi i\vect x\cdot\sum_{k=1}^{d}n_{k}\rec{\vect u_{k}}}.
\end{eqnarray*}
\end_inset
\end_layout
\end_inset
\end_layout
\begin_layout Subsubsection
Fourier transform (OK)
\end_layout
\begin_layout Standard
From the Osgood's book https://see.stanford.edu/materials/lsoftaee261/chap8.pdf,
p.
379
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\uoft{\dc{\basis u}}\left(\vect{\xi}\right)=\left|\det\rec{\basis u}\right|\dc{\rec{\basis u}}^{(d)}\left(\vect{\xi}\right).
\]
\end_inset
And consequently, for unitary/angular frequency it is
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{eqnarray}
\uaft{\dc{\basis u}}\left(\vect k\right) & = & \frac{1}{\left(2\pi\right)^{\frac{d}{2}}}\uoft{\dc{\basis u}}\left(\frac{\vect k}{2\pi}\right)\nonumber \\
& = & \frac{\left|\det\rec{\basis u}\right|}{\left(2\pi\right)^{\frac{d}{2}}}\dc{\rec{\basis u}}^{(d)}\left(\frac{\vect k}{2\pi}\right)\nonumber \\
& = & \left(2\pi\right)^{\frac{d}{2}}\left|\det\rec{\basis u}\right|\dc{\recb{\basis u}}\left(\vect k\right)\nonumber \\
& = & \frac{\left|\det\recb{\basis u}\right|}{\left(2\pi\right)^{\frac{d}{2}}}\dc{\recb{\basis u}}\left(\vect k\right).\label{eq:Dirac comb uaFt}
\end{eqnarray}
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Note Note
status open
\begin_layout Plain Layout
On the third line, we used the stretch theorem, getting
\begin_inset Formula
\[
\dc{\recb{\basis u}}\left(\vect k\right)=\dc{2\pi\rec{\basis u}}\left(\vect k\right)=\left(2\pi\right)^{-d}\dc{\rec{\basis u}}\left(\frac{\vect k}{2\pi}\right)
\]
\end_inset
\end_layout
\end_inset
\end_layout
\begin_layout Subsubsection
Convolution
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\left(f\ast\dc{\basis u}\right)(\vect x)=\sum_{\vect t\in\basis u\ints^{d}}f(\vect x-\vect t)
\]
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Note Note
status open
\begin_layout Plain Layout
So, from the stretch theorem
\begin_inset Formula $\uoft{(f(A\vect x))}=\frac{1}{\left|\det A\right|}\uoft{f\left(A^{-T}\vect{\xi}\right)}=\left|\det A^{-T}\right|\uoft{f\left(A^{-T}\vect{\xi}\right)}$
\end_inset
\end_layout
\begin_layout Plain Layout
From
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:Dirac comb factorisation"
\end_inset
and
\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:1D Dirac comb Ft ordinary freq"
\end_inset
\begin_inset Formula
\[
\uoft{\dc{\basis u}}(\vect{\xi})=\prod_{i=1}^{d}\left\Vert \rec{\vect u_{i}}\right\Vert \dc{}\left(\vect x\cdot\rec{\vect u_{i}}\right).
\]
\end_inset
\end_layout
\end_inset
\end_layout
\end_body
\end_document