# Exact closed forms for the transmittance of electromagnetic waves in one-
dimensional anisotropic periodic media

J. C. Torres-Guzmán<EMAIL_ADDRESS>A. Díaz de Anda  J. Arriaga
Instituto de Física, Benemérita Universidad Autónoma de Puebla, P.O. Box J-48
72570, Puebla, Puebla, Mexico

###### Abstract

In this work, we obtain closed expressions for the transfer matrix and the
transmittance of electromagnetic waves propagating in finite 1D anisotropic
periodic stratified media with an arbitrary number of cells. By invoking the
Cayley-Hamilton theorem on the transfer matrix for the electromagnetic field
in a periodic stratified media formed by $N$ cells, we obtain a fourth-degree
recursive relation for the matrix coefficients that defines the so-called
Tetranacci Polynomials. In the symmetric case, corresponding to a unit-cell
transfer matrix with a characteristic polynomial where the coefficients of the
linear and cubic terms are equal, closed expressions for the solutions to the
recursive relation, known as symmetric Tetranacci Polynomials, have recently
been derived, allowing us to write the transfer matrix and transmittance in a
closed form. We show as sufficient conditions that the $4\times 4$
differential propagation matrix of each layer in the binary unit cell,
$\mathbb{\Delta}$, $a)$ has eigenvalues of the form $\pm~{}p_{1}$,
$\pm~{}p_{2}$, with $p_{1}\neq p_{2}$, and $b)$ its off-diagonal $2\times 2$
block matrices possess the same symmetric structure in both layers. Otherwise,
the recursive relations are still solvable for any $4\times 4$-matrix and
provide an algorithm to compute the $N$-th power of the transfer matrix
without carrying out explicitly the matrix multiplication of $N$ matrices. We
obtain analytical expressions for the dispersion relation and transmittance,
in closed form, for two finite periodic systems: the first one consists of two
birefringent uniaxial media with their optical axis perpendicular to the
$z$-axis, and the second consists of two isotropic media subject to an
external magnetic field oriented along the $z$-axis and exhibiting the Faraday
effect. Our formalism applies also to lossy media, magnetic anisotropy or
optical activity.

###### keywords:

Anisotropic Optical Media, Cayley-Hamilton Theorem, Transfer matrix method,
Tetranacci polynomials, Faraday effect.

††journal: Journal of Physics A: Mathematical and Theoretical

##  1 Introduction

Unlike travelling or standing waves, wave propagation in infinite periodic
media allows frequencies that fall into continuous bands, separated by
forbidden gaps. This fact, first noted by Kronig and Penney, is the foundation
for the modern theory of solids [1]. The band structure is at the core of
solid-state physics and has been studied in many different systems, quantum
and classical. The systems range from the nanometric scale in photonics to
oceanographic waves, covering a range in size of more than a dozen orders of
magnitude.

In the study of infinite periodic systems, the problem is dramatically
simplified thanks to Bloch’s theorem [2]. However, in practice only a few
systems can be considered infinite periodic. Therefore, finite periodic
systems are more difficult to analize, since Bloch’s theorem does not apply.
Fortunately, the finite periodic case, involving $2\times 2$-matrices, can be
solved analytically for an arbitrary number of unit cells, $N$, first
demonstrated in the optical case by Abèles [3] in 1950 and extended to the
quantum context [4] in 1981 and later by others [5, 6]. Thus, for one-
dimensional waves in locally periodic (finite periodic) media [7, 8, 9, 10,
11, 12, 13], analytical expressions for quantities like the transmittance are
available in terms of the Chevychev polynomials of the second kind [14].

Since then, light propagation in modern optical systems has become a central
issue. Thin films, periodic layered media, and other stratified systems
exhibit novel and attractive features in technological applications as anti-
reflection coatings [15, 16], x-ray mirrors [17, 18], transistors [19, 20],
thin-film photovoltaics [21, 22], photonic crystals [23], among others. In
particular, layered optical periodic media finds application in the
development of photonic topological insulators [24, 25, 26, 27, 28, 29, 30],
as they possess robust edge states, which persist in the presence of
perturbations and make them ideal for the transfer or storage of energy or
information. Even realizations of robust edge states in one-dimensional
structures [31, 32, 33, 34, 35], have been reported recently.

Many of the above active and passive optical devices require the epitaxial
growth of anisotropic thin films. Therefore, the design of these devices
relies on powerful mathematical and numerical techniques to adequately
describe light propagation in anisotropic layered media. In this direction and
from a theoretical point of view, the Abèles formalism was extended to a
$4\times 4$-matrix technique [36] to describe the propagation and reflection
by stratified anisotropic dielectric media, and later extended to include
magnetic anisotropy and optical activity [37]. In this formalism, the transfer
matrix of an anisotropic layer is expressed in terms of its eigenvalues and
eigenvectors. However, the solution of the associated eigenvalue problem is
left out. This eigenvalue problem is then solved and incorporated into a new
formalism to describe the electromagnetic propagation in non-magnetic
arbitrarily birefringent layered media [38, 39]. By taking advantage of
Bloch’s theorem and by carrying out the matrix multiplication of the whole
layered system, this new formalim allowed the study of the band structure and
Bloch waves in a periodic birefringent media. These works have led to the
development of several numerical methods to describe light propagation in
anisotropic stratified systems [40, 41, 42]. Recently, a different approach
has been derived [43] to study isotropic and anisotropic alternating layers
with one or two anisotropic defect layers in the framework of the temporal
coupled-mode theory. Despite these commendable efforts, closed-form analytical
expressions, e.g. for transmittance, are not yet available.

Among the pioneering efforts on the formulation of the transfer matrix in the
context of finite elastic periodic systems, we should mention the work of Lin
and McDaniel [44], who used the Cayley-Hamilton (C-H) theorem [7, 45, 46, 47],
to obtain an expression for the $4\times 4$-transfer matrix of $N$ cells in
the propagation of bending oscillations. Unfortunately, this expression is
invalid for degenerate eigenvalues of the unit-cell matrix [48]. We have
recently shown [49] by a similar procedure as that of Ref. [44], that the
$N$-th power of a $4\times 4$-matrix can be reduced to the calculation of a
set of four functions, recursively defined, instead of multiplying the
matrices, as was the original goal in Ref. [44]. Such recursive relations
define the Tetranacci Polynomials (TP) [50] and are solvable for an arbitrary
number of cells, $N$ [49]. Furthermore, we have shown that in free
oscillations of Timoshenko-Ehrenfest beams [49], the eigenvalues of the
transfer matrix corresponding to each layer composing the unit cell are non-
degenerate. Hence, we obtained the recursive relationship for symmetric TP,
which allows us to use the recently derived closed-form expressions obtained
in the context of a finite Kitaev chain [51, 52, 53] to write the transfer
matrix of $N$ cells in a closed-form. Moreover, these expressions account for
the degeneracy of the eigenvalues of the unit-cell transfer matrix [48]. In
this work, we obtain the sufficient conditions the optical $4\times
4$-differential propagation matrix of each layer in the unit cell must fulfil
to obtain closed forms, in terms of symmetric TP’s, for the transmittance in
finite 1D anisotropic periodic stratified media.

The structure of the paper is as follows: In Section 2, we review the
quadrivector formulation to describe the propagation of waves in optical
stratified media. In Section 3, we obtain a set of four recursive relations to
calculate the $N$-th power of the transfer matrix by using the C-H theorem. In
Section 4, we discuss the restrictions on the $4\times 4$ differential
propagation matrix to fulfil the symmetry condition that allows the use of
symmetric TP’s. In Section 5, we obtain the dispersion relation without using
explicitly Bloch’s theorem. In Section 6, we obtain closed forms for the
$N$-th power of the transfer matrix and the transmittance of the finite
periodic system. In Section 7, we give our conclusions.

##  2 $4\times 4$-Matrix formulation

We are interested in a locally periodic anisotropic media consisting of $N$
cells. First, we review the electro-magnetic field in each layer. The specific
nature of each layer can be different in the more general case. In special
cases, exact results can be obtained. Thus, a constitutive relation between
the vectors $\mathbf{C}=(D_{x},D_{y},D_{z},B_{x},B_{y},B_{z})$ and
$\mathbf{G}=(E_{x},E_{y},E_{z},H_{x},H_{y},H_{z})$, where $E_{j}$, $D_{j}$,
$B_{j}$ and $H_{j}$, $j=x,y,z,$ are the cartesian components of the
electromagnetic-field vectors, can be generally expressed as

$\mathbf{C}=\mathbf{M}\mathbf{G},$ (1)

where the matrix $\mathbf{M}$, is the so-called optical matrix and carries all
the information about the anisotropic optical properties of the medium [37].
The optical matrix can be partitioned as

$\mathbf{M}=\left(\begin{array}[]{cc}\mathbf{\epsilon}&\mathbf{\rho}\\\
\mathbf{\rho}^{\prime}&\mathbf{\mu}\end{array}\right),$ (2)

where as usual, $\mathbf{\epsilon}$ and $\mathbf{\mu}$ are the permitivity and
permeability tensors, respectively, and $\mathbf{\rho}$ and
$\mathbf{\rho}^{\prime}$ are the optical rotation tensors [37]. Consider the
propagation of a monochromatic plane wave of frequency $\omega$ obliquely
incident int the $xz$-plane from an isotropic medium ($z<0$) onto a periodic
anisotropic layered media of $N$ cells ($z>0$) stratified along the $z$-axis,
in which $\mathbf{M}$ is a function only of $z$. In this case, the
$x$-component of the propagation vector, $k_{x}$, is constant and there is no
$y$ component in any layer of the structure,

$k_{x}=\frac{\omega}{c}n_{0}\sin\theta_{0},$ (3)

where $c$, $n_{0}$, and $\theta_{0}$ are the free-space wave velocity, the
refractive index of the isotropic ambient and the angle of incidence,
respectively. Then, the Maxwell’s equations in an anisotropic layer can be
cast into a $4\times 4$ matrix form [37]

$\frac{\partial\bm{\Psi}}{\partial z}=i\omega\mathbb{\Delta}\bm{\Psi},$ (4)

where

$\bm{\Psi}=\left(\begin{array}[]{c}E_{x}\\\ H_{y}\\\ E_{y}\\\
-H_{x}\end{array}\right),$ (5)

is the generalized quadrivector field, and the $4\times 4$-matrix
$\mathbb{\Delta}$ is function of the elements of the $6\times 6$-optical
matrix $\mathbf{M}$. The details to construct $\mathbb{~{}}\Delta$ from
$\mathbf{M}$ can be found in Refs. [37, 54]. In the general case of a
stratified anisotropic structure, $\mathbf{M}$ is some arbitrary function of
$z$ and the Eq. (4) does not, in general, have an analytical solution. In the
special case when $\mathbf{M}$ is independent of $z$, Eq. (4) is directly
integrable to yield

$\bm{\Psi}(z+\ell)=\mathbf{T}(\ell)\bm{\Psi}(z),$ (6)

where $\mathbf{T}(\ell)=e^{i\omega\ell\mathbb{\Delta}}$, which can be expanded
in a power series and the sum carried out analytically in a few cases only.
Analytical expressions for $\mathbf{T}(\ell)$ can be determined by analyzing
the four plane-wave solutions of Eq. (4). We address this case in Section 4.
Here, we pause the $4\times 4$-matrix formulation to introduce the C-H
theorem, upon which our formalism relies to obtain the transfer matrix of an
anisotropic periodic media.

##  3 The Cayley-Hamilton theorem and the Tetranacci polynomials recurrence
relation

Consider a locally periodic anisotropic media consisting of $N$ cells. Each
unitary cell consists of two anisotropic layers of thickness $\ell_{1}$ and
$\ell_{2}$, respectively. Let us choose the unitary cell such that the
frontiers are at $z=z_{0}$ and $z=z_{2}$, and the interface between the two
layers is at $z=z_{1}$, where $z_{0}<z_{1}<z_{2}$. Then the EM-quadrivector,
$\bm{\Psi}$, at $z=z_{1}$ can be written in terms of its value at $z=z_{0}$,

$\bm{\Psi}(z_{1})=\mathbf{T}_{1}(\ell_{1})\bm{\Psi}(z_{0}),$ (7)

where $\mathbf{T}_{1}(\ell_{1})$ is the transfer matrix for the first layer.
Similarly,

