diff --git a/lepaper/arrayscat.lyx b/lepaper/arrayscat.lyx
index a7548ab..c318f56 100644
--- a/lepaper/arrayscat.lyx
+++ b/lepaper/arrayscat.lyx
@@ -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
diff --git a/lepaper/finite.lyx b/lepaper/finite.lyx
index 8a63351..97c63f1 100644
--- a/lepaper/finite.lyx
+++ b/lepaper/finite.lyx
@@ -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
diff --git a/lepaper/infinite.lyx b/lepaper/infinite.lyx
index 141050b..296ba20 100644
--- a/lepaper/infinite.lyx
+++ b/lepaper/infinite.lyx
@@ -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