of the tensors $ \underline{f}$ and $\underline{g}$, respectively.

 Canonical quantization
=======================

Following the standard approach, we choose the Coulomb gauge $\vec{k}\cdot\underline{\vec{A}}(\vec{k},t)=0$ to quantize electromagnetic field. In this gauge the vector potential $\vec{A}$ is a purely transverse field and can be decomposed along the unit polarization vectors $\vec{e}_{\lambda\vec{k}}\hspace {1.00
cm}\lambda=1,2$ which are orthogonal to each other and to the wave vector $\vec{k}$. $$\label{c16}
\underline{\vec{A}}(\vec{k},t)=\sum_{ \lambda=1}^2\
\underline{A}_\lambda( \vec{k},t)\vec{e}_{\lambda\vec{k}}$$ Although the vector potential is purely transverse, but the dynamical fields $\underline{\vec{X}}_\omega$ and $\underline{\vec{Y}}_\omega$ may have both transverse and longitudinal parts and can be expanded along the three mutually orthogonal unit vectors $\vec{e}_{\lambda\vec{k}}\hspace {0.50 cm}\lambda=1,2$ and $\vec{e}_{3\vec{k}}=\hat{k}=\frac{\vec{k}}{|\vec{k}|}$ as $$\begin{aligned}
\label{c17}
\underline{\vec{X}}_\omega(\vec{k},t)&=&\sum_{ \lambda=1}^3
\underline{X}_{\omega\lambda}( \vec{k},t)\vec{e}_{\lambda\vec{k}}\nonumber\\
\underline{\vec{Y}}_\omega(\vec{k},t)&=&\sum_{ \lambda=1}^3
\underline{Y}_{\omega\lambda}( \vec{k},t)\vec{e}_{\lambda\vec{k}}\end{aligned}$$ Furthermore in Coulomb gauge the Euler- Lagrange equation (\[c12\]) can be used to eliminate the extra degree of freedom $\underline{\varphi}(\vec{k},t)$ from the Larangian of the system $$\label{c18}
\underline{\varphi}(\vec{k},t)=-\frac{\imath\vec{k}\cdot\underline{\vec{P}}(\vec{k},t)}{\varepsilon_0
|\vec{k}|^2 }$$ The total Lagrangian (\[c7\])-(\[c10\]) can now be rewritten in terms of the independent dynamical variables $\underline{A}_\lambda\hspace{0.50 cm}\lambda=1,2$ and $
\underline{\vec{X}}_{\omega\lambda} \ ,\
\underline{\vec{Y}}_{\omega\lambda}  \hspace{0.50
cm}\lambda=1,2,3 $ which constitute completely the coordinates of the total system $$\begin{aligned}
\label{c19}
 \underline{L}(t)&=&\int_0^\infty d\omega\ \int'd^3k
\sum_{\lambda=1}^3\left(|\dot{\underline{\vec{X}}}_{\omega\lambda}|^2-\omega^2|\underline{\vec{X}}_{\omega\lambda}|^2
 +|\dot{\underline{\vec{Y}}}_{\omega\lambda}|^2-\omega^2|\underline{\vec{Y}}_{\omega\lambda}|^2\right)\nonumber\\
& +& \int'd^3k
\sum_{\lambda=1}^2\left(\varepsilon_0|\dot{\underline{{A}}_\lambda}|^2
-\frac{|\vec{k}\ \underline{A}_{\lambda}|^2}{\mu_0} \right)
-\frac{1}{\varepsilon_0} \int' d^3k\
\frac{|\imath\vec{k}\cdot\underline{\vec{P}}|^2}{|\vec{k}|^2}\nonumber\\
& +&\int' d^3k 