$\bm{\Psi}(z_{2})=\mathbf{T}_{2}(\ell_{2})\bm{\Psi}(z_{1}).$ (8)

At $z=z_{1}$, we demand continuity in the parallel components of the EM field,
i.e. continuity on $\bm{\Psi}$. Thus, the transfer matrix for the unitary
cell, $\mathbf{T}_{c}$, is given by

$\mathbf{T}_{c}=\mathbf{T}_{2}(\ell_{2})\mathbf{T}_{1}(\ell_{1}).$ (9)

The transfer matrix for $N$ cells is then given by $\mathbf{T}_{c}^{N}$. Now
we invoke the C-H theorem [7, 46, 47] to calculate $\mathbf{T}_{c}^{N}$. It
states that any matrix, $\mathbf{T}_{c}$, obeys its own characteristic
polynomial equation

$\mathbf{T}_{c}^{4}-\eta\mathbf{T}_{c}^{3}+\zeta\mathbf{T}_{c}^{2}-\tau\mathbf{T}_{c}+\delta_{0}\mathbf{I}=\mathbb{0},$
(10)

where

$\displaystyle\eta$ $\displaystyle=$ $\displaystyle\mbox{Tr}(\mathbf{T}_{c}),$
(11) $\displaystyle\zeta$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\eta^{2}-\mbox{Tr}(\mathbf{T}_{c}^{2})\right),$
(12) $\displaystyle\tau$ $\displaystyle=$
$\displaystyle\frac{1}{6}\left(\eta^{3}-3\eta\mbox{Tr}(\mathbf{T}_{c}^{2})+2\mbox{Tr}(\mathbf{T}_{c}^{3})\right),$
(13) $\displaystyle\delta_{0}$ $\displaystyle=$
$\displaystyle\mbox{det}(\mathbf{T}_{c}),$ (14)

where Tr and det denote the trace and the determinant of a matrix,
respectively. $\mathbf{I}$ and $\mathbb{0}$ are the identity and null $4\times
4$-matrices, respectively. As in Refs. [44, 49], this equation can be then
used to calculate any power, $N$, of $\mathbf{T}_{c}$

$\mathbf{T}_{c}^{N}=\alpha(N-1)\mathbf{T}_{c}^{3}-\beta(N-2)\mathbf{T}_{c}^{2}+\gamma(N-3)\mathbf{T}_{c}-\delta(N-4)\mathbf{I},\quad\mbox{for
}N\geq 5.$ (15)

where the functions $\alpha(n)$, $\beta(n),$ $\gamma(n)$, and $\delta(n)$, are
recursively defined by [49],

$\displaystyle\alpha(n)$ $\displaystyle=$
$\displaystyle\eta\alpha(n-1)-\beta(n-2),$ (16) $\displaystyle\beta(n-1)$
$\displaystyle=$ $\displaystyle\zeta\alpha(n-1)-\gamma(n-3),$ (17)
$\displaystyle\gamma(n-2)$ $\displaystyle=$
$\displaystyle\tau\alpha(n-1)-\delta(n-4),$ (18) $\displaystyle\delta(n-3)$
$\displaystyle=$ $\displaystyle\delta_{0}\alpha(n-1).$ (19)

Equation (15) can be demonstrated by mathematical induction. Eqs. (16)-(19)
are valid for $n\geq 4$ with initial values $\alpha(3)=\eta$,
$\beta(2)=\zeta$, $\gamma(1)=\tau$ and $\delta(0)=\delta_{0}$, corresponding
to the case $N=4$. The powers, $N=2,3$, have to be done by matrix
multiplication, but this is computationally cheap. Substitution of Eqs.
(17)-(19) into Eq. (16) yields an Eq. for $\alpha(n)$ solely [49]

$\alpha(n)=\eta\alpha(n-1)-\zeta\alpha(n-2)+\tau\alpha(n-3)-\delta_{0}\alpha(n-4),$
(20)

valid for $n\geq 7$. Such recurrence relationship defines the Tetranacci
polynomials [50], with coefficients $\eta,\zeta,\tau$, and $\delta_{0}$. With
the initial values, $\alpha(3)$, $\beta(2)$, $\gamma(1)$ and $\delta(0)$, and
Eqs. (16)-(19), we can obtain the other seminal values, $\alpha(j)$,
$j=4,5,6,$ to iterate Eq. (20) and obtain $\alpha(N-1)$. Then, the values
$\alpha(n)$, with $n=N-1,N-2,N-3,N-4$ can be substituted back into Eqs.
(16)-(19) once again, to obtain the coefficients $\alpha(N-1)$, $\beta(N-2)$,
$\gamma(N-3)$, and $\delta(N-4)$, required in Eq. (15) to finally obtain
$\mathbf{T}_{c}^{N}$. Therefore, the computational cost is significantly
reduced since no matrix multiplication is involved beyond the one implicit in
$\mathbf{T}_{c}^{2}$ and $\mathbf{T}_{c}^{3}$. Also, the recurrence
relationship of Eq. (20) is valid even for degenerate eigenvalues of the
transfer matrix $\mathbf{T}_{c}$ [48]. This iterative algorithm is one of the
methods implemented in this work and we will refer to it as Method A [49].

##  4 The traceless matrix condition

Further simplifications to calculate the $N$-th power of the transfer matrix
$\mathbf{T}$ can be achieved if the matrix, $\mathbb{\Delta}$, in addition to
being constant and independent of $z$, has a characteristic polynomial that is
an even function [46, 47, 49]. In other words, if $p$ is an eigenvalue of
$\mathbb{\Delta}$, then $-p$ is also an eigenvalue. Hence, its eingenvalues
are $(p_{1},p_{2},-p_{1},-p_{2})$, with $p_{1}\neq~{}p_{2}$. Although this
condition seems arbitrary and very restrictive, it is not. This condition is
consistent with having identical waves propagating in both directions $\pm z$,
a situation encountered in many physical systems beyond the particular cases
addressed here. Consider the ansatz for the solution of Eq. (4),

$\bm{\Psi}(z)=\bm{\Psi}_{\ell}(0)e^{ik_{\ell}z},\quad\ell=1,2,3,4,$ (21)

that accounts for four distinct plane-wave solutions. The substitution of Eq.
(21) into Eq. (4) yields the matrix eigenvalue equation,

$(k_{z}\mathbf{I}-\omega\mathbb{\Delta})\bm{\Psi}(0)=0,$ (22)

whose eigenvalues, $k_{\ell}$ (we have dropped the subindex $z$ to simplify
nomenclature), are proportional to the roots, $p_{\ell}$, of the
characteristic polinomial of $\mathbb{\Delta}$, with $\omega$ the
proportionality constant, since $k_{\ell}$ may be analytically obtained by
solving the quartic polynomial equation,

$\mbox{det}[k_{z}\mathbf{I}-\omega\mathbb{\Delta}]=0.$ (23)

Therefore $k_{3}=-k_{1}$ and $k_{4}=-k_{2}$, corresponding to two waves
propagating in both directions $\pm z$.

According to the discussion above, the even parity of the characteristic
polynomial implies that both its linear and cubic coefficients,
$\eta_{\mathbb{\Delta}}$ and $\tau_{\mathbb{\Delta}}$, respectively, vanish,
i.e. $\eta_{\mathbb{\Delta}}=\tau_{\mathbb{\Delta}}=0$. In terms of traces,
according to Eqs. (11)-(13) applied to $\mathbb{\Delta}$, we have Tr
$\mathbb{\Delta}=\mbox{ Tr}~{}\mathbb{\Delta}^{3}=0$, i.e. both
$\mathbb{\Delta}$ and $\mathbb{\Delta}^{3}$ matrices are traceless [46, 47,
49]. The corresponding non-degeneracy condition, $p_{1}\neq~{}p_{2}$, has two
implications: the first one is that the independent term, $\mbox{det
}\mathbb{\Delta}$, must be different from zero, i.e. $|\mathbb{\Delta}|\neq
0$, to avoid the trivial and degenerate solution, $p=0$. The second
implication is that according to Eq. (12) applied to $\mathbb{\Delta}$,
$(\mbox{Tr~{}}\mathbb{\Delta}^{2})^{2}\neq 16\mbox{ det~{}}\mathbb{\Delta}$.

We show in Appendix A, that two types of $\mathbb{\Delta}$-matrices can be
distinguished for optical systems that fulfill the above conditions according
to the symmetry of its off-diagonal $2\times 2$ block matrices. If the two
constituents of the unit cell have the same symmetric structure, then we have
$\eta=\tau$, i. e. the characteristic polynomial, Eq. (10), is symmetric. We
demonstrate this in Appendix B, by noticing from the reciprocal polynomial
that if $\eta=\tau$, then
$\mbox{Tr~{}}\mathbf{T}_{c}=\mbox{Tr~{}}\mathbf{T}_{c}^{-1}$ [49]. Under the
category of type I, we can identify the case of the uniaxial orthorhombic
crystal with its optical axis perpendicular to the axis of propagation and
rotated from the plane of incidence [54, 55], or the Faraday rotator of light
subjected to a magnetic field in the $z$ direction [37, 56]. Moreover, under
the category of type II, we can find the Drude model for an isotropic
optically active medium [37, 57], or other similar models for the same
phenomenon [58, 59, 60]. Unfortunately, as we show in Appendix B, if the
constituents of the unit cell have different symmetric structure, then
$\eta\neq\tau$. Even so, we can have closed expressions for the transmittance
of periodic optical systems consisting of uniaxial crystals, Faraday rotators
or combinations of these. The same applies for periodic optical systems
consisting of layers with natural optical activity.

The general solution for $\bm{\Psi}$ is a linear combination of the
eigenvectors, $\vec{V}_{j}$, $j=1,2,3,4$, contained as columns in the matrix
$\mathbb{V}=(\vec{V}_{1},\vec{V}_{2},\vec{V}_{3},\vec{V}_{4})$,

$\bm{\Psi}(z)=\mathbb{V}\mathbb{\lambda}(z)\vec{A},$ (24)

where

$\mathbb{\lambda}(z)=\mbox{diag}(e^{ik_{1}z},e^{ik_{2}z},e^{-ik_{1}z},e^{-ik_{2}z}),$
(25)

i.e., $\mathbb{\lambda}(z)$ is a diagonal matrix, and
$\vec{A}=(A_{1},A_{2},A_{3},A_{4})^{T}$, where $A_{j}$, $j=1,...,4$, are
arbitrary constants. If $\bm{\Psi}$ is known at $z=z_{0}$ in a layer, then the
value of $\bm{\Psi}$ at any $z$ value in the same layer, can be obtained by
solving for $\vec{A}$ from
$\bm{\Psi}(z_{0})=\mathbb{V}\mathbb{\lambda}(z_{0})\vec{A}$ and then
substituting into Eq. (24) to obtain

$\bm{\Psi}(z)=\mathbb{V}\mathbb{\lambda}(z-z_{0})\mathbb{V}^{-1}\bm{\Psi}(z_{0}),$
(26)

which provides the transfer matrix, $\mathbf{T}$, for the quadrivector
$\bm{\Psi}$ in the layer in question [37],

$\mathbf{T}(z-z_{0})=\mathbb{V}\mathbb{\lambda}(z-z_{0})\mathbb{V}^{-1}.$ (27)

From Eqs. (25) and (27), we also obtain $\mbox{det}(\mathbf{T})=1$ [63]. We
emphasize that Eq. (27) is valid only if $k_{1}\neq~{}k_{2}$, otherwise
$\mathbb{\Delta}$ might not have four independent eigenvectors. It is
precisely for isotropic media when the case, $k_{1}=~{}k_{2}$, might occur,
and therefore we exclude this possibility. However, from a practical point of
view, we can recover the case of an isotropic layer to any degree of
approximation by taking the anisotropy arbitrarily small, such that the above
expression is still valid to avoid discontinuous solutions [42, 61, 62].

Notice also, the transfer matrix, $\mathbf{T}$, is the exponential of the
$\mathbb{\Delta}$ matrix [63]

$\mathbf{T}(z-z_{0})=e^{i\omega\mathbb{\Delta}(z-z_{0})},$ (28)

and therefore, the transfer matrix of the unit cell, $\mathbf{T}_{c}$, is then
given by

$\mathbf{T}_{c}(\ell)=e^{i\omega\mathbb{\Delta}_{2}\ell_{2}}e^{i\omega\mathbb{\Delta}_{1}\ell_{1}},$
(29)

