Infinite systems WIP (how to numerical solutions)

Former-commit-id: dd3118e392b58eb9ae2260581e93028561571245
This commit is contained in:
Marek Nečada 2019-07-31 10:10:45 +03:00
parent d5b6f8f5d4
commit f167360c1e
3 changed files with 233 additions and 24 deletions

View File

@ -802,6 +802,17 @@ literal "true"
\end_inset
\end_layout
\begin_layout Standard
\begin_inset CommandInset include
LatexCommand include
filename "symmetries.lyx"
literal "true"
\end_inset
\end_layout
\begin_layout Standard

View File

@ -290,8 +290,8 @@ outgoing
, respectively, defined as follows:
\begin_inset Formula
\begin{align*}
\vswfrtlm 1lm\left(k\vect r\right) & =j_{l}\left(kr\right)\vsh 1lm\left(\uvec r\right),\\
\vswfrtlm 2lm\left(k\vect r\right) & =\frac{1}{kr}\frac{\ud\left(krj_{l}\left(kr\right)\right)}{\ud\left(kr\right)}\vsh 2lm\left(\uvec r\right)+\sqrt{l\left(l+1\right)}\frac{j_{l}\left(kr\right)}{kr}\vsh 3lm\left(\uvec r\right),
\vswfrtlm1lm\left(k\vect r\right) & =j_{l}\left(kr\right)\vsh1lm\left(\uvec r\right),\\
\vswfrtlm2lm\left(k\vect r\right) & =\frac{1}{kr}\frac{\ud\left(krj_{l}\left(kr\right)\right)}{\ud\left(kr\right)}\vsh2lm\left(\uvec r\right)+\sqrt{l\left(l+1\right)}\frac{j_{l}\left(kr\right)}{kr}\vsh3lm\left(\uvec r\right),
\end{align*}
\end_inset
@ -299,8 +299,8 @@ outgoing
\begin_inset Formula
\begin{align*}
\vswfouttlm 1lm\left(k\vect r\right) & =h_{l}^{\left(1\right)}\left(kr\right)\vsh 1lm\left(\uvec r\right),\\
\vswfouttlm 2lm\left(k\vect r\right) & =\frac{1}{kr}\frac{\ud\left(krh_{l}^{\left(1\right)}\left(kr\right)\right)}{\ud\left(kr\right)}\vsh 2lm\left(\uvec r\right)+\sqrt{l\left(l+1\right)}\frac{h_{l}^{\left(1\right)}\left(kr\right)}{kr}\vsh 3lm\left(\uvec r\right),\\
\vswfouttlm1lm\left(k\vect r\right) & =h_{l}^{\left(1\right)}\left(kr\right)\vsh1lm\left(\uvec r\right),\\
\vswfouttlm2lm\left(k\vect r\right) & =\frac{1}{kr}\frac{\ud\left(krh_{l}^{\left(1\right)}\left(kr\right)\right)}{\ud\left(kr\right)}\vsh2lm\left(\uvec r\right)+\sqrt{l\left(l+1\right)}\frac{h_{l}^{\left(1\right)}\left(kr\right)}{kr}\vsh3lm\left(\uvec r\right),\\
& \tau=1,2;\quad l=1,2,3,\dots;\quad m=-l,-l+1,\dots,+l,
\end{align*}
@ -326,9 +326,9 @@ vector spherical harmonics
\begin_inset Formula
\begin{align*}
\vsh 1lm\left(\uvec r\right) & =\frac{1}{\sqrt{l\left(l+1\right)}}\nabla\times\left(\vect r\ush lm\left(\uvec r\right)\right)=\frac{1}{\sqrt{l\left(l+1\right)}}\nabla\ush lm\left(\uvec r\right)\times\vect r,\\
\vsh 2lm\left(\uvec r\right) & =\frac{1}{\sqrt{l\left(l+1\right)}}r\nabla\ush lm\left(\uvec r\right),\\
\vsh 3lm\left(\uvec r\right) & =\uvec r\ush lm\left(\uvec r\right).
\vsh1lm\left(\uvec r\right) & =\frac{1}{\sqrt{l\left(l+1\right)}}\nabla\times\left(\vect r\ush lm\left(\uvec r\right)\right)=\frac{1}{\sqrt{l\left(l+1\right)}}\nabla\ush lm\left(\uvec r\right)\times\vect r,\\
\vsh2lm\left(\uvec r\right) & =\frac{1}{\sqrt{l\left(l+1\right)}}r\nabla\ush lm\left(\uvec r\right),\\
\vsh3lm\left(\uvec r\right) & =\uvec r\ush lm\left(\uvec r\right).
\end{align*}
\end_inset
@ -452,7 +452,7 @@ noprefix "false"
\end_inset
inside a ball
\begin_inset Formula $\openball 0{R^{>}}$
\begin_inset Formula $\openball0{R^{>}}$
\end_inset
with radius
@ -470,7 +470,7 @@ noprefix "false"
\end_inset
to have a complete basis of the solutions in the volume
\begin_inset Formula $\openball 0{R^{>}}\backslash B_{0}\left(R\right)$
\begin_inset Formula $\openball0{R^{>}}\backslash B_{0}\left(R\right)$
\end_inset
.
@ -496,7 +496,7 @@ The single-particle scattering problem at frequency
\end_inset
and let the whole volume
\begin_inset Formula $\openball 0{R^{>}}\backslash B_{0}\left(R\right)$
\begin_inset Formula $\openball0{R^{>}}\backslash B_{0}\left(R\right)$
\end_inset
be filled with a homogeneous isotropic medium with wave number
@ -532,7 +532,7 @@ If there was no scatterer and
\end_inset
due to sources outside
\begin_inset Formula $\openball 0R$
\begin_inset Formula $\openball0R$
\end_inset
would remain.
@ -793,7 +793,7 @@ literal "true"
.
Let the field in
\begin_inset Formula $\openball 0{R^{>}}\backslash B_{0}\left(R\right)$
\begin_inset Formula $\openball0{R^{>}}\backslash B_{0}\left(R\right)$
\end_inset
have expansion as in
@ -812,7 +812,7 @@ noprefix "false"
\end_inset
to
\begin_inset Formula $\openball 0{R^{>}}\backslash B_{0}\left(R\right)$
\begin_inset Formula $\openball0{R^{>}}\backslash B_{0}\left(R\right)$
\end_inset
via by electromagnetic radiation is
@ -864,8 +864,8 @@ In many scattering problems considered in practice, the driving field is
with expansion coefficients
\begin_inset Formula
\begin{eqnarray}
\rcoeffptlm{}1lm\left(\vect k,\vect E_{0}\right) & = & 4\pi i^{l}\vshD 1lm\left(\uvec k\right),\nonumber \\
\rcoeffptlm{}2lm\left(\vect k,\vect E_{0}\right) & = & -4\pi i^{l+1}\vshD 2lm\left(\uvec k\right).\label{eq:plane wave expansion}
\rcoeffptlm{}1lm\left(\vect k,\vect E_{0}\right) & = & 4\pi i^{l}\vshD1lm\left(\uvec k\right),\nonumber \\
\rcoeffptlm{}2lm\left(\vect k,\vect E_{0}\right) & = & -4\pi i^{l+1}\vshD2lm\left(\uvec k\right).\label{eq:plane wave expansion}
\end{eqnarray}
\end_inset
@ -1251,7 +1251,17 @@ noprefix "false"
\end_inset
can then be solved using standard numerical linear algebra methods.
can then be solved using standard numerical linear algebra methods (typically,
by LU factorisation of the
\begin_inset Formula $\left(I-T\trops\right)$
\end_inset
matrix at a given frequency, and then solving with Gauss elimination for
as many different incident
\begin_inset Formula $\rcoeffinc$
\end_inset
vectors as needed).
\end_layout
\begin_layout Standard

View File

@ -324,9 +324,9 @@ noprefix "false"
As in the case of a finite system, eq.
can be written in a shorter block-matrix form,
\begin_inset Formula
\[
\left(I-WT\right)\outcoeffp{\vect 0}\left(\vect k\right)=\rcoeffincp{\vect 0}\left(\vect k\right)
\]
\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
@ -343,8 +343,196 @@ noprefix "false"
can be used to calculate electromagnetic response of the structure to external
quasiperiodic driving field most notably a plane wave.
However, if one sets the right the right-hand side to zero, it can also
be used to find electromagnetic lattice modes
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.
\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
@ -617,8 +805,8 @@ reference "eq:W sum in reciprocal space"
\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}}(\vect0\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\vect0)}\right)\left(\vect k-\vect K\right)\label{eq:W Long definition}
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
@ -656,7 +844,7 @@ CHECK THE FOLLOWING EXPRESSION FOR CORRECT FUNCTION ARGUMENTS
\begin_inset Formula
\begin{equation}
\sigma_{\nu}^{\mu}\left(\vect k\right)=\sum_{\vect n\in\ints^{d}\backslash\left\{ \vect0\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}
\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