\left[\left(-
\sum_{\lambda=1}^2\underline{\dot{A}}_\lambda
\vec{e}_{\lambda\vec{k}}\right)\cdot \underline{\vec{P}}^*+ \left(
\imath\vec{k}\times\sum_{\lambda=1}^2\underline{A}_\lambda
\vec{e}_{\lambda\vec{k}}\right)\cdot
\underline{\vec{M}}^*+h.c\right]\nonumber\\
&&\end{aligned}$$ where the polarizations $\underline{\vec{P}}$ and $
\underline{\vec{M}}$ have been defined previously in terms of the dynamical variables of the $E$ and $M$ fields in equations (\[c13\]) and (\[c13.1\]). The Lagrangian (\[c19\]) can now be used to define the canonical conjugate momenta of the system. For any wave vector $\vec{k}$ in half space $k_z\geq 0$ this momenta are defined as $$\begin{aligned}
\label{c23}
-\underline{D}_\lambda (\vec{k},t)&=&\frac{\delta
\underline{L}}{\delta \left(\underline{\dot{A}}^*_\lambda
(\vec{k},t)\right)}=
\varepsilon_0\dot{A}_\lambda (\vec{k},t)-\vec{e}_{\lambda\vec{k}}\cdot\underline{\vec{P}}(\vec{k},t)\nonumber\\
\underline{Q}_{\omega\lambda}(\vec{k},t)&=&\frac{\delta
\underline{L}}{\delta\left(
\underline{\dot{X}}^*_{\omega\lambda}(\vec{k},t)\right)}=
\underline{\dot{X}}_{\omega\lambda}(\vec{k},t)\nonumber\\
\underline{\Pi}_{\omega\lambda}(\vec{k},t)&=&\frac{\delta
\underline{L}}{\delta
\left(\underline{\dot{Y}}^*_{\omega\lambda}(\vec{k},t)\right)}=
\underline{\dot{Y}}_{\omega\lambda}(\vec{k},t)\end{aligned}$$ The total system are quantized canonically in a standard method by imposing equal-time commutation relations between the coordinates of the system and their conjugates variables as follows $$\begin{aligned}
\label{c24}
\left[ \underline{A}^*_\lambda ( \vec{k},t) \ ,\
-\underline{D}_{\lambda'}( \vec{k'} ,
t)\right]&=&\imath\hbar\delta_{\lambda\lambda'}\delta(\vec{k}-\vec{k'})\nonumber\\
\left[ \underline{X}^*_{\omega\lambda }( \vec{k},t) \ ,\
\underline{Q}_{\omega'\lambda'}( \vec{k'} ,
t)\right]&=&\imath\hbar\delta_{\lambda\lambda'}\delta(\omega-\omega')\delta(\vec{k}-\vec{k'})\nonumber\\
\left[ \underline{Y}^*_{\omega\lambda }( \vec{k},t) \ ,\
\underline{\Pi}_{\omega'\lambda'}( \vec{k'} ,
t)\right]&=&\imath\hbar\delta_{\lambda\lambda'}\delta(\omega-\omega')\delta(\vec{k}-\vec{k'})\nonumber\\
&&\end{aligned}$$ Using the Lagrangian (\[c19\]) and the conjugates momenta introduced in (\[c23\]) the Hamiltonian of the total system can be written in the form $$\begin{aligned}
\label{c25}
H(t)&=&\int' d^3k \sum_{\lambda=1}^2\left(
\frac{|\underline{D}_{\lambda}-\vec{e}_{\lambda\vec{k}}\cdot
\underline{\vec{P}}|^2}{\varepsilon_0}+\frac{|\vec{k}\underline{A}_\lambda
|^2}{\mu_0}\right)+\int' d^3k\frac{|\imath\vec{k}\cdot\underline{\vec{P}}|^2}{\varepsilon_0|\vec{k}|^2}\nonumber\\
&-&\int' d^3k \left[\left(
\imath\vec{k}\times\sum_{\lambda=1}^2\underline{A}_{\lambda}\vec{e}_{\lambda\vec{k}}\right)\cdot
\underline{\vec{M}}^*+h.c\right]
\nonumber\\
&+&\int_0^\infty d\omega\int' d^3k\sum_{\lambda=1}^3\left( |\
\underline{Q}_{\omega \lambda}|^2+\omega^2|\ \underline{X}_{\omega
\lambda}|^2 + |\ \underline{\Pi}_{\omega \lambda}|^2+\omega^2|\
\underline{Y}_{\omega\lambda}|^2\right)\nonumber\\
&&\end{aligned}$$

Maxwell equations
-----------------