which is equal to the exponential of the sum of the matrix exponents only if
$\mathbb{\Delta}_{1}$ and $\mathbb{\Delta}_{2}$ commute. However, this
condition is generally not satisfied and then the transfer matrix,
$\mathbf{T}_{c}$, can be computed by means of the Lie product formula [63] and
taking a large number of terms to approximate it. Alternatively, if
$\mathbb{\Delta}_{1}$ and $\mathbb{\Delta}_{2}$ are sufficiently small
matrices, then $\mathbf{T}_{c}=e^{i\omega\mathbf{Q}}$, where $\mathbf{Q}$ may
be computed as a series in commutators of $\ell_{1}\mathbb{\Delta}_{1}$ and
$\ell_{2}\mathbb{\Delta}_{2}$ by using the Baker-Campbell-Hausdorff formula
[64]. This approach is used in Ref. [65] to investigate the effects of
chirality in the Casimir torque.

##  5 The dispersion relation

When the symmetric condition $\eta=\tau$ is substituted into Eq. (20), the
recursion relation for the symmetric Tetranacci polynomials is obtained,

$\alpha(n)=\eta\alpha(n-1)-\zeta\alpha(n-2)+\eta\alpha(n-3)-\alpha(n-4),\quad
n\in\mathbb{Z},$ (30)

in terms of its initial values $\alpha(j)$, $j=-2,-1,0,1$, and complex
coefficients, $\eta$ and $\zeta$. They can be expressed as [51, 52, 53]

$\alpha(n)=Ae^{i\phi_{1}n}+Be^{-i\phi_{1}n}+Ce^{i\phi_{2}n}+De^{-i\phi_{2}n},$
(31)

where

$\cos(\phi_{\nu})=s_{\nu}/2,\quad\nu=1,2,$ (32)
$s_{\nu}=\frac{\eta\pm\sqrt{\eta^{2}-4(\zeta-2)}}{2},$ (33)

and provided $s_{\nu}^{2}\neq 4$ and $s_{1}\neq s_{2}$.

The coefficients $A,B,C$ and $D$, are given by the initial values $\alpha(j)$.
In Ref. [53], the values $j=-2,-1,0,1$, are used as seminal values while in
this work we use $j=3,4,5,6$, since we can always define $\alpha(n+5)=\xi(n)$,
with both $\alpha$ and $\xi$, satisfiying Eq. (30).

Notice that by substituting the ansatz, $\alpha(n)\propto w^{n}$, into Eq.
(30), we arrive at the characteristic equation, (see Eq. (10)), with
$\eta=\tau$,

$w^{4}-\eta w^{3}+\zeta w^{2}-\eta w+1=0.$ (34)

Hence, the eigenvalues of $\mathbf{T}_{c}$, are

$w_{\mu}=e^{\pm i\phi_{1}},e^{\pm i\phi_{2}},\quad\mu=1,...,4,$ (35)

where $\phi_{1}$ and $\phi_{2}$ can be complex quantities. They correspond to
two waves propagating in both directions $\pm\hat{z}$ with phases $\phi_{1}$
and $\phi_{2}$, respectively.

We consider two examples of binary periodic optical system of type I only, we
will investigate type II binary cells in a future work. The first one has a
unit cell composed of two birefringent uniaxial media and the second has a
unit cell composed by two isotropic media subject to an external magnetic
field oriented along the $z$-direction, exhibiting Faraday effect. For
simplicity, we consider lossless media. Under these assumptions, it turns out
that $\eta$ and $\zeta$ are real valued [49] for both systems, then $s_{\nu}$,
$\nu=1,2$, are also real for $\eta^{2}\geq 4(\zeta-2)$, and become complex if
$\eta^{2}<4(\zeta-2)$. Notice on the one hand, that in the regime of real
$s_{\nu}$, if $|s_{\nu}|\leq 2$, then $\phi_{\nu}$ is a real quantity and
describes propagating oscillations in both directions $\pm\hat{z}$; and if
$|s_{\nu}|>2$, then $\phi_{\nu}$ is pure imaginary and describes exponentially
decaying or increasing oscillations. On the other hand, in the regime of
complex values of $s_{\nu}$, we obtain from Eq. (33), $s_{2}=s_{1}^{*}$,
therefore $\phi_{2}=\pm\phi_{1}^{*}$. Also, according to Eq. (35), both
exponentials, $e^{\pm i\phi_{2}}$, are eigenvalues of $\mathbf{T}_{c}$,
therefore we can arbitrarily set $\phi_{2}=\phi_{1}^{*}$.

Figure 1: Dispersion relation for a unit cell composed by two birefringent
uniaxial media with their optical axis perpendicular to the $z$-axis, with
ordinary refraction indexes, $n_{o}^{(1)}=1.6$ and $n_{o}^{(2)}=1.1$, and
extraordinary refraction indexes $n_{e}^{(1)}=1.9$ and $n_{e}^{(2)}=1.4$. The
corresponding filling fractions are $f_{1}=0.4$ and $f_{2}=0.6$, respectively.
We have assumed both the isotropic ambient and the substrate to have the same
refraction index, $n_{0}=n_{2}=1$ and for simplicity, normal incidence is
assumed, $\theta_{0}=0$. We have set the angles of the optical axes with
respect to $x$-axis, $\psi_{1}$ and $\psi_{2}$, of the first and second layer,
$\psi_{1}=0$ and $\psi_{2}=\pi/4$, respectively, to allow the so-called mode
coupling.

In Figure 1, we show the dispersion relationship, $\phi_{\nu}=k_{\nu}\ell$,
$\nu=1,2$, where $\ell$, is the unit cell length, $\ell=\ell_{1}+\ell_{2}$,
calculated from Eqs. (32), (33) and (35), for a binary cell consisting of
birefringent uniaxial media, with their optical axis perpendicular to the
$z$-axis. By convention, the imaginary part of $k_{\nu}\ell$ is represented by
negative abscissa values while the real part is represented by positive ones.
The interval of frequencies where $|s_{\nu}|\leq 2$, either for $\nu=1$ or
$\nu=2$ or both, defines the bands observed in Figure 1 ($k_{\nu}\ell>0$).
Similarly, the interval of frequencies where $|s_{\nu}|>2$, for both
$\nu=1,2$, defines the gap observed in Figure 1 ($k_{\nu}\ell<0$) around the
value $\omega^{\prime}=0.38$, where $\omega^{\prime}$ is the normalized
frequency $\omega^{\prime}=\frac{\omega\ell}{2\pi c}$. Around
$\omega^{\prime}=0.7$ and just above $\omega^{\prime}=1$, both branches shown
in Figure 1, overlap. Furthermore, the corresponding wavenumbers $k_{1}$ and
$k_{2}$ are complex, which corresponds to attenuating or increasing
propagating waves, and seem to be degenerate, i. e. $k_{1}=k_{2}$. This
scenario corresponds to the regime of complex $s_{\nu}$ and according to the
discussion in the above paragraph, $k_{2}=k_{1}^{*}$. Thus, despite the
imaginary part of both $k_{1}$ and $k_{2}$ have opposite signs, they are
plotted with the same sign in Figure 2. From Eq. (33), the edges of these
complex bands, occur when $\eta^{2}=4(\zeta-2)$, clearly observable in Figure
1, where the wavenumbers also become degenerate. The complex band around
$\omega^{\prime}=0.7$, connects two typical bands, as well as does the complex
band just above $\omega^{\prime}=1$. Finally, notice that as a thumb rule,
when one of the wavenumbers has an imaginary part in a particular interval of
frequency, the maximum value of its imaginary part is proportional to the size
of the interval of frequency. This thumb rule can be helpful in understanding
the partial transmission, as we will see in the next section.

Figure 2: Dispersion relation for a unit cell composed by two isotropic media
that exhibit the Faraday effect when an external magnetic field is applied
along the $z$-axis, with Faraday rotation parameters,
$\gamma_{1}/\epsilon_{v}=0.36$ and $\gamma_{2}/\epsilon_{v}=0.001$,
respectively. The corresponding refraction indexes of the first and second
layer are, $n_{o}^{(1)}=1.47$ and $n_{o}^{(2)}=1.7$, with filling fractions
$f_{1}=0.4$ and $f_{2}=0.6$, respectively. Again, we have assumed both the
isotropic ambient and the substrate to have the same refraction index,
$n_{0}=n_{2}=1$, and for simplicity, normal incidence is assumed,
$\theta_{0}=0$.

In Figure 2, we show the dispersion relation for our second example of binary
cell, consisting of istropic media, that exhibit the Faraday effect when an
external magnetic field is applied along the $z$-axis. A similar structure is
observed as in our first example. However, a closer look reveals the absence
of complex bands for the set of parameters chosen and for the range of
frequencies shown. Notice while in our first example, the complex bands
connected two typical bands, now these typical bands are just joint at
frequencies where $\eta^{2}=4(\zeta-2)$, i. e. the wavenumbers become
degenerate, $k_{1}=k_{2}$, at these frequencies.

##  6 Reflectance and transmittance

We can use the closed forms for the symmetric Tetranacci polynomials in Eq.
(30) [51, 52, 53],

$\begin{array}[]{cc}\xi(n)=&\frac{1}{s_{1}-s_{2}}\left[\varphi_{2}(n)\left(\xi_{-2}-s_{1}\xi_{-1}+\xi_{0}\right)-\varphi_{1}(n)\left(\xi_{-2}-s_{2}\xi_{-1}+\xi_{0}\right)\right.\\\
&\left.+\varphi_{1}(n+1)\left(\xi_{-1}-s_{2}\xi_{0}+\xi_{1}\right)-\varphi_{2}(n+1)\left(\xi_{-1}-s_{1}\xi_{0}+\xi_{1}\right)\right],\end{array}$
(36)

without degeneration, $s_{1}\neq s_{2}$, to calculate any power of the
transfer matrix $\mathbf{T}_{c}$, where we have employed the notation of Ref.
[53], with $\xi(j)=\xi_{j}$, $j=-2,-1,0,1$. $\varphi_{\nu}(n)$, $\nu=1,2,$
denotes the generalized Fibonacci polynomial [66, 67, 68], defined by,

$\varphi_{\nu}(n+1)=s_{\nu}\varphi_{\nu}(n)-\varphi_{\nu}(n-1),\quad
n\in\mathbb{Z},$ (37)

with $s_{\nu}$ given by Eq. (33) and initial values $\varphi_{\nu}(0)=0$,
$\varphi_{\nu}(1)=1$. By defining $\alpha(n+5)=\xi(n)$, we obtain
$\alpha(3)=\xi_{-2}$, $\alpha(4)=\xi_{-1}$, $\alpha(5)=\xi_{0}$ and
$\alpha(6)=\xi_{1}$. Then we can substitute the values $n=N-6$, $n=N-7$,
$n=N-8$, and $n=N-9$, in Eq. (36), to obtain $\alpha(n)$ for $n=N-1,N-2,N-3$
and $N-4$. The susbstitution of these $\alpha$-values in Eq. (17)-(19) yields
the coefficients $\beta(N-2)$, $\gamma(N-3)$ and $\delta(N-4)$, and therefore,
the $N$-power of $\mathbf{T}_{c}$. We refer to the method that uses the closed
forms of TP of Eq. (36) as Method B. Finally, we can calculate the
transmittance or reflectance as usual [54], from the relation,

$\Psi_{t}=\mathbf{T}_{c}^{N}(\Psi_{i}+\Psi_{r})$ (38)

where,

$\begin{array}[]{rcl}\Psi_{i}&=&\left(E_{ip}\cos\theta_{0},\frac{E_{ip}}{Z_{0}},E_{is},\frac{E_{is}}{Z_{0}}\cos\theta_{0}\right)^{T},\\\
\Psi_{r}&=&\left(-E_{rp}\cos\theta_{0},\frac{E_{rp}}{Z_{0}},E_{rs},-\frac{E_{rs}}{Z_{0}}\cos\theta_{0}\right)^{T},\\\
\Psi_{t}&=&\left(E_{tp}\cos\theta_{2},\frac{E_{tp}}{Z_{2}},E_{ts},\frac{E_{ts}}{Z_{2}}\cos\theta_{0}\right)^{T},\end{array}$
(39)

