\documentclass{article}
\usepackage[top=30pt,bottom=30pt,left=48pt,right=46pt]{geometry}

\usepackage[utf8]{inputenc}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{pdflscape}

\newcommand{\BG}{\bar{g}}
\newcommand{\BKN}{\bar{k}_n}

\newcommand{\CT}{\cos{(\theta^-)}}
\newcommand{\CTB}{\cos{(\Bar{\theta}^-)}}
\newcommand{\CTT}{\cos{(2\theta^-)}}
\newcommand{\CST}{\cos^2{(\theta^-)}}
\newcommand{\CSTB}{\cos^2{(\Bar{\theta}^-)}}
\newcommand{\ST}{\sin{(\theta^-)}}
\newcommand{\STB}{\sin{(\Bar{\theta}^-)}}
\newcommand{\STT}{\sin{(2\theta^-)}}
\newcommand{\SST}{\sin^2{(\theta^-)}}
\newcommand{\SSTB}{\sin^2{(\Bar{\theta}^-)}}


\newcommand{\CTH}{\cos{(\theta)}} 
\newcommand{\CTPH}{\cos{(\phi})}
\newcommand{\STH}{\sin{(\theta)}} 
\newcommand{\STPH}{\sin{(\phi})}



\newcommand{\CTHPM}{\cos{(\theta^\pm)}} 
\newcommand{\STHPM}{\sin{(\theta^\pm)}} 

\newcommand{\CTP}{\cos{(\theta^+)}}
\newcommand{\CTPB}{\cos{(\Bar{\theta}^+)}}
\newcommand{\CTTP}{\cos{(2\theta^+)}}
\newcommand{\CSTP}{\cos^2{(\theta^+)}}
\newcommand{\CSTPB}{\cos^2{(\Bar{\theta}^+)}}
\newcommand{\STP}{\sin{(\theta^+)}}
\newcommand{\STPB}{\sin{(\Bar{\theta}^+)}}
\newcommand{\STTP}{\sin{(2\theta^+)}}
\newcommand{\SSTP}{\sin^2{(\theta^+)}}
\newcommand{\SSTPB}{\sin^2{(\Bar{\theta}^+)}}

\newcommand{\CP}{\cos{(\phi)}}
\newcommand{\CSPH}{\cos^2{(\phi)}}
\newcommand{\SP}{\sin{(\phi)}}
\newcommand{\SSPH}{\sin^2{(\phi)}}

