#LyX 2.4 created this file. For more info see https://www.lyx.org/ \lyxformat 583 \begin_document \begin_header \save_transient_properties true \origin unavailable \textclass article \use_default_options true \maintain_unincluded_children false \language english \language_package default \inputencoding utf8 \fontencoding auto \font_roman "default" "TeX Gyre Pagella" \font_sans "default" "default" \font_typewriter "default" "default" \font_math "auto" "auto" \font_default_family default \use_non_tex_fonts false \font_sc false \font_roman_osf true \font_sans_osf false \font_typewriter_osf false \font_sf_scale 100 100 \font_tt_scale 100 100 \use_microtype false \use_dash_ligatures true \graphics default \default_output_format default \output_sync 0 \bibtex_command default \index_command default \float_placement class \float_alignment class \paperfontsize default \spacing single \use_hyperref true \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 default \use_geometry false \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 \use_minted 0 \use_lineno 0 \index Index \shortcut idx \color #008000 \end_index \secnumdepth 3 \tocdepth 3 \paragraph_separation indent \paragraph_indentation default \is_math_indent 0 \math_numbering_side default \quotes_style english \dynamic_quotes 0 \papercolumns 1 \papersides 1 \paperpagestyle default \tablestyle 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 Section Infinite periodic systems \begin_inset FormulaMacro \newcommand{\dlv}{\vect a} \end_inset \begin_inset FormulaMacro \newcommand{\rlv}{\vect b} \end_inset \end_layout \begin_layout Standard Although large finite systems are where MSTMM excels the most, there are several reasons that makes its extension to infinite lattices (where periodic boundary conditions might be applied) desirable as well. Other methods might be already fast enough, but MSTMM will be faster in most cases in which there is enough spacing between the neighboring particles. MSTMM works well with any space group symmetry the system might have (as opposed to, for example, FDTD with cubic mesh applied to a honeycomb lattice), which makes e.g. application of group theory in mode analysis quite easy. \begin_inset Note Note status open \begin_layout Plain Layout Topology anoyne? \end_layout \end_inset And finally, having a method that handles well both infinite and large finite system gives a possibility to study finite-size effects in periodic scatterer arrays. \end_layout \begin_layout Subsection Notation \end_layout \begin_layout Standard TODO Fourier transforms, Delta comb, lattice bases, reciprocal lattices etc. \end_layout \begin_layout Subsection Formulation of the problem \end_layout \begin_layout Standard Let us have a linear system of compact EM scatterers on a homogeneous background as in Section \begin_inset CommandInset ref LatexCommand eqref reference "subsec:Multiple-scattering" plural "false" caps "false" noprefix "false" \end_inset , but this time, the system shall be periodic: let there be a \begin_inset Formula $d$ \end_inset -dimensional ( \begin_inset Formula $d$ \end_inset can be 1, 2 or 3) lattice embedded into the three-dimensional real space, with lattice vectors \begin_inset Formula $\left\{ \dlv_{i}\right\} _{i=1}^{d}$ \end_inset , and let the lattice points be labeled with an \begin_inset Formula $d$ \end_inset -dimensional integar multiindex \begin_inset Formula $\vect n\in\ints^{d}$ \end_inset , so the lattice points have cartesian coordinates \begin_inset Formula $\vect R_{\vect n}=\sum_{i=1}^{d}n_{i}\vect a_{i}$ \end_inset . There can be several scatterers per unit cell with indices \begin_inset Formula $\alpha$ \end_inset from set \begin_inset Formula $\mathcal{P}_{1}$ \end_inset and (relative) positions inside the unit cell \begin_inset Formula $\vect r_{\alpha}$ \end_inset ; any particle of the periodic system can thus be labeled by a multiindex from \begin_inset Formula $\mathcal{P}=\ints^{d}\times\mathcal{P}_{1}$ \end_inset . The scatterers are located at positions \begin_inset Formula $\vect r_{\vect n,\alpha}=\vect R_{\vect n}+\vect r_{\alpha}$ \end_inset and their \begin_inset Formula $T$ \end_inset -matrices are periodic, \begin_inset Formula $T_{\vect n,\alpha}=T_{\alpha}$ \end_inset . In such system, the multiple-scattering problem \begin_inset CommandInset ref LatexCommand eqref reference "eq:Multiple-scattering problem" plural "false" caps "false" noprefix "false" \end_inset can be rewritten as \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} \outcoeffp{\vect n,\alpha}-T_{\alpha}\sum_{\left(\vect m,\beta\right)\in\mathcal{P}\backslash\left\{ \left(\vect n,\alpha\right)\right\} }\tropsp{\vect n,\alpha}{\vect m,\beta}\outcoeffp{\vect m,\beta}=T_{\alpha}\rcoeffincp{\vect n,\alpha}.\quad\left(\vect n,\alpha\right)\in\mathcal{P}\label{eq:Multiple-scattering problem periodic} \end{equation} \end_inset \end_layout \begin_layout Standard Due to periodicity, we can also write \begin_inset Formula $\tropsp{\vect n,\alpha}{\vect m,\beta}=\tropsp{\alpha}{\beta}\left(\vect R_{\vect m}-\vect R_{\vect n}\right)=\tropsp{\alpha}{\beta}\left(\vect R_{\vect m-\vect n}\right)=\tropsp{\vect 0,\alpha}{\vect m-\vect n,\beta}$ \end_inset . Assuming quasi-periodic right-hand side with quasi-momentum \begin_inset Formula $\vect k$ \end_inset , \begin_inset Formula $\rcoeffincp{\vect n,\alpha}=\rcoeffincp{\vect 0,\alpha}\left(\vect k\right)e^{i\vect k\cdot\vect R_{\vect n}}$ \end_inset , the solutions of \begin_inset CommandInset ref LatexCommand eqref reference "eq:Multiple-scattering problem periodic" plural "false" caps "false" noprefix "false" \end_inset will be also quasi-periodic according to Bloch theorem, \begin_inset Formula $\outcoeffp{\vect n,\alpha}=\outcoeffp{\vect 0,\alpha}\left(\vect k\right)e^{i\vect k\cdot\vect R_{\vect n}}$ \end_inset , and eq. \begin_inset CommandInset ref LatexCommand eqref reference "eq:Multiple-scattering problem periodic" plural "false" caps "false" noprefix "false" \end_inset can be rewritten as follows \begin_inset Formula \begin{align} \outcoeffp{\vect 0,\alpha}\left(\vect k\right)e^{i\vect k\cdot\vect R_{\vect n}}-T_{\alpha}\sum_{\left(\vect m,\beta\right)\in\mathcal{P}\backslash\left\{ \left(\vect n,\alpha\right)\right\} }\tropsp{\vect n,\alpha}{\vect m,\beta}\outcoeffp{\vect 0,\beta}\left(\vect k\right)e^{i\vect k\cdot\vect R_{\vect m}} & =T_{\alpha}\rcoeffincp{\vect 0,\alpha}\left(\vect k\right)e^{i\vect k\cdot\vect R_{\vect n}},\nonumber \\ \outcoeffp{\vect 0,\alpha}\left(\vect k\right)-T_{\alpha}\sum_{\left(\vect m,\beta\right)\in\mathcal{P}\backslash\left\{ \left(\vect n,\alpha\right)\right\} }\tropsp{\vect 0,\alpha}{\vect m-\vect n,\beta}\outcoeffp{\vect 0,\beta}\left(\vect k\right)e^{i\vect k\cdot\vect R_{\vect m-\vect n}} & =T_{\alpha}\rcoeffincp{\vect 0,\alpha}\left(\vect k\right),\nonumber \\ \outcoeffp{\vect 0,\alpha}\left(\vect k\right)-T_{\alpha}\sum_{\left(\vect m,\beta\right)\in\mathcal{P}\backslash\left\{ \left(\vect 0,\alpha\right)\right\} }\tropsp{\vect 0,\alpha}{\vect m,\beta}\outcoeffp{\vect 0,\beta}\left(\vect k\right)e^{i\vect k\cdot\vect R_{\vect m}} & =T_{\alpha}\rcoeffincp{\vect 0,\alpha}\left(\vect k\right),\nonumber \\ \outcoeffp{\vect 0,\alpha}\left(\vect k\right)-T_{\alpha}\sum_{\beta\in\mathcal{P}}W_{\alpha\beta}\left(\vect k\right)\outcoeffp{\vect 0,\beta}\left(\vect k\right) & =T_{\alpha}\rcoeffincp{\vect 0,\alpha}\left(\vect k\right),\label{eq:Multiple-scattering problem unit cell} \end{align} \end_inset so we reduced the initial scattering problem to one involving only the field expansion coefficients from a single unit cell, but 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 m\in\ints^{d}}\left(1-\delta_{\alpha\beta}\right)\tropsp{\vect 0,\alpha}{\vect m,\beta}e^{i\vect k\cdot\vect R_{\vect m}},\label{eq:W definition} \end{equation} \end_inset evaluation of which is possible but quite non-trivial due to the infinite lattice sum, so we explain it separately in Sect. \begin_inset CommandInset ref LatexCommand eqref reference "subsec:W operator evaluation" plural "false" caps "false" noprefix "false" \end_inset . \end_layout \begin_layout Standard As in the case of a finite system, eq. can be written in a shorter block-matrix form, \begin_inset Formula \begin{equation} \left(I-WT\right)\outcoeffp{\vect 0}\left(\vect k\right)=\rcoeffincp{\vect 0}\left(\vect k\right)\label{eq:Multiple-scattering problem unit cell block form} \end{equation} \end_inset Eq. \begin_inset CommandInset ref LatexCommand eqref reference "eq:Multiple-scattering problem unit cell" plural "false" caps "false" noprefix "false" \end_inset can be used to calculate electromagnetic response of the structure to external quasiperiodic driving field – most notably a plane wave. However, the non-trivial solutions of the equation with right hand side (i.e. the external driving) set to zero, \begin_inset Formula \begin{equation} \left(I-WT\right)\outcoeffp{\vect 0}\left(\vect k\right)=0,\label{eq:lattice mode equation} \end{equation} \end_inset describes the \emph on lattice modes. \emph default Non-trivial solutions to \begin_inset CommandInset ref LatexCommand eqref reference "eq:lattice mode equation" plural "false" caps "false" noprefix "false" \end_inset exist if the matrix on the left-hand side \begin_inset Formula $M\left(\omega,\vect k\right)=\left(I-W\left(\omega,\vect k\right)T\left(\omega\right)\right)$ \end_inset is singular – this condition gives the \emph on dispersion relation \emph default for the periodic structure. Note that in realistic (lossy) systems, at least one of the pair \begin_inset Formula $\omega,\vect k$ \end_inset will acquire complex values. The solution \begin_inset Formula $\outcoeffp{\vect 0}\left(\vect k\right)$ \end_inset is then obtained as the right \begin_inset Note Note status open \begin_layout Plain Layout CHECK! \end_layout \end_inset singular vector of \begin_inset Formula $M\left(\omega,\vect k\right)$ \end_inset corresponding to the zero singular value. \end_layout \begin_layout Subsection Numerical solution \end_layout \begin_layout Standard In practice, equation \begin_inset CommandInset ref LatexCommand eqref reference "eq:Multiple-scattering problem unit cell block form" plural "false" caps "false" noprefix "false" \end_inset is solved in the same way as eq. \begin_inset CommandInset ref LatexCommand eqref reference "eq:Multiple-scattering problem block form" plural "false" caps "false" noprefix "false" \end_inset in the multipole degree truncated form. \end_layout \begin_layout Standard The lattice mode problem \begin_inset CommandInset ref LatexCommand eqref reference "eq:lattice mode equation" plural "false" caps "false" noprefix "false" \end_inset is (after multipole degree truncation) solved by finding \begin_inset Formula $\omega,\vect k$ \end_inset for which the matrix \begin_inset Formula $M\left(\omega,\vect k\right)$ \end_inset has a zero singular value. A naïve approach to do that is to sample a volume with a grid in the \begin_inset Formula $\left(\omega,\vect k\right)$ \end_inset space, performing a singular value decomposition of \begin_inset Formula $M\left(\omega,\vect k\right)$ \end_inset at each point and finding where the lowest singular value of \begin_inset Formula $M\left(\omega,\vect k\right)$ \end_inset is close enough to zero. However, this approach is quite expensive, for \begin_inset Formula $W\left(\omega,\vect k\right)$ \end_inset has to be evaluated for each \begin_inset Formula $\omega,\vect k$ \end_inset pair separately (unlike the original finite case \begin_inset CommandInset ref LatexCommand eqref reference "eq:Multiple-scattering problem block form" plural "false" caps "false" noprefix "false" \end_inset translation operator \begin_inset Formula $\trops$ \end_inset , which, for a given geometry, depends only on frequency). Therefore, a much more efficient approach to determine the photonic bands is to sample the \begin_inset Formula $\vect k$ \end_inset -space (a whole Brillouin zone or its part) and for each fixed \begin_inset Formula $\vect k$ \end_inset to find a corresponding frequency \begin_inset Formula $\omega$ \end_inset with zero singular value of \begin_inset Formula $M\left(\omega,\vect k\right)$ \end_inset using a minimisation algorithm (two- or one-dimensional, depending on whether one needs the exact complex-valued \begin_inset Formula $\omega$ \end_inset or whether the its real-valued approximation is satisfactory). Typically, a good initial guess for \begin_inset Formula $\omega\left(\vect k\right)$ \end_inset is obtained from the empty lattice approximation, \begin_inset Formula $\left|\vect k\right|=\sqrt{\epsilon\mu}\omega/c_{0}$ \end_inset (modulo lattice points; TODO write this a clean way). A somehow challenging step is to distinguish the different bands that can all be very close to the empty lattice approximation, especially if the particles in the systems are small. In high-symmetry points of the Brilloin zone, this can be solved by factorising \begin_inset Formula $M\left(\omega,\vect k\right)$ \end_inset into irreducible representations \begin_inset Formula $\Gamma_{i}$ \end_inset and performing the minimisation in each irrep separately, cf. Section \begin_inset CommandInset ref LatexCommand ref reference "sec:Symmetries" plural "false" caps "false" noprefix "false" \end_inset , and using the different \begin_inset Formula $\omega_{\Gamma_{i}}\left(\vect k\right)$ \end_inset to obtain the initial guesses for the nearby points \begin_inset Formula $\vect k+\delta\vect k$ \end_inset . \end_layout \begin_layout Subsection Computing the Fourier sum of the translation operator \begin_inset CommandInset label LatexCommand label name "subsec:W operator evaluation" \end_inset \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 $\tropsp{\vect 0,\alpha}{\vect m,\beta}\sim\left|\vect R_{\vect b}\right|^{-1}e^{ik_{0}\left|\vect R_{\vect b}\right|}$ \end_inset that does not in the strict sense converge for any \begin_inset Formula $d>1$ \end_inset -dimensional lattice. \begin_inset Note Note status open \begin_layout Plain Layout \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 \end_layout \end_inset In electrostatics, this problem can be solved with Ewald summation [TODO REF]. Its basic idea is that if what asymptoticaly decays poorly in the direct space, will perhaps decay fast in the Fourier space. We use the same idea here, but the technical details are more complicated 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 \begin_inset FormulaMacro \renewcommand{\basis}[1]{\mathfrak{#1}} \end_inset \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}\\ & = & \frac{\left|\det\rec{\basis u}\right|}{\left(2\pi\right)^{\frac{d}{2}}}\sum_{\vect K\in\recb{\basis u}\ints^{d}}e^{i\left(\vect k-\vect K\right)\cdot\left(-\vect r_{\beta}+\vect r_{\alpha}\right)}\left(\uaft{S(\vect{\bullet}\leftarrow\vect 0)}\right)\left(\vect k-\vect K\right)\nonumber \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 expected to converge nicely. We note that the elements of the translation operators \begin_inset Formula $\tropr,\trops$ \end_inset in \begin_inset CommandInset ref LatexCommand eqref reference "eq:translation operator" plural "false" caps "false" noprefix "false" \end_inset can be rewritten as linear combinations of expressions \begin_inset Formula $\ush{\nu}{\mu}\left(\uvec d\right)j_{n}\left(d\right),\ush{\nu}{\mu}\left(\uvec d\right)h_{n}^{(1)}\left(d\right)$ \end_inset (TODO WRITE THEM EXPLICITLY IN THIS FORM), respectively, hence if we are able evaluate the lattice sums sums \begin_inset Note Note status open \begin_layout Plain Layout CHECK THE FOLLOWING EXPRESSION FOR CORRECT FUNCTION ARGUMENTS \end_layout \end_inset \begin_inset Formula \begin{equation} \sigma_{\nu}^{\mu}\left(\vect k\right)=\sum_{\vect n\in\ints^{d}\backslash\left\{ \vect 0\right\} }e^{i\vect{\vect k}\cdot\vect R_{\vect n}}\ush{\nu}{\mu}\left(\uvec{R_{n}}\right)h_{n}^{(1)}\left(R_{n}\right),\label{eq:sigma lattice sums} \end{equation} \end_inset then by linearity, we can get the \begin_inset Formula $W_{\alpha\beta}\left(\vect k\right)$ \end_inset operator as well. \end_layout \begin_layout Standard TODO ADD MOROZ AND OTHER REFS HERE. \begin_inset CommandInset citation LatexCommand cite key "linton_one-_2009" literal "true" \end_inset offers an exponentially convergent Ewald-type summation method for \begin_inset Formula $\sigma_{\nu}^{\mu}\left(\vect k\right)=\sigma_{\nu}^{\mu(\mathrm{S})}\left(\vect k\right)+\sigma_{\nu}^{\mu(\mathrm{L})}\left(\vect k\right)$ \end_inset . Here we rewrite them in a form independent on the convention used for spherical harmonics (as long as they are complex \begin_inset Note Note status open \begin_layout Plain Layout lepší formulace \end_layout \end_inset ). The short-range part reads (UNIFY INDEX NOTATION) \begin_inset Formula \begin{multline} \sigma_{n}^{m(\mathrm{S})}\left(\vect{\beta}\right)=-\frac{2^{n+1}i}{k^{n+1}\sqrt{\pi}}\sum_{\vect R\in\Lambda}^{'}\left|\vect R\right|^{n}e^{i\vect{\beta}\cdot\vect R}Y_{n}^{m}\left(\vect R\right)\int_{\eta}^{\infty}e^{-\left|\vect R\right|^{2}\xi^{2}}e^{-k/4\xi^{2}}\xi^{2n}\ud\xi\\ +\frac{\delta_{n0}\delta_{m0}}{\sqrt{4\pi}}\Gamma\left(-\frac{1}{2},-\frac{k}{4\eta^{2}}\right)Y_{n}^{m},\label{eq:Ewald in 3D short-range part} \end{multline} \end_inset \begin_inset Note Note status open \begin_layout Plain Layout NEPATŘÍ TAM NĚJAKÁ DELTA FUNKCE K PŮVODNÍMU \begin_inset Formula $\sigma_{n}^{m(0)}$ \end_inset ? \end_layout \end_inset and the long-range part (FIXME, this is the 2D version; include the 1D and 3D lattice expressions as well) \begin_inset Formula \begin{multline} \sigma_{n}^{m(\mathrm{L})}\left(\vect{\beta}\right)=-\frac{i^{n+1}}{k^{2}\mathscr{A}}\sqrt{\pi}2\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!\times\\ \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(\frac{\pi}{2},\phi_{\vect{\beta}_{pq}}\right)\sum_{j=0}^{\left[\left(n-\left|m\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/k\right)^{n-2j}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(n-m\right)-j\right)!\left(\frac{1}{2}\left(n+m\right)-j\right)!}\left(\gamma_{pq}\right)^{2j-1}\label{eq:Ewald in 3D long-range part} \end{multline} \end_inset where \begin_inset Formula $\xi$ \end_inset is TODO, \begin_inset Formula $\beta_{pq}$ \end_inset is TODO, \begin_inset Formula $\Gamma_{j,pq}$ \end_inset is TODO and \begin_inset Formula $\eta$ \end_inset is a real parameter that determines the pace of convergence of both parts. The larger \begin_inset Formula $\eta$ \end_inset is, the faster \begin_inset Formula $\sigma_{n}^{m(\mathrm{S})}\left(\vect{\beta}\right)$ \end_inset converges but the slower \begin_inset Formula $\sigma_{n}^{m(\mathrm{L})}\left(\vect{\beta}\right)$ \end_inset converges. Therefore (based on the lattice geometry) it has to be adjusted in a way that a reasonable amount of terms needs to be evaluated numerically from both \begin_inset Formula $\sigma_{n}^{m(\mathrm{S})}\left(\vect{\beta}\right)$ \end_inset and \begin_inset Formula $\sigma_{n}^{m(\mathrm{L})}\left(\vect{\beta}\right)$ \end_inset . Generally, a good choice for \begin_inset Formula $\eta$ \end_inset is TODO; in order to achieve accuracy TODO, one has to evaluate the terms on TODO lattice points. (I HAVE SOME DERIVATIONS OF THE ESTIMATES IN MY NOTES; SHOULD I INCLUDE THEM?) \end_layout \begin_layout Standard In practice, the integrals in \begin_inset CommandInset ref LatexCommand eqref reference "eq:Ewald in 3D short-range part" plural "false" caps "false" noprefix "false" \end_inset can be easily evaluated by numerical quadrature and the incomplete \begin_inset Formula $\Gamma$ \end_inset -functions using the series TODO and TODO from DLMF. \end_layout \end_body \end_document