where $E_{\mu\nu}$, denotes the amplitude of the incident ($\mu=i$), reflected
($\mu=r$), and transmitted electric field for $\nu$-polarization, $\nu=p,s$.
$\theta_{2}$ is the transmission angle towards the substrate, $Z_{0}$ and
$Z_{2}$, are the isotropic ambient and substrate impedances, respectively. $T$
means transpose. The substitution of Eqs. (39) into Eq. (38) yields four
linear equations in the six components, $E_{\mu\nu}$, ($\mu=i,r,t$;
$\nu=s,p$). $E_{ts}$ can be eliminated from two of these equations while
$E_{tp}$ can in turn be eliminated from the remaining two equations. This
procedure yields two algebraic linear equations for the incident and reflected
fields alone, which can be written as,

$\bm{E}_{r}=\mathbf{r}\bm{E}_{i},$ (40)

where $\bm{E}_{\mu}=(E_{\mu p},E_{\mu s})^{T}$, $\mu=i,r,t$ and $\mathbf{r}$
is the $2\times 2$ complex-amplitude reflection matrix

$\mathbf{r}=\left(\begin{array}[]{cc}r_{pp}&r_{ps}\\\
r_{sp}&r_{ss}\end{array}\right),$ (41)

and is given by

$\mathbf{r}\\!=\\!(\mathcal{P}_{r}^{s}\mathcal{Q}_{r}^{p}-\mathcal{P}_{r}^{p}\mathcal{Q}_{r}^{s})^{-1}\left(\\!\begin{array}[]{cc}\mathcal{P}_{i}^{p}\mathcal{Q}_{r}^{s}-\mathcal{P}_{r}^{s}\mathcal{Q}_{i}^{p}&\mathcal{P}_{i}^{s}\mathcal{Q}_{r}^{s}-\mathcal{P}_{r}^{s}\mathcal{Q}_{i}^{s}\\\
\mathcal{P}_{r}^{p}\mathcal{Q}_{i}^{p}-\mathcal{P}_{i}^{p}\mathcal{Q}_{r}^{p}&\mathcal{P}_{r}^{p}\mathcal{Q}_{i}^{s}-\mathcal{P}_{i}^{s}\mathcal{Q}_{r}^{p}\\\
\end{array}\\!\right)$ (42)

where

$\begin{array}[]{ccl}\mathcal{P}_{\mu}^{p}&=&\pm\cos\theta_{0}(T_{11}-T_{21}Z_{2}\cos\theta_{2})+(T_{12}-T_{22}Z_{2}\cos\theta_{2})/Z_{0},\\\
\mathcal{P}_{\mu}^{s}&=&\pm\frac{\cos\theta_{0}}{Z_{0}}(T_{14}-T_{24}Z_{2}\cos\theta_{2})+T_{13}-T_{23}Z_{2}\cos\theta_{2},\\\
\mathcal{Q}_{\mu}^{p}&=&\pm\cos\theta_{0}(T_{31}\frac{\cos\theta_{2}}{Z_{2}}-T_{41})+(T_{32}\frac{\cos\theta_{2}}{Z_{2}}-T_{42})/Z_{0},\\\
\mathcal{Q}_{\mu}^{s}&=&\pm\frac{\cos\theta_{0}}{Z_{0}}(T_{34}\frac{\cos\theta_{2}}{Z_{2}}-T_{44})+T_{33}\frac{\cos\theta_{2}}{Z_{2}}-T_{43},\end{array}$
(43)

$T_{mn}$ are the matrix elements of the transfer matrix $\mathbf{T}_{c}^{N}$.
In the above equations, the $\mu=i$ and $\mu=r$ values on the left correspond
to the upper (+) and lower (-) signs on the right, respectively. Finally, the
reflectance, $R_{\mu}$, for incident $\mu$-polarized light is given by

$R_{\mu}=|r_{p\mu}|^{2}+|r_{s\mu}|^{2},\quad\mu=p,s.$ (44)

Substitution of $\mathbf{T}_{c}^{N}$ into this equation yield a closed form
for the reflectance. Similarly, by substituting Eq. (40) into the second and
third component of $\Psi_{t}$ of Eq. (38), the reflected field components are
eliminated, yielding two linear equations for the transmitted and incident
fields $\mathbf{E}_{t}$ and $\mathbf{E}_{i}$ only,

$\bm{E}_{t}=\mathbf{t}\bm{E}_{i},$ (45)

where $\mathbf{t}$ is the $2\times 2$ complex-amplitude transmission matrix

$\mathbf{t}=\left(\begin{array}[]{cc}t_{pp}&t_{ps}\\\
t_{sp}&t_{ss}\end{array}\right),$ (46)

where

$\begin{array}[]{ccl}t_{pp}&=&Z_{2}\left[(T_{21}\cos\theta_{0}+\frac{T_{22}}{Z_{0}})(1-r_{pp})+r_{sp}(T_{23}-T_{24}\frac{\cos\theta_{0}}{Z_{0}})\right],\\\
t_{ps}&=&Z_{2}\left[T_{23}+T_{24}\frac{\cos\theta_{0}}{Z_{0}}+r_{ss}(T_{23}-T_{24}\frac{\cos\theta_{0}}{Z_{0}})+r_{ps}(\frac{T_{22}}{Z_{0}}-T_{21}\cos\theta_{0})\right],\\\
t_{sp}&=&T_{31}\cos\theta_{0}+\frac{T_{32}}{Z_{0}}+r_{pp}(-T_{31}\cos\theta_{0}+\frac{T_{32}}{Z_{0}})+r_{sp}(T_{33}-T_{34}\frac{\cos\theta_{0}}{Z_{0}}),\\\
t_{ss}&=&T_{33}+T_{34}\frac{\cos\theta_{0}}{Z_{0}}+r_{ss}(T_{33}-T_{34}\frac{\cos\theta_{0}}{Z_{0}})+r_{ps}(\frac{T_{32}}{Z_{0}}-T_{31}\cos\theta_{0}).\\\
\end{array}$ (47)

Therefore, the transmittance, $T_{\mu}$, for incident $\mu$-polarized light is
given by

$T_{\mu}=\frac{Z_{0}\cos\theta_{2}}{Z_{2}\cos\theta_{0}}\left(|t_{p\mu}|^{2}+|t_{s\mu}|^{2}\right),\quad\mu=p,s.$
(48)

For $\eta$ and $\zeta$ real [49], the condition $s_{1}=s_{2}$, occurs at the
edges of the complex bands, i.e. when $\eta^{2}=4(\zeta-2)$. Thus, the
corresponding expressions for $\xi(n)$ in the degenerate case, $s_{1}=s_{2}$,
with either $s_{1}^{2}\neq 4$ or $s_{1}^{2}=4$, are required. They can be
obtained [49] from preliminary steps shown in [53], and we reproduce them
here: If $s_{1}=s_{2}$ but $s_{1}^{2}\neq 4$, then

$\begin{array}[]{cl}\xi(n)=&\frac{1}{s_{1}^{2}-4}\left[\varphi_{1}(n+2)\left(n\xi_{1}-3(n+1)s_{1}\xi_{0}\right)\right.\\\
&+\varphi_{1}(n+1)\left((1-n)\xi_{-2}+2(n+2)(s_{1}^{2}-1)\xi_{0}\right)\\\
&+\varphi_{1}(n)\left(-(n+2)\xi_{1}+2(n-1)(s_{1}^{2}-1)\xi_{-1}\right)\\\
&\left.+\varphi_{1}(n-1)\left((n+1)\xi_{-2}-3ns_{1}\xi_{-1}\right)\right],\end{array}$
(49)

if $s_{1}=s_{2}$ and $s_{1}^{2}=4$, then

$\begin{array}[]{cc}\xi(n)=&\frac{s_{1}}{12}\left[(n+1)(n+3)\xi_{0}\varphi_{1}(n+2)+n(n+2)(\xi_{1}-2s_{1}\xi_{0})\varphi_{1}(n+1)\right.\\\
&\left.+(n^{2}-1)(2s_{1}\xi_{-1}-\xi_{-2})\varphi_{1}(n)+n(2-n)\xi_{-1}\varphi_{1}(n-1)\right].\end{array}$
(50)

These expressions allow to treat the degeneracy of the $\mathbf{T}_{c}$
eigenvalues [48], that was not possible to do in Ref. [44].

Figure 3: Transmittance for normally incident $p$-polarized light on the
finite periodic optical system with unit cell described in Figure 1 and
consisting of $N=16$ cells. The results obtained with both methods, $A$ and
$B$, are shown. They are indistinguishable.

In Figure 3, we compare both methods, $A$ and $B$, i.e. the iteration of Eq.
(30) up to $n=N-1$, from its seminal values and the one that uses the closed
forms of Eqs. (36)-(37), and (49)-(50), to calculate the transmittance for
normally incident $p$-polarized light on the finite periodic optical system
with unit cell described in Figure 1 and consisting of $N=16$ cells. We can
observe that both results are indistinguishable. Notice that the range of
frequencies where the transmittance is attenuated, corresponds with the gaps
for the infinite system (cyan strips and shown in Figure 1), where both
wavenumbers, $k_{1}$ and $k_{2}$, are complex. The first gap (pink strip) is
different in nature and indicates the interval of frequency where both,
$k_{1}$ and $k_{2}$, are pure imaginary, corresponding to exponentially
decaying or increasing optical field oscillations. Gray strips indicate the
interval of frequencies where one of the wavenumbers is imaginary pure and the
other is real. The rest of the intervals of frequency correspond to real
wavenumbers, i.e. propagating waves. The first obvious observation is that
when we have propagating waves, they can interfer constructively or
destructively, in our case partially destructive, as it can be observed in the
oscillations of the transmittance just below the level of maximum
transmission. Also observe the relatively wide intervals of frequency
indicated by gray areas around a gap, indicative of a relatively large
imaginary pure wavenumber, where the transmittance drops drastically for
frequencies below the gap compared with the transmittance for frequencies
above the gap. This effect might indicative of a destructive/constructive
interference between the different exponentially decaying or increasing
optical field oscillations, and can be studied only in the regime of a finite
periodic optical system. This is one of the main results of this work. In
contrast, this effect fades in the relatively narrow intervals of frequency
indicated by gray areas around a gap, possibly due to the interference
process, just discussed, become negligible.

Figure 4: (Color online) Transmittance of the same system of Figure 1 as
function of the number of cells and for normally incident $s$-polarized light.
Figure 5: Transmittance of the same system of Figure 1, as function of the
number of cells for normally incident $p$-polarized light.

In Figure 4, we show the transmittance for the same periodic system as in
Figure 1, as a function of the number of cells for normally incident
$s$-polarized light. In Figure 5, we show the corresponding result for
normally incident $p$-polarized light. Although the band formation is evident,
notice however that the drop in the transmittance around the emerging gaps, as
a function of the number of cells, has larger velocity for incident
$p$-polarization as the number of cells increases. To analize this subtle
effect, it is important to have a reliable method to precisely calculate the
spectra. For instance, notice that around $\omega^{\prime}=1.6$ and for
$N=32$, the transmittance drops significatively for normally incident
$p$-polarized light while for $s$-polarized light transmittance is still
relatively high. This might be due to the above mentioned effect of
destructive/constructive interference between the different exponentially
decaying or increasing field oscillations, and can be thoroughly investigated
with our formalism.

Figure 6: (Color online) Transmittance for either normally incident $s$ or
$p$-polarized light on the finite periodic optical system with unit cell
described in Figure 2 and consisting of $N=32$ cells.

Finally, in Figure 6, we show transmittance for the same periodic system as in
Figure 2 and consisting of $N=32$ cells, for either normally incident $s$ or
$p$-polarized light. Notice that only for normal incidence and isotropic
media, makes no sense to distinguish between polarizations of light. Therefore
both polarizations yield the same result, as is indeed the case. We can
observe that in the relatively wide gray intervals of frequency, the
transmittance drops significatively but not as low as in the gap (pink area).

##  7 Conclusions