\newcommand{\FXP}{F_{x'}^\rightharpoonup}
\newcommand{\FYP}{F_{y'}^\rightharpoonup}
\newcommand{\FXM}{F_{x'}^\leftharpoonup}
\newcommand{\FYM}{F_{y'}^\leftharpoonup}

\newcommand{\MXPXP}{\tilde{M}_{x'+}^{x'+}}
\newcommand{\MXPYP}{\tilde{M}_{y'+}^{x'+}}
\newcommand{\MXPXM}{\tilde{M}_{x'-}^{x'+}}
\newcommand{\MXPYM}{\tilde{M}_{y'-}^{x'+}}
\newcommand{\MYPXP}{\tilde{M}_{x'+}^{y'+}}
\newcommand{\MYPYP}{\tilde{M}_{y'+}^{y'+}}
\newcommand{\MYPXM}{\tilde{M}_{x'-}^{y'+}}
\newcommand{\MYPYM}{\tilde{M}_{y'-}^{y'+}}
\newcommand{\MXMXP}{\tilde{M}_{x'+}^{x'-}}
\newcommand{\MXMYP}{\tilde{M}_{y'+}^{x'-}}
\newcommand{\MXMXM}{\tilde{M}_{x'-}^{x'-}}
\newcommand{\MXMYM}{\tilde{M}_{y'-}^{x'-}}
\newcommand{\MYMXP}{\tilde{M}_{x'+}^{y'-}}
\newcommand{\MYMYP}{\tilde{M}_{y'+}^{y'-}}
\newcommand{\MYMXM}{\tilde{M}_{x'-}^{y'-}}
\newcommand{\MYMYM}{\tilde{M}_{y'-}^{y'-}}



\newcommand{\ALE}{2\alpha_\varepsilon}
\newcommand{\ALMU}{2\alpha_\mu}
\newcommand{\ALED}{\alpha_\varepsilon}
\newcommand{\ALMUD}{\alpha_\mu}
\newcommand{\ALM}{\alpha^-}
\newcommand{\ALP}{\alpha^+}
\newcommand{\ALMB}{\bar{\alpha}^-}
\newcommand{\ALPB}{\bar{\alpha}^+}
\newcommand{\ALEB}{2\Bar{\alpha}_\varepsilon}
\newcommand{\ALMUB}{2\Bar{\alpha}_\mu}

\newcommand{\OM}{\Omega}
\newcommand{\EB}{\varepsilon_1}
\newcommand{\MB}{\mu_1}
\newcommand{\EYN}{E_{y,n}}
\newcommand{\EXN}{E_{x,n}}
\newcommand{\HYN}{H_{y,n}}
\newcommand{\HXN}{H_{x,n}}
\newcommand{\KN}{(k+n g)}
\newcommand{\OMN}{(\omega+n\Omega)}
\newcommand{\PDT}[1]{\frac{\partial #1}{\partial t}}
\newcommand{\PDX}[1]{\frac{\partial #1}{\partial x}}
\newcommand{\PDY}[1]{\frac{\partial #1}{\partial y}}
\newcommand{\PDZ}[1]{\frac{\partial #1}{\partial z}}
\newcommand{\SQBR}[1]{\bigg[ #1 \bigg]}
\newcommand{\RBR}[1]{\bigg( #1 \bigg)}
\newcommand{\PXP}{P_x^+}
\newcommand{\PXM}{P_x^-}
\newcommand{\PYP}{P_y^+}
\newcommand{\PYM}{P_y^-}

\newcommand{\PPP}{P_+^+}
\newcommand{\PPM}{P_-^+}
\newcommand{\PMP}{P_+^-}
\newcommand{\PMM}{P_-^-}

\renewcommand{\cos}[1]{c #1 }
\renewcommand{\sin}[1]{s #1 }

\title{An Archimedes' Screw for Light - Supplementary Material}
\author{Emanuele Galiffi, Paloma A. Huidobro and J. B. Pendry}
% \date{}

\begin{document}

\maketitle





%%%
%%%

%Having discussed electromagnetic media featuring a linear, isotropic space-time modulation of their dielectric and magnetic parameters, we now put a twist on them, moving on to explore the new realm of chiral space-time metamaterials. In this chapter we present a rich analytic model for what constitutes effectively an \textit{Archimedes' screw for light}, and demonstrate its capabilities for the drag and amplification of light, as well as its potential as a model to study high order exceptional points and topology in electromagnetics.

\section{Analytic solution for $D$ and $B$ for $\omega(k)$}

Let us consider a medium with a uniaxial perturbation of its electromagnetic tensors, which rotates describing a helix in space and time. The tensors are subject the following coordinate transformation:
\begin{align}
    \varepsilon/\varepsilon_1 &= 
    \begin{pmatrix}
    1&0&0\\
    0&1&0\\
    0&0&1
    \end{pmatrix} +
    \begin{pmatrix}
    \CT&-\ST&0\\
    \ST&\CT&0\\
    0&0&1
    \end{pmatrix}
    \begin{pmatrix}
    \ALE&0&0\\
    0&0&0\\
    0&0&0
    \end{pmatrix}
    \begin{pmatrix}
    \CT&\ST&0\\
    -\ST&\CT&0\\
    0&0&1
    \end{pmatrix},
\end{align}{}
and similarly for $\mu$:
\begin{align}
    \mu/\mu_1 &= 
    \begin{pmatrix}
    1&0&0\\
    0&1&0\\
    0&0&1
    \end{pmatrix} +
    \begin{pmatrix}
    \CTP&-\STP&0\\
    \STP&\CTP&0\\
    0&0&1
    \end{pmatrix}
    \begin{pmatrix}
    0&0&0\\
    0&\ALMU&0\\
    0&0&0
    \end{pmatrix}
    \begin{pmatrix}
    \CTP&\STP&0\\
    -\STP&\CTP&0\\
    0&0&1
    \end{pmatrix},
\end{align}{}
where $\theta^\pm = g z - \OM t \pm \phi$ and $\epsilon_1$ and $\mu_1$ are the background permittivity and permeability, and we denote, for brevity, cosines and sines as $c$ and $s$ respectively. The constant $\phi$ results in a dephasing of $2\phi$ between the electric and the magnetic components of the modulation. Note that the modulations of the $x-x$ component of $\varepsilon$ and the $y-y$ of $\mu$ coincide, so that the modulation acts on a specific linear polarisation wherever at a specific point in spacetime (e.g. for $\theta = \phi = 0$, only $x$-polarised waves experience a different refractive index, while $y$-polarised waves are unaffected). Performing the matrix multiplications, we obtain the following forms for the modulated tensors:
\begin{align}
    \varepsilon/\varepsilon_1 &= 
    \begin{pmatrix}
    1&0&0\\
    0&1&0\\
    0&0&1
    \end{pmatrix} + \ALE
    \begin{pmatrix}
    \CST&\CT\ST&0\\
    \CT\ST&\SST&0\\
    0&0&0
    \end{pmatrix}\\
    \mu/\mu_1 &=     
    \begin{pmatrix}
    1&0&0\\
    0&1&0\\
    0&0&1
    \end{pmatrix} + \ALMU
    \begin{pmatrix}
    \SSTP&-\CTP\STP&0\\
    -\CTP\STP&\CSTP&0\\
    0&0&0
    \end{pmatrix}
\end{align}{}

Let us start from the complete Maxwell Equations for a non-dispersive medium, with $\mathbf{D}=\hat{\varepsilon} \mathbf{E}$ and $\mathbf{B}=\hat{\mu} \mathbf{H}$ (where we have already inverted the electromagnetic tensors to express $\mathbf{E}$ as $\hat{\varepsilon}^{-1}\mathbf{D}$ and $\mathbf{H}$ as $\hat{\mu}^{-1}\mathbf{H}$ ). In this instance, we shall restrict ourselves to the case of normal incidence, so that $\frac{\partial}{\partial x} = \frac{\partial}{\partial y} = 0$.

\begin{align}
    \PDT{D_y} &= \frac{1}{\mu_1}\bigg\{-\ALMUB g [-\STTP B_x + \CTTP B_y] + [1+\ALMUB \SSTP] \PDZ{B_x}- \ALMUB \CTP \STP \PDZ{B_y}  \bigg\} \label{eq:me1} \\
    %
    \PDT{D_x} &= \frac{1}{\mu_1}\bigg\{+\ALMUB g [+\CTTP B_x + \STTP B_y] - [1+\ALMUB \CSTP] \PDZ{B_y}+ \ALMUB \CTP \STP \PDZ{B_x}  \bigg\}\label{eq:me2}\\
    %
    \PDT{B_y} &= \frac{1}{\varepsilon_1}\bigg\{-\ALEB g [ -\STT D_x + \CTT D_y ]  - [1+\ALEB \CST] \PDZ{D_x} - \ALEB \CT \ST \PDZ{D_y} \bigg\}\label{eq:me3}\\
    %
    \PDT{B_x} &= \frac{1}{\varepsilon_1}\bigg\{+\ALEB g [+\CTT D_x + \STT D_y ] + [1+\ALEB \SST] \PDZ{D_y} + \ALEB \CT \ST \PDZ{D_x} \bigg\}\label{eq:me4}
\end{align}
where we defined $\bar{\alpha}_{\varepsilon}=-\frac{\alpha_\varepsilon}{1+2\alpha_\varepsilon}$ and $\bar{\alpha}_{\mu}=-\frac{\alpha_\mu}{1+2\alpha_\mu}$. We can rewrite this as:
\begin{align}
    \mu_1\PDT{D_y} &= +\PDZ{B_x} + \ALMUB g B_y- 4\bar{\alpha}_\mu g\CTP [-\STP B_x + \CTP B_y] - \ALMUB \STP [- \STP \PDZ{B_x} +\CTP \PDZ{B_y}]  \label{eq:me1_adj1} \\
    %
    \mu_1\PDT{D_x} &= -\PDZ{B_y} + \ALMUB g B_x + 4\bar{\alpha}_\mu g \STP[-\STP B_x + \CTP B_y ] - \ALMUB \CTP [-\STP \PDZ{B_x}+ \CTP \PDZ{B_y}] \label{eq:me2_adj1}\\
    %
    \varepsilon_1\PDT{B_y} &= - \PDZ{D_x} -\ALEB g D_y +4\bar{\alpha}_\varepsilon g \ST [ +\CT D_x + \ST D_y ]  - \ALEB \CT[+ \CT\PDZ{D_x} + \ST \PDZ{D_y}]  \label{eq:me3_adj1}\\
    %
    \varepsilon_1\PDT{B_x} &= +\PDZ{D_y} - \ALEB g D_x + 4\bar{\alpha}_\varepsilon g \CT [+\CT D_x + \ST D_y ] + \ALEB \ST [+\CT\PDZ{D_x} + \ST \PDZ{D_y}]  \label{eq:me4_adj1}
\end{align}
We now combine the above equations as follows:
\begin{align}
     \STP Eq.(\ref{eq:me1})&+\CTP Eq.(\ref{eq:me2})
     & 
     \CTP Eq.(\ref{eq:me1})&-\STP Eq.(\ref{eq:me2})\\
     \ST Eq.(\ref{eq:me3})&+\CT Eq.(\ref{eq:me4})
     &
     \CT Eq.(\ref{eq:me3})&-\ST Eq.(\ref{eq:me4})
\end{align}{}
and define the screwing coordinates:
\begin{align}
    x'(\theta) &= \CTH x + \STH y &
    y'(\theta) &= -\STH x + \CTH y 
\end{align}
% \begin{figure}
%     \centering
%     \includegraphics[width=\textwidth]{CoordsScrew.png}
%     \caption{Diagram depicting the screwing coordinates used.}
%     \label{fig:screwing_coordinates}
% \end{figure}
as well as the related screwing field components:
\begin{align}
        D_{x'}^+(\theta^+) &= \CTP D_x(\theta^+) + \STP D_y(\theta^+) & D_{y'}(\theta^+) 
        &= -\STP D_x(\theta^+) + \CTP D_y(\theta^+)\\
        B_{x'}^{-}(\theta) &= \CT B_x(\theta^-) + \ST B_y(\theta^-) & 
        B_{y'}^-(\theta^-) &= -\ST B_x(\theta^-) + \CT B_y(\theta^-)
\end{align}
where we need to account for the space-and time-dependence of $\theta^\pm(z,t)$ as we include the trigonometric functions within the derivatives via the chain-rule:
\begin{align}
    \CTH\PDZ{\psi}&=\PDZ{(\CTH \psi)} + g \STH \psi &     \STH\PDZ{\psi}&=\PDZ{(\STH \psi)} - g \CTH \psi\\
    \CTH\PDT{\psi}&=\PDT{(\CTH \psi)} - \OM \STH \psi&    \STH\PDT{\psi}&=\PDT{(\STH \psi)} + \OM \CTH \psi,
\end{align}{}
we thus arrive at the new set of equations:
\begin{align}
      \frac{\partial D_{x'}^{+}}{\partial t} + \OM D_{y'}^+  &= \frac{1}{\mu_1}\bigg\{-\frac{\partial B_{y'}^{+}}{\partial z} - g B_{x'}^+ - 2\bar{\alpha}_\mu \PDZ{B_{y'}^+} + [-\STP \PDX{} + \CTP \PDY{}]B_z \bigg\} \label{eq:coordinate_change1}\\
      %
      \frac{\partial D_{y'}^{+}}{\partial t} - \OM D_{x'}^+ &= \frac{1}{\mu_1}\bigg\{+\frac{\partial B_{x'}^{+}}{\partial z} - g B_{y'}^+ - 2\bar{\alpha}_\mu g  B_{y'}^+ - [+\CTP \PDX{} + \STP \PDY{}]B_z \bigg\}
    \label{eq:coordinate_change2}
    \\
    \frac{\partial B_{x'}^{-}}{\partial t} + \OM B_{y'}^-  &= \frac{1}{\mu_1}\bigg\{+\frac{\partial D_{y'}^{-}}{\partial z} + g D_{x'}^- + 2\bar{\alpha}_\varepsilon g  D_{x'}^- - [-\ST \PDX{}+\CT \PDY{}]D_z\bigg\}
    \label{eq:coordinate_change3}
    \\
    \frac{\partial B_{y'}^{-}}{\partial t} - \OM B_{x'}^-  &= \frac{1}{\mu_1}\bigg\{-\frac{\partial D_{x'}^{-}}{\partial z} +g D_{y'}^- - 2\bar{\alpha}_\varepsilon \frac{\partial D_{x'}^-}{\partial z} + [+\CT \PDX{} + \ST \PDY{}]D_z\bigg\}
    \label{eq:coordinate_change4}
    %\frac{1}{\varepsilon_1}\bigg[ \frac{\partial H_{x'}^-}{} \bigg]
\end{align}
We now want to move to a basis of forward and backward-propagating waves. After taking the combinations:
\begin{align}
    Eq. \ref{eq:coordinate_change1} &+ \sqrt{\frac{\varepsilon_1}{\mu_1}}Eq. \ref{eq:coordinate_change4} & Eq. \ref{eq:coordinate_change1}&- \sqrt{\frac{\varepsilon_1}{\mu_1}}Eq. \ref{eq:coordinate_change4} \\ 
    Eq. \ref{eq:coordinate_change2} &+ \sqrt{\frac{\varepsilon_1}{\mu_1}}Eq. \ref{eq:coordinate_change3} & Eq. \ref{eq:coordinate_change2}&- \sqrt{\frac{\varepsilon_1}{\mu_1}}Eq. \ref{eq:coordinate_change3},
\end{align}
rescaling the $B$ fields as $\bar{B} = \sqrt{\varepsilon_1/\mu_1}B$, and defining the sums and differences between the electric and the magnetic modulations $\bar{\alpha}^\pm = \bar{\alpha}_\varepsilon \pm \bar{\alpha}_\mu$ we can use the identities: 
\begin{align}
    \bar{\alpha}_\varepsilon D_i + \bar{\alpha}_\mu \bar{B}_j &= \frac{1}{2}[(\bar{\alpha}_\varepsilon+\bar{\alpha}_\mu)(D_i+\bar{B}_j)+(\bar{\alpha}_\varepsilon-\bar{\alpha}_\mu)(D_i-\bar{B}_j)]=\frac{1}{2}[\bar{\alpha}^+(D_i+\bar{B}_j)+\bar{\alpha}^-(D_i-\bar{B}_j)] \\
    \bar{\alpha}_\varepsilon D_i - \bar{\alpha}_\mu \bar{B}_j &= \frac{1}{2}[(\bar{\alpha}_\varepsilon-\bar{\alpha}_\mu)(D_i+\bar{B}_j)+(\bar{\alpha}_\varepsilon+\bar{\alpha}_\mu)(D_i-\bar{B}_j)]=\frac{1}{2}[\bar{\alpha}^-(D_i+\bar{B}_j)+\bar{\alpha}^+(D_i-\bar{B}_j)]
\end{align}
to obtain:
\begin{align}
     \PDT{}(D_{x'}^{+}+\bar{B}_{y'}^{-}) + \OM (D_{y'}^+ - \bar{B}_{x'}^{-})  &= c_1\bigg\{-\PDZ{}(D_{x'}^{-}+\bar{B}_{y'}^{+}) +g (D_{y'}^{-} - \bar{B}_{x'}^+) - \PDZ{} \bigg[ \bar{\alpha}^+(D_{x'}^-+\bar{B}_{y'}^+)+\bar{\alpha}^-(D_{x'}^--\bar{B}_{y'}^+) \bigg] 
     \bigg\} \label{eq:for_back1}\\
    %
    \PDT{}(D_{x'}^{+}-\bar{B}_{y'}^{-}) + \OM (D_{y'}^+ + \bar{B}_{x'}^{-})  &= c_1\bigg\{+\PDZ{} (D_{x'}^{-}-\bar{B}_{y'}^{+}) -g (D_{y'}^{-} + \bar{B}_{x'}^+) + \PDZ{} \bigg[\bar{\alpha}^-(D_{x'}^-+\bar{B}_{y'}^+)+\bar{\alpha}^+(D_{x'}^--\bar{B}_{y'}^+)\bigg]\bigg\} \label{eq:for_back2}\\
    %
    \PDT{}(D_{y'}^{+}+\bar{B}_{x'}^{-}) - \OM (D_{x'}^+ - \bar{B}_{y'}^{-})  &= c_1\bigg\{+\PDZ{} (D_{y'}^{-}+\bar{B}_{x'}^{+}) + g (D_{x'}^{-} - \bar{B}_{y'}^+) + g [\bar{\alpha}^-(D_{x'}^-+\bar{B}_{y'}^+)+\bar{\alpha}^+(D_{x'}^--\bar{B}_{y'}^+)] \bigg\} \label{eq:for_back3}\\
    %
    \PDT{}(D_{y'}^{+}-\bar{B}_{x'}^{-}) - \OM (D_{x'}^+ + \bar{B}_{y'}^{-})  &= c_1\bigg\{ -\PDZ{}(D_{y'}^{-}-\bar{B}_{x'}^{+}) - g (D_{x'}^{-} + \bar{B}_{y'}^+) - g[\bar{\alpha}^+(D_{x'}^-+\bar{B}_{y'}^+)+\bar{\alpha}^- (D_{x'}^--\bar{B}_{y'}^+)]\bigg\} \label{eq:for_back4}
\end{align}
Which allowed us to successfully remove all the $\theta$-dependent terms from the equations. We are now left with 8 combinations of fields (due to the presence of both $\pm\phi$ terms), but only 4 equations. In order to reduce the number of variables, we can recognise that the individual terms above can be expanded as:
\begin{align}
    D_{x'}^{+}+\bar{B}_{y'}^{-} &= [\CTP (+D_x) + \CT (+\bar{B}_y)] + [\STP (+D_y) + \ST (-\bar{B}_x)]\\ 
    D_{x'}^{+}-\bar{B}_{y'}^{-} &= [\CTP (+D_x) + \CT (-\bar{B}_y)] + [\STP (+D_y) + \ST (+\bar{B}_x)] \\
    D_{y'}^{+}+\bar{B}_{x'}^{-} &= [\STP (-D_x) + \ST (+\bar{B}_y)] + [\CTP (+D_y) + \CT (+\bar{B}_x)] \\
    D_{y'}^{+}-\bar{B}_{x'}^{-} &= [\STP (-D_x) + \ST (-\bar{B}_y)] + [\CTP (+D_y) + \CT (-\bar{B}_x)] \\
    D_{x'}^{-}+\bar{B}_{y'}^{+} &= [\CT (+D_x) + \CTP (+\bar{B}_y)] + [\ST (+D_y) + \STP(-\bar{B}_x)] \\
    D_{x'}^{-}-\bar{B}_{y'}^{+} &= [\CT (+D_x) + \CTP (-\bar{B}_y)] + [\ST (+D_y) + \STP(+\bar{B}_x)] \\
    D_{y'}^{-}+\bar{B}_{x'}^{+} &= [\ST (-D_x) + \STP (+\bar{B}_y)] + [\CT (+D_y) + \CTP(-\bar{B}_x)] \\
    D_{y'}^{-}-\bar{B}_{x'}^{+} &= [\ST (-D_x) + \STP (-\bar{B}_y)] + [\CT (+D_y) + \CTP(+\bar{B}_x)] 
\end{align}
and use the relations:
\begin{align}
    \CTP \Psi +\CT \Phi &= \frac{1}{2} [(\CTP+\CT)(\Psi+\Phi) +(\CTP -\CT)(\Psi-\Phi)] 
    = \CTH \CP (\Psi + \Phi) - \STH\SP (\Psi-\Phi)\\
    \STP \Psi +\ST \Phi &= \frac{1}{2} [(\STP+\ST)(\Psi+\Phi) +(\STP -\ST)(\Psi-\Phi)] 
    = \STH \CP (\Psi + \Phi) + \CTH\SP (\Psi-\Phi),
\end{align}{}
which follow from the the trigonometric identities:
\begin{align}
    \CTP + \CT &= 2\CTH \CP & \CTP - \CT &= -2\STH \SP \\
    \STP + \ST &= 2\STH \CP & \STP - \ST &= +2\CTH \SP
\end{align}{}
to write the equations above only as a function of the forward and backward-propagating fields:
\begin{align}
    F_{x'}^\rightharpoonup &= \CTH (D_x+\bar{B}_y)+\STH [D_y + (-\bar{B}_x)] & F_{y'}^\rightharpoonup &= -\STH (D_x+\bar{B}_y)+\CTH [D_y + (-\bar{B}_x)] \\
    %
    F_{x'}^\leftharpoonup &= \CTH (D_x-\bar{B}_y)+\STH [D_y - (-\bar{B}_x)] & F_{y'}^\leftharpoonup &= -\STH (D_x-\bar{B}_y)+\CTH [D_y - (-\bar{B}_x)]
\end{align}
so that:
\begin{align}
    D_{x'}^+ + \bar{B}_{y'}^- &= \cos{(\phi)} F_{x'}^\rightharpoonup + \sin{(\phi)} F_{y'}^\leftharpoonup &
    %
    D_{x'}^- + \bar{B}_{y'}^+ &= -\sin{(\phi)} F_{y'}^\leftharpoonup + \cos{(\phi)} F_{x'}^\rightharpoonup     \\
    %
    D_{x'}^+ - \bar{B}_{y'}^- &= \cos{(\phi)} F_{x'}^\leftharpoonup + \sin{(\phi)} F_{y'}^\rightharpoonup &
    %
    D_{x'}^- - \bar{B}_{y'}^+ &= -\sin{(\phi)} F_{y'}^\rightharpoonup + \cos{(\phi)} F_{x'}^\leftharpoonup \\
    %
    D_{y'}^+ + \bar{B}_{x'}^- &= -\sin{(\phi)} F_{x'}^\rightharpoonup + \cos{(\phi)} F_{y'}^\leftharpoonup &
    %
    D_{y'}^- + \bar{B}_{x'}^+ &= \cos{(\phi)} F_{y'}^\leftharpoonup + \sin{(\phi)} F_{x'}^\rightharpoonup \\
    %
    D_{y'}^+ - \bar{B}_{x'}^- &= -\sin{(\phi)} F_{x'}^\leftharpoonup + \cos{(\phi)} F_{y'}^\rightharpoonup &
    %
    D_{y'}^- - \bar{B}_{x'}^+ &= \cos{(\phi)} F_{y'}^\rightharpoonup + \sin{(\phi)} F_{x'}^\leftharpoonup 
\end{align}
Substituting the latter into Eqs. \ref{eq:for_back1}-\ref{eq:for_back4}, we derive four coupled equations, which no longer depend on the spatiotemporal variable $\theta$, which has been absorbed into our new set of basis functions. These can be cast into a standard eigenvalue problem by taking:
\begin{align}
     \PDT{}(\cos{(\phi)} F_{x'}^\rightharpoonup + \sin{(\phi)} F_{y'}^\leftharpoonup) &= - \OM (-\sin{(\phi)} F_{x'}^\leftharpoonup + \cos{(\phi)} F_{y'}^\rightharpoonup) +c_1\bigg\{  -\PDZ{}(-\sin{(\phi)} F_{y'}^\leftharpoonup + \cos{(\phi)} F_{x'}^\rightharpoonup) +g (\cos{(\phi)} F_{y'}^\rightharpoonup + \sin{(\phi)} F_{x'}^\leftharpoonup) \nonumber\\
     &- \PDZ{} \bigg[ \bar{\alpha}^+(-\sin{(\phi)} F_{y'}^\leftharpoonup + \cos{(\phi)} F_{x'}^\rightharpoonup )+\bar{\alpha}^-(-\sin{(\phi)} F_{y'}^\rightharpoonup + \cos{(\phi)} F_{x'}^\leftharpoonup) \bigg] \bigg\} \label{eq:for_backSub1}\\
    %%%%%%
    \PDT{}(\cos{(\phi)} F_{x'}^\leftharpoonup + \sin{(\phi)} F_{y'}^\rightharpoonup) &= - \OM (-\sin{(\phi)} F_{x'}^\rightharpoonup + \cos{(\phi)} F_{y'}^\leftharpoonup ) + c_1\bigg\{+\PDZ{} (-\sin{(\phi)} F_{y'}^\rightharpoonup + \cos{(\phi)} F_{x'}^\leftharpoonup) -g (\cos{(\phi)} F_{y'}^\leftharpoonup + \sin{(\phi)} F_{x'}^\rightharpoonup)\nonumber \\& + \PDZ{} \bigg[\bar{\alpha}^-(-\sin{(\phi)} F_{y'}^\leftharpoonup + \cos{(\phi)} F_{x'}^\rightharpoonup )+\bar{\alpha}^+(-\sin{(\phi)} F_{y'}^\rightharpoonup + \cos{(\phi)} F_{x'}^\leftharpoonup )\bigg]\bigg\} \label{eq:for_backSub2}\\
    %%%%%%
    \PDT{}(-\sin{(\phi)} F_{x'}^\rightharpoonup + \cos{(\phi)} F_{y'}^\leftharpoonup) &= + \OM (\cos{(\phi)} F_{x'}^\leftharpoonup + \sin{(\phi)} F_{y'}^\rightharpoonup)  +c_1\bigg\{+\PDZ{} (\cos{(\phi)} F_{y'}^\leftharpoonup + \sin{(\phi)} F_{x'}^\rightharpoonup) + g (-\sin{(\phi)} F_{y'}^\rightharpoonup + \cos{(\phi)} F_{x'}^\leftharpoonup ) \nonumber \\ & + g [\bar{\alpha}^-(-\sin{(\phi)} F_{y'}^\leftharpoonup + \cos{(\phi)} F_{x'}^\rightharpoonup )+\bar{\alpha}^+(-\sin{(\phi)} F_{y'}^\rightharpoonup + \cos{(\phi)} F_{x'}^\leftharpoonup )] \bigg\} \label{eq:for_backSub3}\\
    %%%%%%
    \PDT{}(-\sin{(\phi)} F_{x'}^\leftharpoonup + \cos{(\phi)} F_{y'}^\rightharpoonup) &= + \OM (\cos{(\phi)} F_{x'}^\rightharpoonup + \sin{(\phi)} F_{y'}^\leftharpoonup) + c_1\bigg\{ -\PDZ{}(\cos{(\phi)} F_{y'}^\rightharpoonup + \sin{(\phi)} F_{x'}^\leftharpoonup ) - g (-\sin{(\phi)} F_{y'}^\leftharpoonup + \cos{(\phi)} F_{x'}^\rightharpoonup ) \nonumber \\& - g[\bar{\alpha}^+(-\sin{(\phi)} F_{y'}^\leftharpoonup + \cos{(\phi)} F_{x'}^\rightharpoonup )+\bar{\alpha}^- (-\sin{(\phi)} F_{y'}^\rightharpoonup + \cos{(\phi)} F_{x'}^\leftharpoonup )]  \bigg\} \label{eq:for_backSub4}
\end{align}
which we conveniently simplify by taking:
\begin{align}
    \cos{(\phi)}Eq.\ref{eq:for_backSub1}&-\sin{(\phi)}Eq.\ref{eq:for_backSub3} &  \sin{(\phi)}Eq.\ref{eq:for_backSub2}&+\cos{(\phi)}Eq.\ref{eq:for_backSub4} &
    %
    \cos{(\phi)}Eq.\ref{eq:for_backSub2}&-\sin{(\phi)}Eq.\ref{eq:for_backSub4} & \sin{(\phi)}Eq.\ref{eq:for_backSub1}&+\cos{(\phi)}Eq.\ref{eq:for_backSub3} 
\end{align}
leaving us with:
\begin{align}
\PDT{\FXP}=  + (c_1 g-\Omega) \FYP +c_1\bigg\{-\PDZ{\FXP} & -\cos{(\phi)}\PDZ{}[\bar{\alpha}^+(-\SP \FYM +\CP \FXP)+\bar{\alpha}^-(-\SP\FYP+\CP\FXM)] \nonumber \\
&- \SP g [\ALM (-\SP \FYM + \CP \FXP) +\ALP(-\SP \FYP+\CP \FXM)]  \bigg\}\\
%
\PDT{\FYP}=-(c_1 g-\Omega)\FXP +c_1\bigg\{-\PDZ{\FYP}&+ \SP \PDZ{}[\ALM(-\SP \FYM+\CP \FXP)+\ALP(-\SP\FYP+\CP \FXM)]\nonumber \\
& -\CP g [\ALP (-\SP \FYM + \CP \FXP)+\ALM(-\SP \FYP+\CP\FXM)] \bigg\}\\
%
\PDT{\FXM}=-(c_1 g+\Omega)\FYM +c_1\bigg\{+\PDZ{\FXM}&+\CP \PDZ{}[\ALM(-\SP \FYM +\CP \FXP)+\ALP(-\SP \FYP+\CP \FXM)] \nonumber\\
&+\SP g[\ALP(-\SP\FYM+\CP\FXP)+\ALM(-\SP \FYP+\CP\FXM)]\bigg\} \\
%
\PDT{\FYM}=  + (c_1 g+\Omega) \FXM +c_1\bigg\{+\PDZ{\FYM}& -\SP\PDZ{}[\bar{\alpha}^+(-\SP \FYM +\CP \FXP)+\bar{\alpha}^-(-\SP\FYP+\CP\FXM)] \nonumber\\&+ \CP g [\ALM (-\SP \FYM + \CP \FXP) +\ALP(-\SP \FYP+\CP \FXM)]  \bigg\}
\end{align}

Inserting Bloch wave solutions: $\Psi = e^{i(kz-\omega t)}\sum_n \psi_n e^{(2n-1)i(gz-\Omega t)}$ yields (defining $\omega_n = \omega+(2n-1)\Omega$ and $k_n = k+(2n-1)g$):
\begin{align}
\omega_n\FXP_n=  +i (c_1 g-\Omega) \FYP_n -c_1\bigg\{- k_n \FXP_n & - k_n\cos{(\phi)}[\bar{\alpha}^+(-\SP \FYM_n +\CP \FXP_n)+\bar{\alpha}^-(-\SP\FYP_n+\CP\FXM_n)] \nonumber \\
&+i g\SP [\ALMB (-\SP \FYM_n + \CP \FXP_n) +\ALPB(-\SP \FYP_n+\CP \FXM_n)]  \bigg\} \label{eq:final_ME_1}\\
%
\omega_n\FYP_n=-i(c_1 g-\Omega)\FXP_n -c_1\bigg\{-k_n\FYP_n&+ k_n\SP [\ALMB(-\SP \FYM_n+\CP \FXP_n)+\ALPB(-\SP\FYP_n+\CP \FXM_n)]\nonumber \\
& +ig\CP [\ALPB (-\SP \FYM_n + \CP \FXP_n)+\ALMB(-\SP \FYP_n+\CP\FXM_n)]  \bigg\}\label{eq:final_ME_2}\\
%
\omega_n\FXM_n=-i(c_1 g+\Omega)\FYM_n -c_1\bigg\{+k_n \FXM_n&+k_n\CP [\ALMB(-\SP \FYM_n +\CP \FXP_n)+\ALPB(-\SP \FYP_n+\CP \FXM_n)] \nonumber\\
&-ig\SP [\ALPB(-\SP\FYM_n+\CP\FXP_n)+\ALMB(-\SP \FYP_n+\CP\FXM_n)] \label{eq:final_ME_3}\\
%
\omega_n\FYM_n=  +i (c_1 g+\Omega) \FXM_n -c_1\bigg\{+k_n \FYM_n& -k_n\SP [\bar{\alpha}^+(-\SP \FYM_n +\CP \FXP_n)+\bar{\alpha}^-(-\SP\FYP_n+\CP\FXM_n)] \nonumber\\&-i g\CP  [\ALMB (-\SP \FYM_n + \CP \FXP_n) +\ALPB(-\SP \FYP_n+\CP \FXM_n)]\bigg\}\label{eq:final_ME_4} 
\end{align}

Abbreviating the notation further with $s=\sin(\phi)$, $c=\cos(\phi)$, $\BG = c_1 g$ and $\BKN = c_1 k_n$ we can rewrite the latter as a standard eigenvalue problem:

\begin{align}
    \omega_n \begin{pmatrix} \mathbf{F}^{\rightharpoonup}_n \\ \mathbf{F}^{\leftharpoonup}_n \end{pmatrix} = 
    \begin{pmatrix} 
    \mathbb{M}^{\rightharpoonup}_{\rightharpoonup}_n& \mathbb{M}^{\rightharpoonup}_{\leftharpoonup}_n \\ \mathbb{M}^{\leftharpoonup}_{\rightharpoonup}_n& \mathbb{M}^{\leftharpoonup}_{\leftharpoonup}_n
    \end{pmatrix}
    \begin{pmatrix} \mathbf{F}^{\rightharpoonup}_n \\ \mathbf{F}^{\leftharpoonup}_n \end{pmatrix} \label{eq:compactEigenproblem}
\end{align}
where the four matrices couple the forward and backward propagating states $\mathbf{F}^{\rightharpoonup}_n = (\FXP_n;\FYP_n)$ and $\mathbf{F}^{\leftharpoonup}_n = (\FXM_n;\FYM_n)$, and read:
\begin{align}
    \mathbb{M}^{\rightharpoonup}_{\rightharpoonup}_n &= \begin{pmatrix} \BKN+ \CP[\ALPB \CP \BKN-i\ALMB \SP \BG] & +i(\BG-\OM)+\SP[-\ALMB\CP\BKN+i\ALPB\SP\BG] \\
    -i(\BG-\OM)+\CP[-\ALMB\SP\BKN-i\ALPB\CP\BG] & \BKN + \SP [\ALPB\SP \BKN + i \ALMB \CP \BG]
    \end{pmatrix}  \\
    %
    \mathbb{M}^{\rightharpoonup}_{\leftharpoonup}_n &= \begin{pmatrix} \CP[\ALMB \CP \BKN -i\ALPB \SP \BG] & \SP[-\ALPB\CP\BKN+i\ALMB\SP\BG] \\
    \CP[-\ALPB \SP \BKN -i \ALMB \CP \BG] & \SP[\ALMB \SP \BKN +i \ALPB \CP \BG]
    \end{pmatrix}  \\
    %
    \mathbb{M}^{\leftharpoonup}_{\rightharpoonup}_n &= \begin{pmatrix} \CP[-\ALMB \CP \BKN +i\ALPB \SP \BG] & \SP[\ALPB\CP \BKN-i\ALMB\SP\BG] \\
    \CP[\ALPB \SP \BKN +i \ALMB \CP \BG] & \SP[-\ALMB \SP \BKN -i \ALPB \CP \BG]
    \end{pmatrix} \\
    %
    \mathbb{M}^{\leftharpoonup}_{\leftharpoonup}_n &= \begin{pmatrix} -\BKN + \CP[-\ALPB \CP \BKN +i\ALMB \SP \BG] & -i(\OM+\BG) + \SP[\ALMB\CP \BKN-i\ALPB\SP\BG] \\
    +i(\OM+\BG) +\CP[\ALMB \SP \BKN + i \ALPB \CP \BG] & -\BKN + \SP[-\ALPB \SP \BKN -i \ALMB \CP \BG]
    \end{pmatrix} 
\end{align}

Note that one peculiarity of this solution strategy is that the ordering of the bands is modified as a result of the coupling between harmonics inherent to the screwing basis vectors. As a result, diagonalizing Eq. \ref{eq:compactEigenproblem} yields all the eigenvalues, but not in increasing order.

\section{Change of basis in matrix form}

Having derived the correct basis which uncouples the bands, we can rewrite the change of basis in a more compact matrix form $\mathbb{F} = \hat{S}\mathbb{G}$ as:
\begin{align}
    \begin{pmatrix}
    F_{x'}^\rightharpoonup \\ F_{y'}^\rightharpoonup \\ F_{x'}^\leftharpoonup \\ F_{y'}^\leftharpoonup
    \end{pmatrix} &= \begin{pmatrix}
    \CTH & \CTH & \STH & \STH \\
    -\STH & -\STH & \CTH & \CTH \\
    \CTH & -\CTH & \STH & -\STH \\
    -\STH & \STH & \CTH & - \CTH
    \end{pmatrix}
    \begin{pmatrix}
    D_x \\ \bar{B}_y \\ D_y \\ -\bar{B}_x
    \end{pmatrix}
\end{align}
Since the matrix above is orthogonal, its inverse corresponds to its transpose, and we can write immediately:
\begin{align}
    \begin{pmatrix}
    D_x \\ \bar{B}_y \\ D_y \\ -\bar{B}_x
    \end{pmatrix}
    &= \begin{pmatrix}
    \CTH & -\STH & \CTH & -\STH \\
    \CTH & -\STH & -\CTH & \STH \\
    \STH & \CTH & \STH & \CTH \\
    \STH & \CTH & -\STH & - \CTH
    \end{pmatrix}
    \begin{pmatrix}
    F_{x'}^\rightharpoonup \\ F_{y'}^\rightharpoonup \\ F_{x'}^\leftharpoonup \\ F_{y'}^\leftharpoonup
    \end{pmatrix}
\end{align}
Thus, we have:
\begin{align}
    D_{x} &= \CTH F_{x'}^\rightharpoonup -\STH F_{y'}^\rightharpoonup+
    \CTH F_{x'}^\leftharpoonup -\STH F_{y'}^\leftharpoonup &
    %
    D_{y} &= \STH F_{x'}^\rightharpoonup +\CTH F_{y'}^\rightharpoonup+
    \STH F_{x'}^\leftharpoonup +\CTH F_{y'}^\leftharpoonup \\
    %
    \bar{B}_{y} &= \CTH F_{x'}^\rightharpoonup -\STH F_{y'}^\rightharpoonup - \CTH F_{x'}^\leftharpoonup +\STH F_{y'}^\leftharpoonup &
    %
    -\bar{B}_x &= \STH F_{x'}^\rightharpoonup + \CTH F_{y'}^\rightharpoonup - \STH F_{x'}^\leftharpoonup -\CTH F_{y'}^\leftharpoonup
\end{align}






\section{Numerical Floquet-Bloch Calculations}

Here we solve for $k(\omega)$, in terms of the $E$ and $H$ fields. Note that, due to a dual symmetry between space and time we can to solve for $\omega(k)$ and for  with the same sets of equations following the substitutions:
\begin{align}
    D&\leftrightarrow E &
    B&\leftrightarrow H &
    \alpha_{\varepsilon/\mu} &\leftrightarrow \bar{\alpha}_{\varepsilon/\mu} =-\frac{\alpha_{\varepsilon/\mu}}{1+2\alpha_{\varepsilon/\mu}} &
    k &\leftrightarrow \omega & g&\leftrightarrow \Omega 
\end{align}
Since we are going to use $E,H$ as our variables, it is useful to write down the time derivatives of the material tensors:
\begin{align}
    \PDT{\varepsilon} &= 2\EB\ALED\OM \begin{pmatrix}
    \STT&-\CTT&0\\
    -\CTT&-\STT&0\\
    0&0&0
    \end{pmatrix} & \PDT{\mu} &=  
    2\MB\ALMUD\OM \begin{pmatrix}
    -\STTP&\CTTP&0\\
    \CTTP&\STTP&0\\
    0&0&0
    \end{pmatrix}
\end{align}{}
so that Maxwell's Equations give us, for normal incidence:
\begin{align}
    \frac{1}{\MB}\bigg(\PDY{E_z} - \PDZ{E_y}\bigg) &= + \ALMU \OM (\STTP H_x - \CTTP H_y) - (1+\ALMU \SSTP) \PDT{H_x} + \ALMU \STP\CTP \PDT{H_y} \\
    %
    -\frac{1}{\MB}\bigg(\PDX{E_z} - \PDZ{E_x}\bigg) &= + \ALMU \OM (-\CTTP H_x - \STTP H_y)  + \ALMU \STP\CTP \PDT{H_x} - (1+\ALMU \CSTP)\PDT{H_y}\\
    %
    \vspace{2pt}
    %
    \frac{1}{\EB}\bigg(\PDY{H_z} - \PDZ{H_y}\bigg) &=  \ALE \OM (\STT E_x - \CTT E_y) + (1+\ALE \CST) \PDT{E_x} + \ALE \ST\CT \PDT{E_y} \\
    %
    -\frac{1}{\EB}\bigg(\PDX{H_z} - \PDZ{H_x}\bigg) &=  \ALE \OM (-\CTT E_x - \STT E_y)  + \ALE \ST\CT \PDT{E_x} + (1+\ALE \SST)\PDT{E_y}
\end{align}{}
and using double-angle formulae:
\begin{align}
    \frac{1}{\MB}\bigg(\PDY{E_z} - \PDZ{E_y}\bigg) &= + \ALMU \OM (\STTP H_x - \CTTP H_y) - [1+\ALMUD (1-\CTTP)] \PDT{H_x} + \ALMUD \STTP \PDT{H_y} \\
    %
    -\frac{1}{\MB}\bigg(\PDX{E_z} - \PDZ{E_x}\bigg) &= + \ALMU \OM (-\CTTP H_x - \STTP H_y)  + \ALMUD \STTP \PDT{H_x} - [1+\ALMUD (1+\CTTP)]\PDT{H_y}\\
    %
    \vspace{2pt}
    %
    \frac{1}{\EB}\bigg(\PDY{H_z} - \PDZ{H_y}\bigg) &= + \ALE \OM (\STT E_x - \CTT E_y) + [1+\ALED (1+\CTT)] \PDT{E_x} + \ALED \STT \PDT{E_y} \\
    %
    -\frac{1}{\EB}\bigg(\PDX{H_z} - \PDZ{H_x}\bigg) &= + \ALE \OM (-\CTT E_x - \STT E_y)  + \ALED \STT \PDT{E_x} + [1+\ALED(1-\CTT)]\PDT{E_y}
\end{align}{}

We now assume ansatz of the form: $\psi \sim e^{i(k z-\omega t)}\sum_n \psi_n e^{2in(gz-\Omega t)}$, so that the respective derivatives yield:
\begin{align}
    \frac{\partial \psi}{\partial t} &= -i e^{i(k z -\omega t)} \sum_n (\omega+2n\Omega) \psi_n e^{2in(gz-\Omega t)} \\
    %
    \frac{\partial \psi}{\partial z} &= i e^{i(k z -\omega t)} \sum_n (k+2ng) \psi_n e^{2in(gz-\Omega t)}
\end{align}
We can thus write Maxwell's equations as:
\begin{align}
     \PDZ{E_y} &= \MB \bigg[ \ALMU \OM (-\STTP H_x + \CTTP H_y)- \ALMUD \STTP \PDT{H_y} + [1+\ALMUD (1-\CTTP)] \PDT{H_x} \bigg] \\
    %
    \PDZ{E_x} &= \MB  \bigg[\ALMU \OM (-\CTTP H_x - \STTP H_y)  + \ALMUD \STTP \PDT{H_x} - [1+\ALMUD (1+\CTTP)]\PDT{H_y}\bigg]\\
    %
    \vspace{2pt}
    %
    \PDZ{H_y}&= \EB\bigg[\ALE \OM (-\STT E_x + \CTT E_y)-\ALED \STT \PDT{E_y} - [1+\ALED( 1+\CTT)] \PDT{E_x} \bigg] \\
    %
    \PDZ{H_x} &= \EB\bigg[\ALE \OM (-\CTT E_x - \STT E_y)  + \ALED \STT \PDT{E_x} + [1+\ALED (1-\CTT)]\PDT{E_y}\bigg]
\end{align}{}
so that substituting in the fields yields the eigenvalue problem for $k$:
\begin{align}
(k+2n'g)E_{y,n'} &= \sum_n \bigg\{
 \ALMUD \OM [(\delta_{n',n+1}-\delta_{n',n-1})\HXN  -i(\delta_{n',n+1}+\delta_{n',n-1})\HYN ] \\ &
-\frac{i}{2}\ALMUD (\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1}) \HYN \\ &
- (1+\ALMUD) (\omega+2n\OM)\delta_{n',n} \HXN \nonumber  + \frac{1}{2}\ALMUD (\omega+2n\OM) (\delta_{n',n+1}+\delta_{n',n-1}) \HXN  \bigg\} \nonumber \\
%
(k+2n'g)E_{x,n'} &= \sum_n \bigg\{ \ALMUD \OM [+(\delta_{n',n+2}-\delta_{n',n-1})\HYN +i(\delta_{n',n+1}+\delta_{n',n-1})\HXN ]\\ & 
+\frac{i}{2}\ALMUD (\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1}) \HXN  \\ &
+  (1+\ALMUD) (\omega+2n\OM)\delta_{n',n} \HYN \nonumber   + \frac{1}{2}\ALMUD (\omega+2n\OM) (\delta_{n',n+1}+\delta_{n',n-1}) \HYN \bigg\} \nonumber \\
%
(k+2n'g)H_{y,n'} &= \sum_n \bigg\{  \ALED \OM [+(\delta_{n',n+1}-\delta_{n',n-1})\EXN  -i(\delta_{n',n+1}+\delta_{n',n-1})\EYN ] \\ &-\frac{i}{2}\ALED (\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1}) \EYN  \\ & 
+ (1+\ALED) (\omega+2n\OM)\delta_{n',n} \EXN \nonumber  + \frac{1}{2}\ALED (\omega+2n\OM) (\delta_{n',n+1}+\delta_{n',n-1}) \EXN \bigg\} \nonumber \\
%
(k+2n'g)H_{x,n'} &= 
\sum_n \bigg\{ \ALED \OM [+(\delta_{n',n+1}-\delta_{n',n-1})\EYN +i(\delta_{n',n+1}+\delta_{n',n-1})\EXN ] \\&
+\frac{i}{2}\ALED (\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1}) \EXN \\&
-  (1+\ALMUD) (\omega+2n\OM)\delta_{n',n} \EYN \nonumber   + \frac{1}{2}\ALED (\omega+2n\OM) (\delta_{n',n+1}+\delta_{n',n-1}) \EYN  \bigg\} \nonumber 
\end{align}

In matrix form, we can write this system as:

\begin{align}
k_n
    \begin{pmatrix}
        \mathbf{E}_y \\ \mathbf{E}_x \\ \mathbf{H}_y \\ \mathbf{H}_x  
    \end{pmatrix}
    &=
    \begin{pmatrix}
        \mathbb{M}^{E_y}_{E_y} & 0 & \mathbb{M}^{E_y}_{H_y} & \mathbb{M}^{E_y}_{H_x} \\
        0 & \mathbb{M}^{E_x}_{E_x} & \mathbb{M}^{E_x}_{H_y} & \mathbb{M}^{E_x}_{H_x} \\
        \mathbb{M}^{H_y}_{E_y} & \mathbb{M}^{H_y}_{E_x} & \mathbb{M}^{H_y}_{H_y} & 0 \\
        \mathbb{M}^{H_x}_{E_y} & \mathbb{M}^{H_x}_{E_x} & 0 & \mathbb{M}^{E_x}_{E_x} \\
    \end{pmatrix}
    \begin{pmatrix}
        \mathbf{E}_y \\ \mathbf{E}_x \\ \mathbf{H}_y \\ \mathbf{H}_x  
    \end{pmatrix}
\end{align}

where the respective $\mathbb{M}$ matrices are tridiagonal, and read:

\begin{align}
    (M^{E_y}_{E_y})_{n',n} &= -2n g \delta_{n,n'} \\
    %
    (M^{E_x}_{E_x})_{n',n} &= -2n g \delta_{n,n'} \\
    %
    (M^{H_y}_{H_y})_{n',n} &= -2n g \delta_{n,n'} \\
    %
    (M^{H_x}_{H_x})_{n',n} &= -2n g \delta_{n,n'} \\
    %
    (M^{E_y}_{H_y})_{n',n} &= -\frac{k_y k_x}{\omega+2n\Omega} \delta_{n',n}+ i\ALMUD\OM (\delta_{n',n+1}+\delta_{n',n-1}) +i \frac{\ALMUD}{2}(\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1}) \\
    %
    (M^{E_y}_{H_x})_{n',n} &= \frac{k_y^2}{\omega+2n\Omega}\delta_{n',n}- \ALMUD\OM (\delta_{n',n+1}-\delta_{n',n-1}) -(1+\ALMUD)(\omega+2n\OM)\delta_{n,n'}- \frac{1}{2}\ALMUD(\omega+2n\OM)(\delta_{n',n+1}+\delta_{n',n-1})\\ 
    %
    (M^{E_x}_{H_y})_{n',n} &= -\frac{k_x^2}{\omega+2n\Omega}\delta_{n',n}- \ALMUD\OM (\delta_{n',n+1}-\delta_{n',n-1}) +(1-\ALMUD)(\omega+2n\OM)\delta_{n,n'}- \frac{1}{2}\ALMUD(\omega+2n\OM)(\delta_{n',n+1}+\delta_{n',n-1})\\
    %
    (M^{E_x}_{H_x})_{n',n} &= \frac{k_x k_y}{\omega+2n\Omega}\delta_{n',n}-i\ALMUD\OM (\delta_{n',n+1}+\delta_{n',n-1}) -i \frac{\ALMUD}{2}(\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1}) \\
    %
    (M^{H_y}_{E_y})_{n',n} &= \frac{k_x k_y}{\omega+2n\Omega}\delta_{n',n}-i \ALED\OM (\delta_{n',n+1}+\delta_{n',n-1}) +i \frac{\ALED}{2}(\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1}) \\
    %
    (M^{H_y}_{E_x})_{n',n} &= -\frac{k_y^2}{\omega+2n\Omega}\delta_{n',n} + \ALED\OM (\delta_{n',n+1}-\delta_{n',n-1}) +(1-\ALED)(\omega+2n\OM)\delta_{n,n'} + \frac{1}{2}\ALED(\omega+2n\OM)(\delta_{n',n+1}+\delta_{n',n-1})\\ 
    %
    (M^{H_x}_{E_y})_{n',n} &= \frac{k_x^2}{\omega+2n\Omega}\delta_{n',n} + \ALED\OM (\delta_{n',n+1}-\delta_{n',n-1}) -(1+\ALMUD)(\omega+2n\OM)\delta_{n,n'}+ \frac{1}{2}\ALED(\omega+2n\OM)(\delta_{n',n+1}+\delta_{n',n-1})\\
    %
    (M^{H_x}_{E_x})_{n',n} &= - \frac{k_x k_y}{\omega+2n\Omega}\delta_{n',n} +i\ALED\OM (\delta_{n',n+1}+\delta_{n',n-1}) -i \frac{\ALED}{2}(\omega+2n\OM)(\delta_{n',n+1}-\delta_{n',n-1})
    \end{align}




\section{Transmission through Archimedes' Screw}

In order to calculate transmission through a finite-length Archimedes' screw we expand the $E$ and $H$ fields into a Floquet-Bloch basis and calculate the wavevector eigenvalues $k$ for a fixed frequency of the impinging wave. Let us assume solutions of the form $\mathbf{E} = E_x \mathbf{x} + E_y \mathbf{y}$, $\mathbf{H} = H_x \mathbf{x} + H_y \mathbf{y}$, where the eigen-solutions for each field are of the form:
\begin{align}
    \psi &= e^{i(k x-\omega t)} \sum_n \Phi_n e^{2in(g x-\Omega t)}
\end{align}{}

In the absence of modulation, the vacuum eigenvalues read, for both $x$ and $y$ polarisations: 
\begin{align}
    k_{vn}^{\pm} = - 2ng \pm (\omega + 2n\Omega),
\end{align}{} and the corresponding eigenvectors are \begin{align}
    \mathbf{v_{n,x}^{\pm}} &= \begin{pmatrix} 0 \\ 1 \\ 1 \\ 0 \end{pmatrix}/\sqrt{2} &
    \mathbf{v_{n,y}^{\pm}} &= 
    \begin{pmatrix} 1 \\ 0 \\ 0 \\-1 
    \end{pmatrix}/\sqrt{2}
\end{align}{}.

Subsequently, we can write the fields at the left (1) and right (2) interfaces of the metamaterial as superpositions of the vacuum eigenvectors:
\begin{align}
    \begin{pmatrix} \mathbf{E_{v,z}^{(1)}} \\ \mathbf{H_{v,y}^{(1)}} \end{pmatrix} &= \mathbf{M^{vinc}} \mathbf{e_{vinc}^{(1)}} + \mathbf{M^{vref}}\mathbf{e_{vref}^{(1)}} & \begin{pmatrix} \mathbf{E_{v,z}^{(2)}} \\ \mathbf{H_{v,y}^{(2)}} \end{pmatrix} = \mathbf{M^{vinc}} \mathbf{e_{vtra}^{(2)}}
\end{align}{}
and inside the metamaterial:
\begin{align}
    \begin{pmatrix} \mathbf{E_{m,z}^{(1)}} \\ \mathbf{H_{m,y}^{(1)}} \end{pmatrix} &= \mathbf{M^{m}} \mathbf{e_{m}} & \begin{pmatrix} \mathbf{E_{m,z}^{(2)}} \\ \mathbf{H_{m,y}^{(2)}} \end{pmatrix} = \mathbf{M^{m}} \mathbf{P}\mathbf{e_{m}}
\end{align}{}
where $\mathbf{M^{vinc}}$ and $\mathbf{M^{vref}}$ are rectangular matrices containing the right- and left-propagating vacuum eigenvectors respectively, $\mathbf{M^m}$ is a square matrix containing all eigenvectors inside of the metamaterial, and $\mathbf{P}$ is a diagonal matrix $P_{mn} = \exp{(ik_m d)} \delta_{mn}$ which propagates each eigenvector from the left to the right interface. The vector $\mathbf{e_{v_inc}}$ contains the amplitudes of the input fields at the left interface, whereas the vectors $\mathbf{e_{v_ref}}$ and $\mathbf{e_{v_tra}}$ contain the unknown reflected and transmitted amplitudes respectively. Applying the continuity of $E_{z}$ and $H_y$ at the two interfaces, we arrive at a matrix equation:
\begin{align}
    \begin{pmatrix} \mathbf{A} & \mathbf{B} \end{pmatrix}\begin{pmatrix} \mathbf{e_{vtr}^{(2)}} \\ \mathbf{e_{vref}^{(2)}} \end{pmatrix} =     \begin{pmatrix} \mathbf{M^{vinc}} & 0 \end{pmatrix}\begin{pmatrix} \mathbf{e_{vinc}^{(1)}} \\ 0 \end{pmatrix} \label{eq:eq_mod_dielectric}
\end{align}{}
where:
\begin{align}
    \mathbf{A} &= \mathbf{M^m} (\mathbf{M}^m \mathbf{P})^{-1}\mathbf{M^{vinc}} & \mathbf{B} = - \mathbf{M^{mvref}}
\end{align}{}
are rectangular matrices, such that their concatenation is square, so that the transmitted and reflected amplitudes can be readily calculated by inverting Eq. \ref{eq:eq_mod_dielectric}.



\end{document}