In this work, the transmittance of a finite periodic anisotropic optical
system, is obtained in a closed form. The only restrictions on the $4\times 4$
differential propagation matrix of each layer composing the binary unit cell,
$\mathbb{\Delta}$, are $a)$ that its eigenvalues are of the form
$\pm~{}p_{1}$, $\pm~{}p_{2}$, with $p_{1}\neq p_{2}$, and $b)$ its off-
diagonal $2\times 2$-block matrices possess the same symmetric structure in
both layers. Then, for an arbitrary number of cells, $N$, the problem is
reduced to calculating the $N$-th power of the transfer matrix of a unit cell.
The Cayley-Hamilton theorem is invoked and a set of recursive relations are
obtained, that leads to the recursive relation defining the Tetranacci
Polynomials. Such recursive relations are solvable for any invertible $4\times
4$-matrix and provide an algorithm to compute the $N$-power of the transfer
matrix avoiding the multiplication of $N$ matrices. Furthermore, by using the
recently closed expressions for symmetric TP, the transmittance of a finite
periodic stratified and anisotropic optical system, fulfilling the above
mentioned conditions, and for an arbitrary number of unit cells, is obtained
in a closed form. The results for two birefringent uniaxial media with their
optical axis perpendicular to the $z$-axis and for two isotropic media subject
to an external magnetic field oriented along the $z$-axis, exhibiting Faraday
effect, are compared for both methods, finding them indistinguishable.
Finally, lossy media with magnetic anisotropy and/or optical activity can be
treated with our formalism.

##  8 Acknowledgements

ADA and JA were supported by CONACyT under project number A1-S-23120. JCTG
acknowledge financial support from CONACyT, Mexico via posdoctoral grant from
the program Estancias Posdoctorales por México.

##  9 Declaration of interests

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  10 Data availability

No datasets were analyzed or generated during the course of the current study.

##  11 Appendix A: The eigenvector matrix

Let us decompose the $\mathbb{\Delta}$-matrix into $2\times 2$ matrices as
follows,

$\mathbb{\Delta}=\left(\begin{array}[]{cc}\mathbb{\Delta}_{11}&\mathbb{\Delta}_{12}\\\
\mathbb{\Delta}_{21}&\mathbb{\Delta}_{22}\end{array}\right).$ (51)

The conditions for an even characteristic polynomial, Tr $\mathbb{\Delta}=0$
and Tr $\mathbb{\Delta}^{3}=0$ [46, 47, 49], can be fulfilled by the following
sufficient conditions,

$\mbox{Tr }\mathbb{\Delta}_{11}=\mbox{Tr }\mathbb{\Delta}_{22}=0,$ (52)

and

$\mbox{Tr
}\left(\mathbb{\Delta}_{11}\mathbb{\Delta}_{12}\mathbb{\Delta}_{21}+\mathbb{\Delta}_{22}\mathbb{\Delta}_{21}\mathbb{\Delta}_{12}\right)=0,$
(53)

respectively, where we have used the first condition in the second. The
corresponding non-degeneracy condition, $p_{1}\neq~{}p_{2}$, derives in two
additional conditions, $|\mathbb{\Delta}|\neq 0$, and $(\mbox{Tr
}\mathbb{\Delta}^{2})^{2}\neq 16\mbox{ det~{}}\mathbb{\Delta}$. In terms of
$2\times 2$-matrices, the latter reads as

$\left(\mbox{Tr~{}}\mathbb{\Delta}_{11}^{2}+\mbox{Tr~{}}\mathbb{\Delta}_{22}^{2}+2\mbox{
Tr~{}}\mathbb{\Delta}_{12}\mathbb{\Delta}_{21}\right)^{2}\neq 16\mbox{
det~{}}\mathbb{\Delta}.$ (54)

An immediate option to fulfill Eq. (53) is that both matrices,
$\mathbb{\Delta}_{11}$ and $\mathbb{\Delta}_{22}$, besides being traceless
according to the condition of Eq. (52), are null matrices. We discard this
option, since these matrices are related to the diagonal of electrical
permitivity and magnetic permeability tensors, which is very unlikely they are
null [37, 54]. A second option is the case where one or both matrices,
$\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$, are null. Notice that in
this case we may find the widely studied isotropic case, in which the
electromagnetic polarizations decouple and we will not consider here. Then, if
one or both matrices, $\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$, are
null and the medium is not isotropic, the non-degeneracy condition of Eq. (54)
may not be fulfilled. In this case,

$\mbox{det }\mathbb{\Delta}=\frac{1}{4}\mbox{
Tr~{}}\mathbb{\Delta}_{11}^{2}\mbox{ Tr~{}}\mathbb{\Delta}_{22}^{2},$ (55)

and the condition of Eq. (54) becomes,

$\left(\mbox{Tr~{}}\mathbb{\Delta}_{11}^{2}+\mbox{Tr~{}}\mathbb{\Delta}_{22}^{2}\right)^{2}\neq
4\mbox{ Tr~{}}\mathbb{\Delta}_{11}^{2}\mbox{ Tr~{}}\mathbb{\Delta}_{22}^{2},$
(56)

which evidently is not fulfilled if $\mbox{
Tr~{}}\mathbb{\Delta}_{11}^{2}=\mbox{ Tr~{}}\mathbb{\Delta}_{22}^{2}$.
Therefore, if one or both matrices, $\mathbb{\Delta}_{12}$ and
$\mathbb{\Delta}_{21}$, are null, then we must require in addition,

$\mbox{ Tr~{}}\mathbb{\Delta}_{11}^{2}\neq\mbox{
Tr~{}}\mathbb{\Delta}_{22}^{2},$ (57)

to guarantee the non-degeneracy condition. We can include this option, along
with the extra condition of Eq. (57), as a special case of a third option,
where both matrices, $\mathbb{\Delta}_{12}\mathbb{\Delta}_{21}$ and
$\mathbb{\Delta}_{21}\mathbb{\Delta}_{12}$, have zero off-diagonal elements,
since their multiplication with a traceless matrix (Eq. (52)) yield a
traceless matrix. This option in turn implies that both matrices,
$\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$, either have zero diagonal
elements or both have zero off-diagonal elements. Then, we distinguish these
two types of matrices as I and II, respectively. According to Eq. (52), then
$\mathbb{\Delta}_{11}$ and $\mathbb{\Delta}_{22}$, are of type I. As another
example, the matrices $\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$,
appearing in free oscillations of Timoshenko-Ehrenfest beams, are of type II
[49]. If one of the matrices, $\mathbb{\Delta}_{12}$ and
$\mathbb{\Delta}_{21}$, is null, then it can be considered of the same type as
the other non-null matrix. If both matrices, $\mathbb{\Delta}_{12}$ and
$\mathbb{\Delta}_{21}$, are null, then they can be considered of any type,
which is convenient to match the type of the other layer in the binary cell.
In the latter two special cases, the medium has to fulfill in addition the
non-degeneracy condition of Eq. (57).

Let us now see the implications on the eigenvector matrix
$\mathbb{V}=(\vec{V}_{1},\vec{V}_{2},\vec{V}_{3},\vec{V}_{4})$. For each
eigenvalue, $p_{j}$ (or $k_{j}$), $j=1,2,3,4$, the corresponding eigenvector
$\vec{V}_{j}$ is obtained from Eq. (22) and therefore, the matrix
$\mathbb{\Delta}$ satisfies

$\mathbb{\Delta}=\mathbb{V}\mathbb{Y}\mathbb{V}^{-1},$ (58)

where $\mathbb{Y}=\mbox{ diag }(p_{1},p_{2},-p_{1},-p_{2})$,
$p_{1}\neq~{}p_{2}$. We emphasize Eq. (58) is valid only if
$p_{1}\neq~{}p_{2}$, otherwise $\mathbb{\Delta}$ might not have four
independent eigenvectors. However, as mentioned in Section 4, from a practical
point of view, we can recover the case of an isotropic layer, in which
$p_{1}=p_{2}$, to any degree of approximation by taking the anisotropy
arbitrarily small, such that the above expression is still valid to avoid
discontinuous solutions [42, 61, 62]. Let us rewrite Eq. (58) in terms of
$2\times 2$ matrices. First, we define

$\mathbb{V}=\left(\begin{array}[]{cc}\mathbb{V}_{11}&\mathbb{V}_{12}\\\
\mathbb{V}_{21}&\mathbb{V}_{22}\end{array}\right),$ (59)

and

$\mathbb{Y}=\left(\begin{array}[]{cc}\mathbb{\upsilon}&\mathbb{0}\\\
\mathbb{0}&-\mathbb{\upsilon}\end{array}\right),$ (60)

where $\mathbb{0}$ stands for the null $2\times 2$-matrix and
$\mathbb{\upsilon}$ is the diagonal $2\times 2$-matrix,
$\mathbb{\upsilon}=\mbox{ diag }(p_{1},p_{2})$. Also, $\mathbb{V}_{\mu\nu}$,
$\mu,\nu=1,2$, are $2\times 2$-matrices. The substitution into Eq. (58) yields
four $2\times 2$-matrix equations,

$\displaystyle\mathbb{\Delta}_{11}\mathbb{V}_{11}+\mathbb{\Delta}_{12}\mathbb{V}_{21}=\mathbb{V}_{11}\mathbb{\upsilon},$
(61)
$\displaystyle\mathbb{\Delta}_{21}\mathbb{V}_{11}+\mathbb{\Delta}_{22}\mathbb{V}_{21}=\mathbb{V}_{21}\mathbb{\upsilon},$
(62)
$\displaystyle\mathbb{\Delta}_{11}\mathbb{V}_{12}+\mathbb{\Delta}_{12}\mathbb{V}_{22}=-\mathbb{V}_{12}\mathbb{\upsilon},$
(63)
$\displaystyle\mathbb{\Delta}_{21}\mathbb{V}_{12}+\mathbb{\Delta}_{22}\mathbb{V}_{22}=-\mathbb{V}_{22}\mathbb{\upsilon}.$
(64)

Notice that Eqs. (61) and (63) are equivalent if a relation exists between
$\mathbb{V}_{11}$ and $\mathbb{V}_{12}$, and between $\mathbb{V}_{22}$ and
$\mathbb{V}_{21}$, which is an expected result, in the sense that if we know
$\vec{V}_{1}$ and $\vec{V}_{2}$ eigenvectors, then we just need to replace
$p_{1}$ by $-p_{1}$ and $p_{2}$ by $-p_{2}$ to obtain the eigenvectors
$\vec{V}_{3}$ and $\vec{V}_{4}$. A similar statement applies to Eqs. (62) and
(64).

For matrices $\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$, of type I,
i.e. with zero diagonal elements (as $\mathbb{\Delta}_{11}$ and
$\mathbb{\Delta}_{22}$), we have

$\mathbb{\sigma}_{z}\mathbb{\Delta}_{\mu\nu}\mathbb{\sigma}_{z}=-\mathbb{\Delta}_{\mu\nu},\quad\mu,\nu=1,2,$
(65)

where $\mathbb{\sigma}_{z}$ is the Pauli matrix,

$\mathbb{\sigma}_{z}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right).$
(66)

Therefore, Eq. (61) can be rewritten as

$-\mathbb{\sigma}_{z}\mathbb{\Delta}_{11}\mathbb{\sigma}_{z}\mathbb{V}_{11}-\mathbb{\sigma}_{z}\mathbb{\Delta}_{12}\mathbb{\sigma}_{z}\mathbb{V}_{21}=\mathbb{\sigma}_{z}^{2}\mathbb{V}_{11}\mathbb{\upsilon},$
(67)

where we have used $\sigma_{z}^{2}=\mathbf{I}$, and $\mathbf{I}$ is the
identity $2\times 2$-matrix. Factorizing a $\mathbb{\sigma}_{z}$-matrix by the
left we obtain,

$-\mathbb{\Delta}_{11}\mathbb{\sigma}_{z}\mathbb{V_{11}}-\mathbb{\Delta}_{12}\mathbb{\sigma}_{z}\mathbb{V_{21}}=\mathbb{\sigma}_{z}\mathbb{V_{11}}\mathbb{\upsilon},$
(68)

and comparing with Eq. (63) yields

$\displaystyle\mathbb{V}_{12}=a\mathbb{\sigma}_{z}\mathbb{V}_{11}$ (69)
$\displaystyle\mathbb{V}_{22}=a\mathbb{\sigma}_{z}\mathbb{V}_{21},$ (70)

where $a$ is non-zero constant. By substituting Eqs. (69)-(70) into Eq. (64)
and multiplying by the left by $\mathbb{\sigma}_{z}$ we obtain

$a\mathbb{\sigma}_{z}\mathbb{\Delta}_{21}\mathbb{\sigma}_{z}\mathbb{V}_{11}+a\mathbb{\sigma}_{z}\mathbb{\Delta}_{22}\mathbb{\sigma}_{z}\mathbb{V}_{21}=-a\mathbb{\sigma}_{z}^{2}\mathbb{V}_{21}\mathbb{\upsilon},$
(71)

by using Eq. (65) again, this Eq. simplifies to

$-a\mathbb{\Delta}_{21}\mathbb{V}_{11}-a\mathbb{\Delta}_{22}\mathbb{V}_{21}=-a\mathbb{V}_{21}\mathbb{\upsilon},$
(72)

finally, by factorizing $-a$ we recover Eq. (62).

For matrices $\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$, of type II,
i.e. with zero off-diagonal elements,

$\mathbb{\sigma}_{z}\mathbb{\Delta}_{\mu\nu}\mathbb{\sigma}_{z}=\left\\{\begin{array}[]{c}-\mathbb{\Delta}_{\mu\nu},\quad\mbox{
for }\mu=\nu,\\\ \mathbb{\Delta}_{\mu\nu},\quad\mbox{ for
}\mu\neq\nu,\end{array}\right.$ (73)

and by analyzing Eqs. (61)-(64) we obtain instead of Eq. (70),

$\mathbb{V_{22}}=-a\mathbb{\sigma}_{z}\mathbb{V_{21}}.$ (74)

Notice that the above results remain unchanged if one or both matrices,
$\mathbb{\Delta}_{12}$ or $\mathbb{\Delta}_{21}$, are null.

The above procedure on $\mathbb{V}$, can be repeated on matrix
$\mathbb{V^{-1}}$,

$\mathbb{V}^{-1}=\left(\begin{array}[]{cc}\mathbb{U}_{11}&\mathbb{U}_{12}\\\
\mathbb{U}_{21}&\mathbb{U}_{22}\end{array}\right),$ (75)

to yield

$\displaystyle\mathbb{U_{21}}=a^{-1}\mathbb{U_{11}}\mathbb{\sigma}_{z}$ (76)
$\displaystyle\mathbb{U_{22}}=\pm a^{-1}\mathbb{U_{12}}\mathbb{\sigma}_{z},$
(77)

where the plus or minus sign is taken for $\mathbb{\Delta}_{12}$,
$\mathbb{\Delta}_{21}$-matrices of the type I or II, respectively. The inverse
of the constant $a$, is required to satisfy the condition that
$\mathbb{V}\mathbb{V}^{-1}=\mathbb{V}^{-1}\mathbb{V}=\mathbf{I}$, where
$\mathbf{I}$, is the identity $4\times 4$-matrix. Therefore, the $\mathbb{V}$
and $\mathbb{V^{-1}}$ matrices, can be expressed in terms only of its two left
and upper $2\times 2$-matrices, respectively,

$\mathbb{V}=\left(\begin{array}[]{cc}\mathbb{V}_{11}&a\mathbb{\sigma}_{z}\mathbb{V}_{11}\\\
\mathbb{V}_{21}&\pm a\mathbb{\sigma}_{z}\mathbb{V}_{21}\end{array}\right),$
(78)
$\mathbb{V}^{-1}=\left(\begin{array}[]{cc}\mathbb{U}_{11}&\mathbb{U}_{12}\\\
a^{-1}\mathbb{U}_{11}\mathbb{\sigma}_{z}&\pm
a^{-1}\mathbb{U}_{12}\mathbb{\sigma}_{z}\end{array}\right).$ (79)

The result of Eq. (78), states that two of the $\mathbb{\Delta}$-eigenvectors
can be written in terms on the other two, which is not surprising as we stated
above, since we simply have to change the sign of $p_{1}$ and $p_{2}$ in the
calculation of the corresponding eigenvectors. Therefore, the constant $a$,
enters in play since eigenvectors are defined upto a constant. Hence, without
lost of generality we can set $a=1$ from now on.

##  12 Appendix B: The symmetric condition $\eta=\tau$

Let us now construct the matrix
$\mathbb{W}=\mathbb{V}_{1}^{-1}\mathbb{V}_{2}$, where the subindexes label the
first or the second anisotropic layer in the unitary cell, by using Eqs. (78)
and (79) with $a=1$, we obtain,

$\mathbb{W}=\left(\begin{array}[]{cc}\mathbb{\Omega}+\mathbb{\Lambda}&\mathbb{\Gamma}+\iota_{2}\mathbb{\Sigma}\\\
\mathbb{\Gamma}+\iota_{1}\mathbb{\Sigma}&\mathbb{\Omega}+\iota_{1}\iota_{2}\mathbb{\Lambda}\end{array}\right)$
(80)

where $\mathbb{\Omega}$, $\mathbb{\Gamma}$, $\mathbb{\Lambda}$ and
$\mathbb{\Sigma}$, are the $2\times 2$-matrices,

$\begin{array}[]{c}\mathbb{\Omega}=\mathbb{U}_{11}^{(1)}\mathbb{V}_{11}^{(2)},\\\
\mathbb{\Lambda}=\mathbb{U}_{12}^{(1)}\mathbb{V}_{21}^{(2)},\\\
\mathbb{\Gamma}=\mathbb{U}_{11}^{(1)}\sigma_{z}\mathbb{V}_{11}^{(2)},\\\
\mathbb{\Sigma}=\mathbb{U}_{12}^{(1)}\sigma_{z}\mathbb{V}_{21}^{(2)},\end{array}$
(81)

and $\iota_{j}$, $j=1,2$, characterize the type of both of the off-diagonal
matrices $\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$ corresponding to
the $j$-th layer. Thus, $\iota_{j}=\pm 1$ for type I and II, respectively. In
Eq. (81), the superscripts label the first or the second anisotropic layer in
the unitary cell.

Notice $\mathbb{W}^{-1}=\mathbb{V}_{2}^{-1}\mathbb{V}_{1}$. Therefore, the
inverse of $\mathbb{W}$, $\mathbb{W}^{-1}$, has also a similar form,

$\mathbb{W}^{-1}=\left(\begin{array}[]{cc}\mathbb{\Omega}^{\prime}+\mathbb{\Lambda}^{\prime}&\mathbb{\Gamma}^{\prime}+\iota_{1}\mathbb{\Sigma}^{\prime}\\\
\mathbb{\Gamma}^{\prime}+\iota_{2}\mathbb{\Sigma}^{\prime}&\mathbb{\Omega}^{\prime}+\iota_{1}\iota_{2}\mathbb{\Lambda}^{\prime}\end{array}\right)$
(82)

where $\mathbb{\Omega}^{\prime}$, $\mathbb{\Gamma}^{\prime}$,
$\mathbb{\Lambda}^{\prime}$ and $\mathbb{\Sigma}^{\prime}$, are the $2\times
2$-matrices,

$\begin{array}[]{c}\mathbb{\Omega}^{\prime}=\mathbb{U}_{11}^{(2)}\mathbb{V}_{11}^{(1)},\\\
\mathbb{\Lambda}^{\prime}=\mathbb{U}_{12}^{(2)}\mathbb{V}_{21}^{(1)},\\\
\mathbb{\Gamma}^{\prime}=\mathbb{U}_{11}^{(2)}\sigma_{z}\mathbb{V}_{11}^{(1)},\\\
\mathbb{\Sigma}^{\prime}=\mathbb{U}_{12}^{(2)}\sigma_{z}\mathbb{V}_{21}^{(1)}.\end{array}$
(83)

Let us consider the trace of the transfer matrix $\mathbf{T}_{c}$

$\mbox{Tr }\mathbf{T}_{c}=\mbox{Tr
}\mathbb{V}_{1}^{-1}\mathbb{V}_{2}\mathbb{\lambda}_{2}(\ell_{2})\mathbb{V}_{2}^{-1}\mathbb{V}_{1}\mathbb{\lambda}_{1}(\ell_{1})=\mbox{
Tr }\mathbb{W}\lambda_{2}(\ell_{2})\mathbb{W}^{-1}\lambda_{1}(\ell_{1}),$ (84)

where $\mathbb{\lambda}_{j}(\ell_{j})$ is the eigenvalue matrix of Eq. (25)
but corresponding to the $j$-th layer in the unit cell. Also, we have made use
of Eqs. (9) and (27), and of the ciclic property of the trace operation. Now,
by rewriting $\mathbb{\lambda}(\ell_{j})$, $j=1,2$

$\lambda(\ell_{j})=\left(\begin{array}[]{cc}\mathbb{K}_{j}&\mathbb{0}\\\
\mathbb{0}&\mathbb{K}_{j}^{-1}\end{array}\right)$ (85)

where $\mathbb{0}$, stands for the null $2\times 2$-matrix, and
$\mathbb{K}_{j}$, is the diagonal $2\times 2$-matrix $\mathbb{K}_{j}=\mbox{
diag }(e^{ik_{1}^{(j)}\ell_{j}},e^{ik_{2}^{(j)}\ell_{j}})$, where
$k_{1}^{(j)}$ and $k_{2}^{(j)}$ are the eigenvalues of Eq. (25) corresponding
to the $j$-th layer in the unit cell. By expanding the matrix multiplication
in Eq. (84) in terms of $2\times 2$-matrices, we obtain

$\begin{array}[]{ll}\mbox{Tr
}\mathbf{T}_{c}=&\mbox{Tr}\left[(\mathbb{\Omega}+\mathbb{\Lambda})\mathbb{K}_{2}(\mathbb{\Omega}^{\prime}+\mathbb{\Lambda}^{\prime})\mathbb{K}_{1}+(\mathbb{\Gamma}+\iota_{2}\mathbb{\Sigma})\mathbb{K}_{2}^{-1}(\mathbb{\Gamma}^{\prime}+\iota_{2}\mathbb{\Sigma}^{\prime})\mathbb{K}_{1}\right.\\\
&\left.+(\mathbb{\Omega}+\iota_{1}\iota_{2}\mathbb{\Lambda})\mathbb{K}_{2}^{-1}(\mathbb{\Omega}^{\prime}+\iota_{1}\iota_{2}\mathbb{\Lambda}^{\prime})\mathbb{K}_{1}^{-1}+(\mathbb{\Gamma}+\iota_{1}\mathbb{\Sigma})\mathbb{K}_{2}(\mathbb{\Gamma}^{\prime}+\iota_{1}\mathbb{\Sigma}^{\prime})\mathbb{K}_{1}^{-1}\right].\end{array}$
(86)

Similarly,

$\mbox{Tr }\mathbf{T}_{c}^{-1}=\mbox{ Tr
}\mathbb{W}\mathbb{\lambda}_{2}^{-1}(\ell_{2})\mathbb{W}^{-1}\mathbb{\lambda}_{1}^{-1}(\ell_{1}),$
(87)

therefore, to obtain Tr $\mathbf{T}_{c}^{-1}$, we just have to replace
$\mathbb{K}_{j}$ by $\mathbb{K}_{j}^{-1}$ and viceversa in Eq. (86),

$\begin{array}[]{ll}\mbox{Tr
}\mathbf{T}_{c}^{-1}=&\mbox{Tr}\left[(\mathbb{\Omega}+\mathbb{\Lambda})\mathbb{K}_{2}^{-1}(\mathbb{\Omega}^{\prime}+\mathbb{\Lambda}^{\prime})\mathbb{K}_{1}^{-1}+(\mathbb{\Gamma}+\iota_{2}\mathbb{\Sigma})\mathbb{K}_{2}(\mathbb{\Gamma}^{\prime}+\iota_{2}\mathbb{\Sigma}^{\prime})\mathbb{K}_{1}^{-1}\right.\\\
&\left.+(\mathbb{\Omega}+\iota_{1}\iota_{2}\mathbb{\Lambda})\mathbb{K}_{2}(\mathbb{\Omega}^{\prime}+\iota_{1}\iota_{2}\mathbb{\Lambda}^{\prime})\mathbb{K}_{1}+(\mathbb{\Gamma}+\iota_{1}\mathbb{\Sigma})\mathbb{K}_{2}^{-1}(\mathbb{\Gamma}^{\prime}+\iota_{1}\mathbb{\Sigma}^{\prime})\mathbb{K}_{1}\right],\end{array}$
(88)

therefore, this operation leaves the rhs of Eq. (86) the same only if
$\iota_{1}$ and $\iota_{2}$ have the same sign, i.e. if the off-diagonal
matrices $\mathbb{\Delta}_{12}$ and $\mathbb{\Delta}_{21}$ in both layers are
of the same type. Only in this case we have Tr $\mathbf{T}_{c}$= Tr
$\mathbf{T}_{c}^{-1}$.

Finally, if the transfer matrix, $\mathbb{T}_{c}$, for a unitary cell composed
by two anisotropic layers satisfies Eq. (10), then its inverse,
$\mathbf{T}_{c}^{-1}$, satisfies

$(\mathbf{T}_{c}^{-1})^{4}-\tau(\mathbf{T}_{c}^{-1})^{3}+\zeta(\mathbf{T}_{c}^{-1})^{2}-\eta\mathbf{T}_{c}^{-1}+\mathbf{I_{4}}=\mathbb{0}.$
(89)

Thus, the condition, $\eta=\tau$, implies Tr $\mathbf{T}_{c}$=Tr
$\mathbf{T}_{c}^{-1}$.

##  13 Appendix C: Examples of types I and II media

The condition, Tr $\mathbb{\Delta}=0$, can be put in explicit form [54] for
the components of the optical matrix, $\mathbf{M}$,

$\begin{array}[]{rcl}0&=&\rho^{\prime}_{21}+(\frac{k_{x}}{\omega}+\rho^{\prime}_{23})F_{1}+\mu_{23}F_{5},\\\
0&=&\rho_{12}+\epsilon_{13}F_{4}+\rho_{13}F_{8},\\\
0&=&\rho^{\prime}_{12}+\rho^{\prime}_{13}F_{2}+\mu_{13}F_{6},\\\
0&=&\rho_{21}+(-\frac{k_{x}}{\omega}+\rho_{23})F_{7}+\epsilon_{23}F_{3},\end{array}$
(90)

where

$\begin{array}[]{rcl}F_{1}&=&(\rho^{\prime}_{31}\rho_{33}-\epsilon_{31}\mu_{33})/G,\\\
F_{2}&=&((\rho^{\prime}_{32}-\frac{k_{x}}{\omega})\rho_{33}-\epsilon_{32}\mu_{33})/G,\\\
F_{3}&=&(\mu_{31}\rho_{33}-\rho_{31}\mu_{33})/G,\\\
F_{4}&=&(\mu_{32}\rho_{33}-(\rho_{32}+\frac{k_{x}}{\omega})\mu_{33})/G,\\\
F_{5}&=&(\rho^{\prime}_{33}\epsilon_{31}-\epsilon_{33}\rho^{\prime}_{31})/G,\\\
F_{6}&=&(\rho^{\prime}_{33}\epsilon_{32}-\epsilon_{33}(\rho^{\prime}_{32}-\frac{k_{x}}{\omega}))/G,\\\
F_{7}&=&(\rho^{\prime}_{33}\rho_{31}-\epsilon_{33}\mu_{31})/G,\\\
F_{8}&=&(\rho^{\prime}_{33}(\rho_{32}+\frac{k_{x}}{\omega})-\epsilon_{33}\mu_{32})/G,\\\
G&=&\epsilon_{33}\mu_{33}-\rho^{\prime}_{33}\rho_{33},\end{array}$ (91)

therefore, we can discard as candidates for our formalism, media with non
isotropic optical activity, as well as with dielectric or permeability tensors
with non-zero elements in their third column/row off the diagonal. Despite
this restriction, there is still many systems suitable for this formalism.
Here we provide some examples.

The first example consists of a birefringent uniaxial medium with its optical
axis perpendicular to the axis of propagation and rotated from the plane of
incidence, in this case, $\mathbf{\rho}=\mathbf{\rho}^{\prime}=\mathbb{0}$ and
$\mathbf{\mu}=\mathbf{\mu}_{v}\mathbf{I}$, where $\mathbf{I}$ is the identity
$3\times 3$-matrix, and the form of the $\mathbb{\Delta}$ matrix is [37, 55]

$\mathbb{\Delta}=\left(\begin{array}[]{cccc}0&\mu_{v}-\frac{k_{x}^{2}}{\omega^{2}\epsilon_{o}}&0&0\\\
\epsilon_{xx}&0&\epsilon_{xy}&0\\\ 0&0&0&\mu_{v}\\\
\epsilon_{yx}&0&\epsilon_{yy}-\frac{k_{x}^{2}}{\omega^{2}\mu_{v}}&0\end{array}\right),$
(92)

which is of the type I, where

$\begin{array}[]{cc}\epsilon_{xx}&=\epsilon_{e}\cos^{2}\psi+\epsilon_{o}\sin^{2}\psi,\\\
\epsilon_{yy}&=\epsilon_{o}\cos^{2}\psi+\epsilon_{e}\sin^{2}\psi,\\\
\epsilon_{xy}&=\epsilon_{yx}=\frac{1}{2}(\epsilon_{o}-\epsilon_{e})\sin{2\psi},\end{array}$
(93)

$\epsilon_{v}$ and $\mu_{v}$ are the vacuum permitivity and permeability,
respectively, $n_{o}=\sqrt{\epsilon_{o}/\epsilon_{v}}$, and
$n_{e}=\sqrt{\epsilon_{e}/\epsilon_{v}}$, are the ordinary and extraordinary
refractive indexes, respectively, and $\psi$ is the angle of rotation of the
optical axis from the $x$-axis. The corresponding eigenvalues are [55, 69],

$k_{\ell}=\pm\frac{\omega J}{c},\pm\frac{\omega I}{c},\quad\ell=1,2,3,4,$ (94)

where,

$J=\sqrt{n_{o}^{2}-n_{0}^{2}\sin^{2}\theta_{0}},$ (95)

and

$I=\sqrt{n_{e}^{2}-n_{0}^{2}\sin^{2}\theta_{0}\epsilon_{xx}/\epsilon_{o}},$
(96)

with eigenvector matrix
$\mathbb{V}=(\vec{V}_{1},\vec{V}_{2},\vec{V}_{3},\vec{V}_{4})$ [55],

$\mathbb{V}=\left(\begin{array}[]{cccc}\sin\psi&\cos\psi&\sin\psi&\cos\psi\\\
\frac{n_{o}^{2}\sin\psi}{JZ_{v}}&\frac{n_{o}^{2}I\cos\psi}{J^{2}Z_{v}}&-\frac{n_{o}^{2}\sin\psi}{JZ_{v}}&-\frac{n_{o}^{2}I\cos\psi}{J^{2}Z_{v}}\\\
\cos\psi&-\frac{n_{o}^{2}\sin\psi}{J^{2}}&\cos\psi&-\frac{n_{o}^{2}\sin\psi}{J^{2}}\\\
\frac{J\cos\psi}{Z_{v}}&-\frac{n_{o}^{2}I\sin\psi}{J^{2}Z_{v}}&-\frac{J\cos\psi}{Z_{v}}&\frac{n_{o}^{2}I\sin\psi}{J^{2}Z_{v}}\end{array}\right),$
(97)

where $Z_{v}$ is the vacuum impedance. The corresponding inverse matrix,
$\mathbb{V}^{-1}$, is

$\mathbb{V}^{-1}=\frac{J^{2}Z_{v}}{2I(n_{o}^{2}\sin^{2}\psi+J^{2}\cos^{2}\psi)}\left(\begin{array}[]{cccc}\frac{n_{o}^{2}I\sin\psi}{J^{2}Z_{v}}&\frac{I\sin\psi}{J}&\frac{I\cos\psi}{Z_{v}}&\frac{I\cos\psi}{J}\\\
\frac{I\cos\psi}{Z_{v}}&\frac{J^{2}\cos\psi}{n_{o}^{2}}&-\frac{I\sin\psi}{Z_{v}}&-\sin\psi\\\
\frac{n_{o}^{2}I\sin\psi}{J^{2}Z_{v}}&-\frac{I\sin\psi}{J}&\frac{I\cos\psi}{Z_{v}}&-\frac{I\cos\psi}{J}\\\
\frac{I\cos\psi}{Z_{v}}&-\frac{J^{2}\cos\psi}{n_{o}^{2}}&-\frac{I\sin\psi}{Z_{v}}&\sin\psi\end{array}\right).$
(98)

Notice both expressions, are in agreement with Eqs. (78) and Eqs. (79). Also
notice when $\psi=0$, the off-diagonal $2\times 2$ block matrices of the
differential propagation matrix, $\mathbb{\Delta}$, become null and the
eigenvectors $\vec{V}_{1}$ and $\vec{V}_{2}$ become orthogonal, therefore the
$\mathbb{V}$ matrix can be written in either form, in this case in the form
corresponding to type I matrices. However, the non-degeneracy condition of Eq.
(57) is still fulfilled as long as $n_{o}\neq~{}n_{e}$.

With little modifications, we have a second example that consists of an
isotropic medium, that in the presence of an external magnetic field along the
$z$-axis, exhibits the Faraday effect, in this case [56],

$\mathbb{\Delta}=\left(\begin{array}[]{cccc}0&\mu_{v}-\frac{k_{x}^{2}}{\omega^{2}\epsilon}&0&0\\\
\epsilon&0&i\gamma&0\\\ 0&0&0&\mu_{v}\\\
-i\gamma&0&\epsilon-\frac{k_{x}^{2}}{\omega^{2}\mu_{v}}&0\end{array}\right),$
(99)

which has the same structure than the birefringent uniaxial medium and
therefore is of the type I. The corresponding eigenvalues are

$k_{\ell}=\frac{\omega p_{\ell}}{c},\quad\ell=1,2,3,4,$ (100)

where

$p_{\ell}=\sqrt{J\left(J\pm\frac{\gamma}{\sqrt{\epsilon\epsilon_{v}}}\right)},\quad\ell=1,2,$
(101) $J=\sqrt{n^{2}-n_{0}^{2}\sin^{2}\theta_{0}},$ (102)

$n=\sqrt{\epsilon/\epsilon_{v}}$ is the isotropic refraction index, and
$p_{\ell+2}=-p_{\ell}$. The eigenvector matrix,
$\mathbb{V}=(\vec{V}_{1},\vec{V}_{2},\vec{V}_{3},\vec{V}_{4})$, is

$\mathbb{V}=\left(\begin{array}[]{cccc}1&1&1&1\\\
\frac{n^{2}p_{1}}{J^{2}Z_{v}}&\frac{n^{2}p_{2}}{J^{2}Z_{v}}&-\frac{n^{2}p_{1}}{J^{2}Z_{v}}&-\frac{n^{2}p_{2}}{J^{2}Z_{v}}\\\
-in/J&in/J&-in/J&in/J\\\
-i\frac{np_{1}}{JZ_{v}}&i\frac{np_{2}}{JZ_{v}}&i\frac{np_{1}}{JZ_{v}}&-i\frac{np_{2}}{JZ_{v}},\end{array}\right).$
(103)

The corresponding inverse matrix, $\mathbb{V}^{-1}$, is

$\mathbb{V}^{-1}=\frac{1}{4}\left(\begin{array}[]{cccc}1&\frac{J^{2}Z_{v}}{n^{2}p_{1}}&iJ/n&i\frac{JZ_{v}}{np_{1}}\\\
1&\frac{J^{2}Z_{v}}{n^{2}p_{2}}&-iJ/n&-i\frac{JZ_{v}}{np_{2}}\\\
1&-\frac{J^{2}Z_{v}}{n^{2}p_{1}}&iJ/n&-i\frac{JZ_{v}}{np_{1}}\\\
1&-\frac{J^{2}Z_{v}}{n^{2}p_{2}}&-iJ/n&i\frac{JZ_{v}}{np_{2}},\end{array}\right).$
(104)

Notice again both expressions, are in agreement with Eqs. (78) and Eqs. (79),
which in this particular case guarantees the breaking of the time reversal
symmetry [70].

A third example consists of an anisotropic magnetic medium, with permeability
tensor

$\mathbf{\mu}=\left(\begin{array}[]{ccc}\mu_{11}&\mu_{12}&0\\\
\mu_{21}&\mu_{22}&0\\\ 0&0&\mu_{33}\end{array}\right),$ (105)

$\mathbf{\epsilon}=\tilde{\epsilon}\mathbf{I}$, where $\tilde{\epsilon}$ is a
scalar and $\mathbf{\rho}=\mathbf{\rho}^{\prime}=\mathbb{0}$, in this case,

$\mathbb{\Delta}=\left(\begin{array}[]{cccc}0&\mu_{22}-\frac{k_{x}^{2}}{\omega^{2}\tilde{\epsilon}}&0&-\mu_{21}\\\
\tilde{\epsilon}&0&0&0\\\ 0&-\mu_{12}&0&\mu_{11}\\\
0&0&\tilde{\epsilon}-\frac{k_{x}^{2}}{\omega^{2}\mu_{33}}&0\end{array}\right),$
(106)

which is also of the type I. We omit its eigenvalues and eigenvectors.

Finally, our last example consists of a medium with isotropic dielectric,
permeability and optical tensors, i.e.
$\mathbf{\epsilon}=\tilde{\epsilon}\mathbf{I}$, where $\tilde{\epsilon}$ is a
scalar, $\mathbf{\mu}=\mu_{v}\mathbf{I}$,
$\mathbf{\rho}^{\prime}=-i\tilde{\alpha}\mathbf{I}$, and
$\mathbf{\rho}=i\tilde{\alpha}\mathbf{I}$, where $\tilde{\alpha}$ is a scalar
[58, 59, 60]. In this case,

$\mathbb{\Delta}=\left(\begin{array}[]{cccc}0&\mu_{v}\left(1-\frac{k_{x}^{2}}{\omega^{2}(\tilde{\epsilon}\mu_{v}-\tilde{\alpha}^{2})}\right)&-i\tilde{\alpha}\left(1+\frac{k_{x}^{2}}{\omega^{2}(\tilde{\epsilon}\mu_{v}-\tilde{\alpha}^{2})}\right)&0\\\
\tilde{\epsilon}&0&0&-i\tilde{\alpha}\\\ i\tilde{\alpha}&0&0&\mu_{v}\\\
0&i\tilde{\alpha}\left(1+\frac{k_{x}^{2}}{\omega^{2}(\tilde{\epsilon}\mu_{v}-\tilde{\alpha}^{2})}\right)&\tilde{\epsilon}\left(1-\frac{k_{x}^{2}}{\omega^{2}(\tilde{\epsilon}\mu_{v}-\tilde{\rho}\tilde{\rho}^{\prime})}\right)&0\end{array}\right),$
(107)

which is of the type II. The corresponding eigenvalues are

$k_{\ell}=\frac{\omega p_{\ell}}{c},\quad\ell=1,2,3,4,$ (108)

where

$p_{\ell}=\sqrt{\left(\tilde{\alpha}c\pm
n\right)^{2}-n_{0}^{2}\sin^{2}\theta_{0}},\quad\ell=1,2,$ (109)

$n=\sqrt{\epsilon/\epsilon_{v}}$ is the isotropic refraction index, and
$p_{\ell+2}=-p_{\ell}$, with eigenvector matrix,
$\mathbb{V}=(\vec{V}_{1},\vec{V}_{2},\vec{V}_{3},\vec{V}_{4})$,

$\mathbb{V}=\left(\begin{array}[]{cccc}1&1&1&1\\\
\frac{n(n+\tilde{\alpha}c)}{Z_{v}p_{1}}&\frac{n(n-\tilde{\alpha}c)}{Z_{v}p_{2}}&-\frac{n(n+\tilde{\alpha}c)}{Z_{v}p_{1}}&-\frac{n(n-\tilde{\alpha}c)}{Z_{v}p_{2}}\\\
i\frac{(n+\tilde{\alpha}c)}{p_{1}}&i\frac{(\tilde{\alpha}c-n)}{p_{2}}&-i\frac{(n+\tilde{\alpha}c)}{p_{1}}&-i\frac{(\tilde{\alpha}c-n)}{p_{2}}\\\
in/Z_{v}&-in/Z_{v}&in/Z_{v}&-in/Z_{v}\end{array}\right).$ (110)

The corresponding inverse matrix, $\mathbb{V}^{-1}$, is

$\mathbb{V}^{-1}=\frac{1}{4}\left(\begin{array}[]{cccc}1&\frac{Z_{v}p_{1}}{n(n+\tilde{\alpha}c)}&-i\frac{p_{1}}{n+\tilde{\alpha}c}&-iZ_{v}/n\\\
1&\frac{Z_{v}p_{2}}{n(n-\tilde{\alpha}c)}&-i\frac{p_{2}}{\tilde{\alpha}c-n}&iZ_{v}/n\\\
1&-\frac{Z_{v}p_{1}}{n(n+\tilde{\alpha}c)}&i\frac{p_{1}}{n+\tilde{\alpha}c}&-iZ_{v}/n\\\
1&-\frac{Z_{v}p_{2}}{n(n-\tilde{\alpha}c)}&i\frac{p_{2}}{\tilde{\alpha}c-n}&iZ_{v}/n\end{array}\right).$
(111)

Notice once again, that both expressions, are in agreement with Eqs. (78) and
Eqs. (79), which in this particular case and unlike the Faraday rotator,
guarantees the time reversal symmetry [37, 57, 70].

## References

  * [1] Kroning R de L and Penney W G 1931 Proc. R. Soc. London Ser. A 130 499
  * [2] Ashcroft N W and Mermin N D 1976 Solid State Physics (Saunders) Chap. 8.
  * [3] Abelès F 1950 Ann. Phys. (Paris) 5 706
  * [4] Cvetič M and Pičman L 1981 J. Phys. A 14 379
  * [5] Erdös P and Herndon R C 1982 Adv. Phys. 31 65
  * [6] Kiang D 1974 Am. J. Phys. 42 785
  * [7] Griffiths D J and Steinke C A, Am. J. Phys. 69 2
  * [8] Yariv A and Yeh P 1984 Optical Waves in Crystals: Propagation and Control of Laser Radiation (Wiley)
  * [9] Yeh P, Yariv A and Hong C-S 1977 J. Opt. Soc. Am. 67 423
  * [10] Bandelow U and Leonhardt U 1993 Opt. Commun. 101 92
  * [11] Lekner J 1994 J. Opt. Soc. Am. A 11 2892
  * [12] Bendickson J M, Dowling J P and Scalora M 1996 Phys. Rev. E 53 4107
  * [13] Kazanskiy V B and Podlozny V V 1999 Microwave Opt. Technol. Lett. 21 299
  * [14] Abromowitz M and Stegun I A 1965 Handbook of Mathematical Functions (Dover) 777, 782
  * [15] Hiller J, Mendelsohn J D and Rubner M F 2002 Nat. Mater. 1 59
  * [16] Xi J-Q, Schubert M F, Kim J, Schubert E F, Chen M, Lin S-Y, Lie W and Smart J A 2007 Nat. Photonics 1 176
  * [17] Boher P and Houdy P 1992 Ann. Phys. (Paris) 17 141
  * [18] Spiller E 1994, Soft X-Ray Optics (SPIE) Chap. 7
  * [19] Nomura K 2003 Science 300 1269
  * [20] Street R A 2009 Adv. Mater. 21 2007
  * [21] Chopra K L, Paulson P D and Dutta V 2004 Prog. Photovoltaics 12 69
  * [22] Shin B, Gunawan O, Zhu Y, Bojarczuk N A, Chey S J and Guha S 2013 Prog. Photovoltaics 21 72
  * [23] Joannopoulos J D, Meade R D and Winn J N 1995 Photonic Crystals: Molding the Flow of Light (Princeton U.P.)
  * [24] Hasan M Z and Kane C L 2010 Rev. Mod. Phys. 82 3045
  * [25] Kane C L and Mele E. J. 2005 Phys. Rev. Lett. 95 146802
  * [26] Lu L, Joannopoulos J D and Soljaić M 2014 Nat. Phot 8 821829
  * [27] Khanikaev A B, Mousavi S H, Tse W-K, Kargarian M, MacDonald A H and Shvets G 2013 Nat. Mat. 12 233239
  * [28] Schomerus H 2013 Opt. Lett. 380 1912
  * [29] Rechtsman M C, Zeuner J M, Plotnik Y, Lumer Y, Podolsky D, Dreisow F, Nolte S, Segev M and Szameit A 2013 Nature 496 196200
  * [30] Weimann S, Kremer M, Plotnik Y, Lumer Y, Nolte S, Makris K G, Segev M, Rechtsman M C and Szameit A 2017 Nat. Mat. 16 433438
  * [31] Fefferman C L, Lee-Thorp J P and Weinstein M I 2014 PNAS 111 87598763
  * [32] Lee-Thorp J P, Vukićević I, Xu X, Yang J, Fefferman C L, Wong C W and Weinstein M I 2016 Phys. Rev. A 93 033822
  * [33] Levy E and Akkermans E 2017 Eur. Phys. J. Spec. Top. 226 1563
  * [34] Xiao M, Zhang Z Q and Chan C T 2014 Phys. Rev. X 4 021017
  * [35] Choi K H, Ling C W, Lee K F, Tsang Y H and Fung K H 2016 Opt. Lett. 41(7) 1644
  * [36] Teitler S and Henvis B W 1970 J. Opt. Soc. Am. 60 830
  * [37] Berreman D W 1972 J. Opt. Soc. Am. 62 502
  * [38] Yeh P 1979 J. Opt. Soc. Am. 69 742
  * [39] Yeh P 1980 Surf. Sc. 96 41
  * [40] Lin-Chung P J and Teitler S 1984 J. Opt. Soc. Am. A 1 703
  * [41] Li Z-M, Sullivan B T and Parsons R R 1988 Appl. Opt. 27 1334
  * [42] Passler N C and Paarmann A 2017 J. Opt. Soc. Am. B 34 2128
  * [43] Nabol S V, Pankin P S, Maksimov D N and Timofeev I V 2022 Phys. Rev. B 106 245403
  * [44] Lin Y K and McDaniel T J 1969 J. Eng. Ind. 91 1133
  * [45] Lang S 1987 Linear Algebra (Springer) 3rd ed. 241.
  * [46] Barut A O, Zeni J R and Laufer A 1994 J. Phys. A: Math. Gen. 27 5239
  * [47] Barut A O, Zeni J R and Laufer A 1994 J. Phys. A: Math. Gen. 27 6799
  * [48] Laufer A 1997 J. Phys. A: Math. Gen. 30 5455
  * [49] Torres-Guzmán J C, Díaz-de-Anda A, Martínez-Argüello A and Arriaga J 2024 SSRN 4717999
  * [50] Soykan Y 2020 Int. J. Adv. Appl. Math. and Mech. 8(1), 15
  * [51] Leumer N, Marganska M, Muralidharan B and Grifoni M 2020 J. Phys.: Condens. Matter 32 445502
  * [52] Leumer N, Grifoni M, Muralidharan B and Marganska M 2021 Phys. Rev B 103 165432
  * [53] Leumer N 2023 J. Phys. A: Math. Theor. 56 435202
  * [54] Azzam R M A and Bashara N M 1977 Ellipsometry and polarized light (Amsterdam:North-Holland) 350
  * [55] Sosnowski T P 1972 Opt. Commun. 4 408
  * [56] Born M 1933 Optik (Springer)
  * [57] Drude P 1965 The theory of optics (Dover) Ch. VI
  * [58] Treyakov S A, Sihvola A H, Sochava A A and Simovski C R 1998 J. Electromagn. Waves Appl. 12 481
  * [59] Oldano C and Becchi M 1999 Pramana-J. Phys. 53 131
  * [60] Wade J, Hilfiker J N, Brandt J R et al. 2020 Nat. Commun. 11 6137
  * [61] Xu W, Wood L T and Golding T D 2000 Phys. Rev. B 61 1740
  * [62] Zhang S and Wyrowski F 2015 Proc. SPIE 9630 96300A
  * [63] Hall B C (2015) Lie groups, Lie algebras, and representations: An elementary introduction, Graduate Texts in Mathematics (Springer) Section 2
  * [64] Higham N J 2008 Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics) pp. xx+425.
  * [65] Broer W and Podgornik R 2023 Phys. Rev. A 108 (1) 012814
  * [66] Webb W A and Parberry E A 1969 Fibonacci Q. 7(5) 457
  * [67] Hoggatt V Jr. and Long C T 1974 Fibonacci Q. 12(2) 113
  * [68] Özvatan M and Pashaev O K 2017 arXiv:1707.09151
  * [69] Lekner J 1991 J. Phys.: Condens. Matter 3 6121
  * [70] Haldane F D M and Raghu S 2008 Phys. Rev. Lett. 100 013904