# Canonical quantization of electromagnetic field in the presence of nonlinear
anisotropic magnetodielectric medium with spatial-temporal dispersion

M. Amooshahi1<EMAIL_ADDRESS>
1 Faculty of science, University of Isfahan ,Hezar Jarib Ave., Isfahan,Iran

###### Abstract

Modeling a nonlinear anisotropic magnetodielectric medium with spatial-
temporal dispersion by two continuum collections of three dimensional harmonic
oscillators, a fully canonical quantization of the electromagnetic field is
demonstrated in the presence of such a medium. Some coupling tensors of
various ranks are introduced that couple the magnetodielectric medium with the
electromagnetic field. The polarization and magnetization fields of the medium
are defined in terms of the coupling tensors and the oscillators modeling the
medium. The electric and magnetic susceptibility tensors of the medium are
obtained in terms of the coupling tensors. It is shown that the electric field
satisfy an integral equation in frequency domain. The integral equation is
solved by an iteration method and the electric field is found up to an
arbitrary accuracy.

PACS number

12.20.Ds

##  1 Introduction

The quantum features of the electromagnetic field can be influenced by the
magnetodielectric media. For example the Casimir effect [1] \- [10] and the
spontaneous emission rate of an initially excited atom [11]-[14] are modified
in the presence of magnetodielectric media. There is a canonical quantization
approach for the macroscopic electromagnetic field in a linear absorbing
dielectric based on the damped polarization model[15]. In this scheme the
electric polarization field of the medium is appeared in the Lagrangian of the
total system as a part of the degrees of freedom of the medium. The other
parts of the degrees of freedom of the absorbing dielectric are related to the
dynamical variables of a reservoir coupled to the polarization field in order
to inclusion the absorption property of the medium. The Hamiltonian of the
total system is diagonalized in two steps. The first step is diagonalization
of the polarization-reservoir part and the second step is diagonalization of
the total Hamiltonian. This method has been generalized to inhomogeneous
dielectrics [16]and anisotropic magnetodielectric media [17].
There is a new fully canonical quantization of the electromagnetic field in
the presence of linear absorbing anisotropic magnetodielectric media in which
the medium is modeled by two continuum collections of three dimensional space
dependent harmonic oscillators[18],[19]. One of the collection describes the
electric property and the another collection explains the magnetic property of
the magnetodielectric medium. In contrast to the damped polarization model it
is not needed the electric and magnetic polarization fields to be appeared in
the Lagrangian of the total system. The two collections of the harmonic
oscillators solely constitute the degrees of freedom of the medium. In fact
the space dependent harmonic oscillators are able to describe both absorption
and polarization properties of the medium. This means that in this approach
the electric and magnetic polarization fields are defined in terms of the
harmonic oscillators modeling the medium and the coupling tensors which couple
the medium and electromagnetic field. Among the advantages of this approach is
that this approach dose not need the processes of the diagonalization of the
Hamiltonians of the different parts of the total system and is simply
generalizable to the nonlinear anisotropic magnetodielectric medium with
spatial-temporal dispersion. The aim of the present work is generalization of
the mentioned scheme for a nonlinear anisotropic magnetodielectric medium with
spatial-temporal dispersion. The canonical quantization of a three dimensional
system in the presence of a nonlinear anisotropic absorbing environment has
been done previously and the effect of nonlinearity of the environment on the
spontaneous emission of an initially excited two level atom has been
investigated[20]

##  2 The Lagrangian of the total system

Modeling the nonlinear anisotropic magnetodielectric medium by two continuum
collections three dimensional harmonic oscillators, the Lagrangian of the
total system, that is the electromagnetic field plus the magnetodielectric
medium, is introduced as the sum of three parts

$L(t)=L_{em}(t)+L_{res}(t)+L_{int}(t)$ (1)

The first part is the Lagrangian of the electromagnetic field which as usual
can be written as

$L_{em}(t)=\int d^{3}r\left[\frac{1}{2}\varepsilon_{0}{\bf E}^{2}({\bf
r},t)+\frac{{\bf B}^{2}({\bf r},t)}{2\mu_{0}}\right]$ (2)

where ${\bf E}=-\vec{\nabla\varphi}-\frac{\partial{\bf A}}{\partial t}$ and
${\bf B}=\nabla\times{\bf A}$ and ${\bf A},\varphi$ are respectively the
vector and scalar potential which constitute the dynamical variables related
to of the electromagnetic field. The second part in (1) is the Lagrangian of
the magnetodielectric medium. If we denote the two continuum collections of
three dimensional harmonic oscillators modeling the medium by ${\bf
X}_{\omega}({\bf r},t)$ and ${\bf Y}_{\omega}({\bf r},t)$, with continuous
parameter $\omega$, then the Lagrangian of the medium is written as

$\displaystyle L_{res}(t)$ $\displaystyle=$ $\displaystyle\int
d^{3}r\int_{0}^{\infty}d\omega\left[\frac{1}{2}\dot{{\bf X}}_{\omega}^{2}({\bf
r},t)-\frac{1}{2}\omega^{2}{\bf X}^{2}_{\omega}({\bf r},t)\right]$ (3)
$\displaystyle+$ $\displaystyle\int
d^{3}r\int_{0}^{\infty}d\omega\left[\frac{1}{2}\dot{{\bf Y}}^{2}_{\omega}({\bf
r},t)-\frac{1}{2}\omega^{2}{\bf Y}^{2}_{\omega}({\bf r},t)\right]$

It will be shown that the two collections ${\bf X}_{\omega}({\bf r},t)$ and
${\bf Y}_{\omega}({\bf r},t)$ describe, respectively, the electric and
magnetic properties of the medium. That is the electric polarization field of
the medium is defined in terms of ${\bf X}_{\omega}({\bf r},t)$ and the
magnetic polarization field of the medium is expressed in terms of ${\bf
Y}_{\omega}({\bf r},t)$.
The third term in (1) is the interaction part of the electromagnetic field and
the magnetodielectric medium that for a nonlinear anisotropic medium with
spatial-temporal dispersion is proposed as

$\displaystyle L_{int}(t)=\int_{0}^{\infty}d\omega_{1}\int d^{3}r\int
d^{3}r_{1}f^{(1)}_{ii_{1}}(\omega_{1},{\bf r},{\bf r}_{1})E^{i}({\bf
r},t)X^{i_{1}}_{\omega_{1}}({\bf r}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}r\int d^{3}r_{1}\int
d^{3}r_{2}[f^{(2)}_{ii_{1}i_{2}}(\omega_{1},\omega_{2},{\bf r},{\bf
r}_{1},{\bf r}_{2})$ $\displaystyle\times E^{i}({\bf
r},t)X^{i_{1}}_{\omega_{1}}({\bf r}_{1},t)X^{i_{2}}_{\omega_{2}}({\bf
r}_{2},t)]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int_{0}^{\infty}d\omega_{3}\int
d^{3}r\int d^{3}r_{1}\int d^{3}r_{2}\int
d^{3}r_{3}[f^{(3)}_{ii_{1}i_{2}i_{3}}(\omega_{1},\omega_{2},\omega_{3},{\bf
r},{\bf r}_{1},{\bf r}_{2},{\bf r}_{3})$ $\displaystyle\times E^{i}({\bf
r},t)X^{i_{1}}_{\omega_{1}}({\bf r}_{1},t)X^{i_{2}}_{\omega_{2}}({\bf
r}_{2},t)X^{i_{3}}_{\omega_{3}}({\bf r}_{3},t)]+...$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int d^{3}r\int
d^{3}r_{1}g^{(1)}_{ii_{1}}(\omega_{1},{\bf r},{\bf r}_{1})B^{i}({\bf
r},t)Y^{i_{1}}_{\omega_{1}}({\bf r}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}r\int d^{3}r_{1}\int
d^{3}r_{2}[g^{(2)}_{ii_{1}i_{2}}(\omega_{1},\omega_{2},{\bf r},{\bf
r}_{1},{\bf r}_{2})$ $\displaystyle\times B^{i}({\bf
r},t)Y^{i_{1}}_{\omega_{1}}({\bf r}_{1},t)Y^{i_{2}}_{\omega_{2}}({\bf
r}_{2},t)]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int_{0}^{\infty}d\omega_{3}\int
d^{3}r\int d^{3}r_{1}\int d^{3}r_{2}\int
d^{3}r_{3}[g^{(3)}_{ii_{1}i_{2}i_{3}}(\omega_{1},\omega_{2},\omega_{3},{\bf
r},{\bf r}_{1},{\bf r}_{2},{\bf r}_{3})$ $\displaystyle\times B^{i}({\bf
r},t)Y^{i_{1}}_{\omega_{1}}({\bf r}_{1},t)Y^{i_{2}}_{\omega_{2}}({\bf
r}_{2},t)Y^{i_{3}}_{\omega_{3}}({\bf r}_{3},t)]+...$

where the tensors ${\bf f}^{(1)},{\bf f}^{(2)},{\bf f}^{(3)},...$ and ${\bf
g}^{(1)},{\bf g}^{(2)},{\bf g}^{(3)},...$ are the coupling tensors between the
electromagnetic field and the nonlinear magnetodielectric medium. These
coupling tensors play an important role in this quantization approach. It will
be shown that the electric polarization field is obtained in terms of the
coupling tensors ${\bf f}^{(1)},{\bf f}^{(2)},{\bf f}^{(3)},...$ and the
harmonic oscillators ${\bf X}_{\omega}({\bf r},t)$. Also the electric
susceptibilities tensors of the medium are expressed in terms of the coupling
tensors ${\bf f}^{(1)},{\bf f}^{(2)},{\bf f}^{(3)},...$. Therefore these
tensors are an indication for the ability of the polarization of the medium.
Similarly the magnetic polarization field of the magnetodielectric medium is
defined in terms of the coupling tensors ${\bf g}^{(1)},{\bf g}^{(2)},{\bf
g}^{(3)},...$ and the harmonic oscillators ${\bf Y}_{\omega}({\bf r},t)$. As
well the magnetic susceptibility tensors of the medium are written in terms of
the coupling tensors ${\bf g}^{(1)},{\bf g}^{(2)},{\bf g}^{(3)},...$ and
therefore these tensors explain the ability of magnetization of the medium.
The Lagrangian (1)-(2) is a nonlocal one and there is some problems for a
canonical quantization using this Lagrangian. The easiest way to overcome the
difficulties related to a nonlocal Lagrangian is to work in the reciprocal
space. Therefore we write the classical fields and the coupling tensors
appeared in the Lagrangian (1)-(2) in terms of their spatial fourier
transforms. For example for the field ${\bf X}_{\omega}({\bf r},t)$ we write

${\bf X}_{\omega}({\bf r},t)=\frac{1}{(2\pi)^{\frac{3}{2}}}\int
d^{3}k{\bf\underline{X}}_{\omega}({\bf k},t)e^{i{\bf k}\cdot{\bf r}}$ (5)

Similarly the components of the coupling tensors ${\bf f}^{(n)}$ and ${\bf
g}^{(n)}$ can be written in terms of their Fourier transforms as

$\displaystyle
f^{(n)}_{ii_{1}i_{2}...i_{n}}(\omega_{1},\omega_{2},...\omega_{n},{\bf r},{\bf
r}_{1},...,{\bf r}_{n})$ $\displaystyle=\frac{1}{\frac{3(n+1)}{2}}\int
d^{3}k\int d^{3}k_{1}...\int
d^{3}k_{n}\underline{f}^{(n)}_{ii_{1}...i_{n}}(\omega_{1},...,\omega_{n},{\bf
k},{\bf k}_{1},...,{\bf k}_{n})e^{i{\bf k}\cdot{\bf r}-i{\bf k}_{1}\cdot{\bf
r}_{1}-...-i{\bf k}_{n}\cdot{\bf r}_{n}}$ $\displaystyle
g^{(n)}_{ii_{1}i_{2}...i_{n}}(\omega_{1},\omega_{2},...\omega_{n},{\bf r},{\bf
r}_{1},...,{\bf r}_{n})$ $\displaystyle=\frac{1}{\frac{3(n+1)}{2}}\int
d^{3}k\int d^{3}k_{1}...\int
d^{3}k_{n}\underline{g}^{(n)}_{ii_{1}...i_{n}}(\omega_{1},...,\omega_{n},{\bf
k},{\bf k}_{1},...,{\bf k}_{n})e^{i{\bf k}\cdot{\bf r}-i{\bf k}_{1}\cdot{\bf
r}_{1}-...-i{\bf k}_{n}\cdot{\bf r}_{n}}$

Since the fields appeared in Lagrangian (1)-(2)are real valued , we have

$\underline{{\bf X}}_{\omega}({\bf k},t)=\underline{{\bf X}}^{*}_{\omega}({\bf
k},t)$ (7)

for the field ${\bf X}_{\omega}({\bf r},t)$ and the other fields in the
Lagrangian (1)-(2). Similar relations for the real valued coupling tensors
${\bf f}^{(n)}$ and ${\bf g}^{(n)}$ are as follows

$\displaystyle\underline{f}^{(n)}_{ii_{1}...i_{n}}(\omega_{1},...,\omega_{n},{\bf
k},{\bf k}_{1},...,{\bf k}_{n})$ $\displaystyle=$
$\displaystyle(\underline{f}^{(n)})^{*}_{ii_{1}...i_{n}}(\omega_{1},...,\omega_{n},-{\bf
k},-{\bf k}_{1},...,-{\bf k}_{n})$
$\displaystyle\underline{g}^{(n)}_{ii_{1}...i_{n}}(\omega_{1},...,\omega_{n},{\bf
k},{\bf k}_{1},...,{\bf k}_{n})$ $\displaystyle=$
$\displaystyle(\underline{g}^{(n)})^{*}_{ii_{1}...i_{n}}(\omega_{1},...,\omega_{n},-{\bf
k},-{\bf k}_{1},...,-{\bf k}_{n})$

Because of the relations (7) and (2) the independent dynamical variables can
be taken into account by restricting the integration domain to the half space
$k_{z}\geq 0$ in the reciprocal space. Substituting the real valued fields and
the coupling tensors ${\bf f}^{(n)}$ and ${\bf g}^{(n)}$ in the total
Lagrangian (1)-(2) by their Fourier transforms and using the relation (7) for
${\bf X}_{\omega}({\bf r},t)$, and similar relations for the other fields and
also using (2),one can rewrite the Lagrangian of the total system as

$\displaystyle
L(t)=\int^{\prime}d^{3}k\left[\varepsilon_{0}|\underline{\dot{{\bf A}}}({\bf
k},t)|^{2}+\varepsilon_{0}k^{2}|\underline{\varphi}({\bf
k},t)|^{2}-\frac{|{\bf k}\times\underline{{\bf A}}({\bf
k},t)|^{2}}{\mu_{0}}\right]$
$\displaystyle+\varepsilon_{0}\int^{\prime}d^{3}k\left[-i{\bf
k}\cdot\underline{\dot{A}}({\bf k},t)\underline{\varphi}^{*}({\bf
k},t)+H.C\right]$
$\displaystyle+\int_{0}^{\infty}d\omega\int^{\prime}d^{3}k\left[|\underline{\dot{{\bf
X}}}_{\omega}({\bf k},t)|^{2}-\omega^{2}|\underline{{\bf X}}_{\omega}({\bf
k},t)|^{2}\right]$
$\displaystyle+\int_{0}^{\infty}d\omega\int^{\prime}d^{3}k\left[|\underline{\dot{{\bf
Y}}}_{\omega}({\bf k},t)|^{2}-\omega^{2}|\underline{{\bf Y}}_{\omega}({\bf
k},t)|^{2}\right]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\left[\underline{f}^{(1)}_{ij}(\omega_{1},{\bf
k},{\bf k_{1}})\underline{E}^{*^{i}}({\bf
k},t)\underline{X}^{j}_{\omega_{1}}({\bf k_{1}},t)+H.C\right]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\left[\underline{f}^{(1)}_{ij}(\omega_{1},-{\bf
k},{\bf k_{1}})\underline{E}^{i}({\bf k},t)\underline{X}^{j}_{\omega_{1}}({\bf
k_{1}},t)+H.C\right]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{f}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},{\bf k}_{2})$ $\displaystyle\times\underline{E}^{*^{i}}({\bf
k},t)\underline{X}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{X}^{k}_{\omega_{2}}({\bf k}_{2},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{f}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},-{\bf k}_{2})$ $\displaystyle\times\underline{E}^{*^{i}}({\bf
k},t)\underline{X}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{X}^{*^{k}}_{\omega_{2}}({\bf k}_{2},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{f}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},-{\bf k}_{1},{\bf k}_{2})$ $\displaystyle\times\underline{E}^{*^{i}}({\bf
k},t)\underline{X}^{*^{j}}_{\omega_{1}}({\bf
k}_{1},t)\underline{X}^{k}_{\omega_{2}}({\bf k}_{2},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{f}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},-{\bf k}_{1},-{\bf k}_{2})$ $\displaystyle\times\underline{E}^{*^{i}}({\bf
k},t)\underline{X}^{*^{j}}_{\omega_{1}}({\bf
k}_{1},t)\underline{X}^{*^{k}}_{\omega_{2}}({\bf k}_{2},t)+H.C]+....$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}[\underline{g}^{(1)}_{ij}(\omega_{1},{\bf
k},{\bf k}_{1})\underline{B}^{*^{i}}({\bf
k},t)\underline{Y}^{j}_{\omega_{1}}({\bf k}_{1},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}[\underline{g}^{(1)}_{ij}(\omega_{1},-{\bf
k},{\bf k}_{1})\underline{B}^{i}({\bf k},t)\underline{Y}^{j}_{\omega_{1}}({\bf
k}_{1},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{g}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},{\bf k}_{2})$ $\displaystyle\times\underline{B}^{*^{i}}({\bf
k},t)\underline{Y}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{Y}^{k}_{\omega_{2}}({\bf k}_{2},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{g}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},-{\bf k}_{2})$ $\displaystyle\times\underline{B}^{*^{i}}({\bf
k},t)\underline{Y}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{Y}^{*^{k}}_{\omega_{2}}({\bf k}_{2},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{g}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},-{\bf k}_{1},{\bf k}_{2})$ $\displaystyle\times\underline{B}^{*^{i}}({\bf
k},t)\underline{Y}^{*^{j}}_{\omega_{1}}({\bf
k}_{1},t)\underline{Y}^{k}_{\omega_{2}}({\bf k}_{2},t)+H.C]$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int^{\prime}d^{3}k\int^{\prime}d^{3}k_{1}\int^{\prime}d^{3}k_{2}[\underline{g}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},-{\bf k}_{1},-{\bf k}_{2})$ $\displaystyle\times\underline{B}^{*^{i}}({\bf
k},t)\underline{Y}^{*^{j}}_{\omega_{1}}({\bf
k}_{1},t)\underline{Y}^{*^{k}}_{\omega_{2}}({\bf k}_{2},t)+H.C]+...$

where $\int^{\prime}d^{3}k$ denote the integration over the half space
$k_{z}\geq 0$ and $\underline{{\bf E}}({\bf k},t)=-\underline{\dot{{\bf
A}}}({\bf k},t)-i{\bf k}\underline{\varphi}({\bf k},t)$ and $\underline{{\bf
B}}({\bf k},t)=i{\bf k}\times\underline{{\bf A}}({\bf k},t)$.In the Lagrangian
(2) the three points … indicate the other terms which are dependent on the
coupling tensors more than the third rank. Hereafter we apply the symbols
$\int^{\prime}d^{3}k$ and $\int d^{3}k$ for the integration over the half
space$k_{z}\geq 0$ and integration over the total reciprocal space,
respectively.

##  3 Classical Euler-Lagrange equations

The Lagrangian (2) does not involve the space derivatives of the dynamical
variables and one can easily obtain the classical Euler-Lagrange equations for
the Fourier transforms of the fields appeared in the Lagrangian of the total
system. The classical Euler-Lagrange equation for the scalar potential
$\underline{\varphi}({\bf k},t)$ and vector potential $\underline{{\bf
A}}({\bf k},t)$ are respectively as follows

$\displaystyle\frac{d}{dt}(\frac{\delta
L}{\delta\underline{\dot{\varphi}}^{*}({\bf k},t)})-\frac{\delta
L}{\delta\underline{\varphi}^{*}({\bf k},t)}=0$
$\displaystyle\Rightarrow-\varepsilon_{0}i{\bf k}\cdot\underline{\dot{{\bf
A}}}({\bf k},t)+\varepsilon_{0}k^{2}\underline{\varphi}({\bf k},t)=-i{\bf
k}\cdot\underline{{\bf P}}({\bf k},t)$ (10)
$\displaystyle\frac{d}{dt}(\frac{\delta
L}{\delta\underline{\dot{A}}_{i}^{*}({\bf k},t)})-\frac{\delta
L}{\delta\underline{A}_{i}^{*}({\bf k},t)}=0\hskip 56.9055pti=1,2,3$
$\displaystyle\Rightarrow\mu_{0}\varepsilon_{0}\underline{\ddot{{\bf A}}}({\bf
k},t)+\mu_{0}\varepsilon_{0}i{\bf k}\underline{\dot{\varphi}}({\bf k},t)-{\bf
k}\times({\bf k}\times\underline{{\bf A}}({\bf k},t))$
$\displaystyle=\mu_{0}\underline{\dot{{\bf P}}}({\bf k},t)+i\mu_{0}{\bf
k}\times\underline{{\bf M}}({\bf k},t)$

where $\underline{{\bf P}}({\bf k},t)$ and $\underline{{\bf M}}({\bf k},t)$
are the Fourier transforms of the electric and magnetic polarization densities
of the nonlinear magnetodielectric medium, respectively. The components of the
polarization densities are defined in terms of the coupling tensors
$\underline{{\bf f}}^{(n)},\underline{{\bf g}}^{(n)}$ and the the harmonic
oscillators modeling the medium as

$\displaystyle\underline{P}_{i}({\bf k},t)=\int_{0}^{\infty}d\omega_{1}\int
d^{3}k_{1}\underline{f}^{(1)}_{ij}(\omega_{1},{\bf k},{\bf
k}_{1})\underline{X}^{j}_{\omega_{1}}({\bf k}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}k_{1}\int d^{3}k_{2}\underline{f}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},{\bf k}_{2})\underline{X}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{X}^{k}_{\omega_{2}}({\bf k}_{2},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int_{0}^{\infty}d\omega_{3}\int
d^{3}k_{1}\int d^{3}k_{2}\int
d^{3}k_{3}\underline{f}^{(3)}_{ijkl}(\omega_{1},\omega_{2},\omega_{3},{\bf
k},{\bf k}_{1},{\bf k}_{2},{\bf k}_{3})$
$\displaystyle\times\underline{X}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{X}^{k}_{\omega_{2}}({\bf
k}_{2},t)\underline{X}^{l}_{\omega_{3}}({\bf k}_{3},t)+....$
$\displaystyle\underline{M}_{i}({\bf k},t)=\int_{0}^{\infty}d\omega_{1}\int
d^{3}k_{1}\underline{g}^{(1)}_{ij}(\omega_{1},{\bf k},{\bf
k}_{1})\underline{Y}^{j}_{\omega_{1}}({\bf k}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}k_{1}\int d^{3}k_{2}\underline{g}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},{\bf k}_{2})\underline{Y}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{Y}^{k}_{\omega_{2}}({\bf k}_{2},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int_{0}^{\infty}d\omega_{3}\int
d^{3}k_{1}\int d^{3}k_{2}\int
d^{3}k_{3}\underline{g}^{(3)}_{ijkl}(\omega_{1},\omega_{2},\omega_{3},{\bf
k},{\bf k}_{1},{\bf k}_{2},{\bf k}_{3})$
$\displaystyle\times\underline{Y}^{j}_{\omega_{1}}({\bf
k}_{1},t)\underline{Y}^{k}_{\omega_{2}}({\bf
k}_{2},t)\underline{Y}^{l}_{\omega_{3}}({\bf k}_{3},t)+....$

##  4 Canonical quantization

In order to quantize the electromagnetic field canonically, one should be
eliminated the extra degrees of freedom from the total Lagrangian (2). For
this purpose we apply the coulomb gauge and eliminate the scalar potential
$\underline{\varphi}({\bf k},t)$ from the total Lagrangian using the equation
(3). In reciprocal space the coulomb gauge is in the form ${\bf
k}\cdot\underline{{\bf A}}=0$. Applying the coulomb gauge and using the
equation (3) one can obtain the spatial Fourier transform of the scalar
potential as

$\underline{\varphi}(\vec{k},t)=-\frac{\imath\vec{k}\cdot\underline{\vec{P}}(\vec{k},t)}{\varepsilon_{0}k^{2}}$
(14)

Let ${\bf e}_{{\bf k}\lambda},\lambda=1,2$ and ${\bf e}_{{\bf
k}3}=\hat{k}=\frac{{\bf k}}{k}$ be three mutual orthogonal unit vectors for
each vector ${\bf k}$. Then from the coulomb gauge ${\bf
k}\cdot\underline{{\bf A}}=0$, we are admissible expand the vector potential
along the orthogonal unit vector ${\bf e}_{{\bf k}\lambda},\lambda=1,2$ as

$\underline{\vec{A}}(\vec{k},t)=\sum_{\lambda=1}^{2}\
\underline{A}_{\lambda}(\vec{k},t)\vec{e}_{\lambda\vec{k}}$ (15)

where $\underline{A}_{\lambda}(\vec{k},t),\lambda=1,2$ are the new dynamical
variables of electromagnetic field. Also we introduce the new dynamical
variables $\underline{X}_{\omega\lambda}({\bf k},t)$ and
$\underline{Y}_{\omega\lambda}({\bf k},t)$ for the nonlinear magnetodielectric
medium as the coefficients of the expansion of the fields $\underline{{\bf
X}}_{\omega}({\bf k},t)$ and $\underline{{\bf Y}}_{\omega}({\bf k},t)$ in
terms of the orthogonal unit vectors ${\bf e}_{{\bf k}\lambda},\lambda=1,2,3$,
that is

$\displaystyle\underline{\vec{X}}_{\omega}(\vec{k},t)$ $\displaystyle=$
$\displaystyle\sum_{\lambda=1}^{3}\underline{X}_{\omega\lambda}(\vec{k},t)\vec{e}_{\lambda\vec{k}}$
$\displaystyle\underline{\vec{Y}}_{\omega}(\vec{k},t)$ $\displaystyle=$
$\displaystyle\sum_{\lambda=1}^{3}\underline{Y}_{\omega\lambda}(\vec{k},t)\vec{e}_{\lambda\vec{k}}$
(16)

Substituting the scalar potential $\underline{\varphi}({\bf k},t)$ from (14)
in the Lagrangian (2)and using the expansions (15) and (4), one can rewrite
the Lagrangian of the total system in terms of the new coordinates of the
system as

$\displaystyle L(t)$ $\displaystyle=$ $\displaystyle\int_{0}^{\infty}d\omega\
\int^{\prime}d^{3}k\sum_{\lambda=1}^{3}\left(|\dot{\underline{X}}_{\omega\lambda}|^{2}-\omega^{2}|\underline{X}_{\omega\lambda}|^{2}+|\dot{\underline{Y}}_{\omega\lambda}|^{2}-\omega^{2}|\underline{Y}_{\omega\lambda}|^{2}\right)$
$\displaystyle+$
$\displaystyle\int^{\prime}d^{3}k\sum_{\lambda=1}^{2}\left(\varepsilon_{0}|\dot{\underline{{A}}_{\lambda}}|^{2}-\frac{k^{2}|\underline{A}_{\lambda}|^{2}}{\mu_{0}}\right)-\frac{1}{\varepsilon_{0}}\int^{\prime}d^{3}k\
\frac{|{\bf k}\cdot\underline{{\bf P}}|^{2}}{k^{2}}$ $\displaystyle+$
$\displaystyle\int^{\prime}d^{3}k\left[\left(-\sum_{\lambda=1}^{2}\underline{\dot{A}}_{\lambda}{\bf
e}_{\lambda{\bf k}}\right)\cdot\underline{{\bf P}}^{*}({\bf k},t)+\left(i{\bf
k}\times\sum_{\lambda=1}^{2}\underline{A}_{\lambda}{\bf e}_{\lambda{\bf
k}}\right)\cdot\underline{{\bf M}}^{*}({\bf k},t)+H.c\right]$

where the polarization densities $\underline{{\bf P}}({\bf
k},t),\underline{{\bf M}}({\bf k},t)$ have previously been defined in the
relations (3) and (3). The relation (4) give us the Lagrangian of the total
system in which the extra degrees of freedom have been eliminated and can be
used for a canonical quantization of the electromagnetic field together with
the nonlinear anisotropic magnetodielectric medium. The canonical conjugate
momenta of the total system for any vector ${\bf k}$ in the half space
$k_{z}\geq 0$ can be defined as

$\displaystyle-\underline{D}_{\lambda}(\vec{k},t)$ $\displaystyle=$
$\displaystyle\frac{\delta\underline{L}}{\delta\left(\underline{\dot{A}}^{*}_{\lambda}(\vec{k},t)\right)}=\varepsilon_{0}\dot{\underline{A}}_{\lambda}({\bf
k},t)-{\bf e}_{\lambda{\bf k}}\cdot\underline{{\bf P}}(\vec{k},t)$
$\displaystyle\underline{Q}_{\omega\lambda}(\vec{k},t)$ $\displaystyle=$
$\displaystyle\frac{\delta\underline{L}}{\delta\left(\underline{\dot{X}}^{*}_{\omega\lambda}(\vec{k},t)\right)}=\underline{\dot{X}}_{\omega\lambda}(\vec{k},t)$
$\displaystyle\underline{\Pi}_{\omega\lambda}(\vec{k},t)$ $\displaystyle=$
$\displaystyle\frac{\delta\underline{L}}{\delta\left(\underline{\dot{Y}}^{*}_{\omega\lambda}(\vec{k},t)\right)}=\underline{\dot{Y}}_{\omega\lambda}(\vec{k},t)$
(18)

Now, following the standard way for a canonical quantization, we impose the
following commutation relations between the conjugate dynamical variables of
the total system

$\displaystyle\left[\underline{A}^{*}_{\lambda}(\vec{k},t)\ ,\
-\underline{D}_{\lambda^{\prime}}(\vec{k^{\prime}},t)\right]$ $\displaystyle=$
$\displaystyle\imath\hbar\delta_{\lambda\lambda^{\prime}}\delta(\vec{k}-\vec{k^{\prime}})$
$\displaystyle\left[\underline{X}^{*}_{\omega\lambda}(\vec{k},t)\ ,\
\underline{Q}_{\omega^{\prime}\lambda^{\prime}}(\vec{k^{\prime}},t)\right]$
$\displaystyle=$
$\displaystyle\imath\hbar\delta_{\lambda\lambda^{\prime}}\delta(\omega-\omega^{\prime})\delta(\vec{k}-\vec{k^{\prime}})$
$\displaystyle\left[\underline{Y}^{*}_{\omega\lambda}(\vec{k},t)\ ,\
\underline{\Pi}_{\omega^{\prime}\lambda^{\prime}}(\vec{k^{\prime}},t)\right]$
$\displaystyle=$
$\displaystyle\imath\hbar\delta_{\lambda\lambda^{\prime}}\delta(\omega-\omega^{\prime})\delta(\vec{k}-\vec{k^{\prime}})$

The Hamiltonian of the total system can be written as usual way as

$\displaystyle
H(t)=\int^{\prime}d^{3}k\sum_{\lambda=1}^{2}\left[-\underline{D}_{\lambda}\underline{\dot{A}}^{*}_{\lambda}-\underline{D}^{*}_{\lambda}\underline{\dot{A}}_{\lambda}\right]$
$\displaystyle+\int_{0}^{\infty}d\omega\int^{\prime}d^{3}k\sum_{\lambda=1}^{2}\left[\underline{Q}_{\omega\lambda}\underline{\dot{X}}^{*}_{\omega\lambda}+\underline{Q}^{*}_{\omega\lambda}\underline{\dot{X}}_{\omega\lambda}\right]$
$\displaystyle+\int_{0}^{\infty}d\omega\int^{\prime}d^{3}k\sum_{\lambda=1}^{2}\left[\underline{\Pi}_{\omega\lambda}\underline{\dot{Y}}^{*}_{\omega\lambda}+\underline{\Pi}^{*}_{\omega\lambda}\underline{\dot{Y}}_{\omega\lambda}\right]-L(t)$
(20)

Using the definitions of conjugate momenta given in (4), the Hamiltonian of
the system in terms of the coordinates of the system and their conjugate
momenta becomes

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\int^{\prime}d^{3}k\sum_{\lambda=1}^{2}\left(\frac{|\underline{D}_{\lambda}-{\bf
e}_{\lambda{\bf k}}\cdot\underline{{\bf
P}}|^{2}}{\varepsilon_{0}}+\frac{k^{2}|\underline{A}_{\lambda}|^{2}}{\mu_{0}}\right)+\frac{1}{\varepsilon_{0}}\int^{\prime}d^{3}k\frac{|{\bf
k}\cdot\underline{{\bf P}}|^{2}}{k^{2}}$ $\displaystyle-$
$\displaystyle\int^{\prime}d^{3}k\left[\left(i{\bf
k}\times\sum_{\lambda=1}^{2}\underline{A}_{\lambda}{\bf e}_{\lambda{\bf
k}}\right)\cdot\underline{{\bf M}}^{*}+H.c\right]$ $\displaystyle+$
$\displaystyle\int_{0}^{\infty}d\omega\int^{\prime}d^{3}k\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)$

It can be easily shown that the Heisenberg equation for
$\underline{A}_{\lambda}({\bf k},t),\lambda=1,2$ leads to the relation
$\underline{{\bf D}}({\bf k},t)=\varepsilon_{0}\underline{{\bf E}}({\bf
k},t)+\underline{{\bf P}}({\bf k},t)$ [18], where $\underline{{\bf
D}}=\sum_{\lambda=1}^{2}\underline{D}_{\lambda}{\bf e}_{{\bf k}\lambda}$ is
the displacement field and $\underline{{\bf E}}=-\underline{\dot{{\bf
A}}}-\frac{\hat{{\bf k}}(\hat{{\bf k}}\cdot\underline{{\bf
P}})}{\varepsilon_{0}}$ is the Fourier transform of the electric field. Also
the Heisenberg equation for $\underline{D}_{\lambda},\lambda=1,2$ give us the
Maxwell equation $\underline{\dot{{\bf D}}}({\bf k},t)=i{\bf
k}\times\underline{{\bf H}}({\bf k},t)$ in the reciprocal space , where
$\underline{{\bf H}}=\frac{i{\bf k}\times\underline{{\bf
A}}}{\mu_{0}}-\underline{{\bf M}}$ is the magnetic induction field [18].

##  5 The constitutive relations of the medium

One can obtain the constitutive relations of the nonlinear magnetodielectric
medium using the definitions of polarization and magnetization densities of
the medium given by (3),(3) and the Heisenberg equations for the dynamical
variables $\underline{X}_{\omega\lambda}({\bf
k},t),\underline{Y}_{\omega\lambda}({\bf k},t),\lambda=1,2,3$.
Straightforwardly, one can show that the combination of the Heisenberg
equations of the conjugate variables $\underline{X}_{\omega\lambda}({\bf
k},t),\underline{Q}_{\omega\lambda}({\bf k},t),\lambda=1,2,3$ leads to the
following equation of motion for the components of the field $\underline{{\bf
X}}({\bf k},t)$

$\displaystyle\underline{\ddot{X}}_{\omega i}({\bf
k},t)+\omega^{2}\underline{X}_{\omega i}({\bf k},t)=\int
d^{3}k_{1}\underline{f}^{{\dagger}^{(1)}}_{ij}(\omega,{\bf k}_{1},{\bf
k})\underline{E}^{j}({\bf k}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{2}\int d^{3}k_{1}\int
d^{3}k_{2}\underline{f}^{*^{(}2)}_{jik}(\omega,\omega_{2},{\bf k}_{1},{\bf
k},{\bf k}_{2})\underline{E}^{j}({\bf k}_{1},t)\underline{X}^{*^{k}}({\bf
k}_{2},t)$ $\displaystyle+\int_{0}^{\infty}d\omega_{1}\int d^{3}k_{1}\int
d^{3}k_{2}\underline{f}^{*^{(}2)}_{kji}(\omega_{1},\omega,{\bf k}_{2},{\bf
k}_{1},{\bf k})\underline{E}^{k}({\bf k}_{2},t)\underline{X}^{*^{j}}({\bf
k}_{1},t)+...$

where the completeness relation $\sum_{\lambda=1}^{3}{\bf e}_{{\bf k}\lambda
i}{\bf e}_{{\bf k}\lambda j}=\delta_{ij}$ has been used , three points …
indicate the terms containing the coupling tensors more than the third rank
and $\underline{{\bf E}}=-\underline{\dot{{\bf A}}}-\frac{\hat{{\bf
k}}(\hat{{\bf k}}\cdot\underline{{\bf P}})}{\varepsilon_{0}}$ is the Fourier
transform of the electric field. Similarly the combination of the Heisenberg
equations of the conjugate variables
$\underline{Y}_{\omega\lambda},\underline{\Pi}_{\omega\lambda},\lambda=1,2,3$
give us

$\displaystyle\underline{\ddot{Y}}_{\omega i}({\bf
k},t)+\omega^{2}\underline{Y}_{\omega i}({\bf k},t)=\int
d^{3}k_{1}\underline{g}^{{\dagger}^{(1)}}_{ij}(\omega,{\bf k}_{1},{\bf
k})\underline{B}^{j}({\bf k}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{2}\int d^{3}k_{1}\int
d^{3}k_{2}\underline{g}^{*^{(}2)}_{jik}(\omega,\omega_{2},{\bf k}_{1},{\bf
k},{\bf k}_{2})\underline{B}^{j}({\bf k}_{1},t)\underline{Y}^{*^{k}}({\bf
k}_{2},t)$ $\displaystyle+\int_{0}^{\infty}d\omega_{1}\int d^{3}k_{1}\int
d^{3}k_{2}\underline{g}^{*^{(}2)}_{kji}(\omega_{1},\omega,{\bf k}_{2},{\bf
k}_{1},{\bf k})\underline{B}^{k}({\bf k}_{2},t)\underline{Y}^{*^{j}}({\bf
k}_{1},t)+...$

where $\underline{{\bf B}}({\bf k},t)=i{\bf k}\times\underline{{\bf A}}({\bf
k},t)$ is the Fourier transform of the magnetic field. In (5) and (5)the
components of the tensors $\underline{{\bf f}}^{{\dagger}^{(1)}}$ and
$\underline{{\bf g}}^{{\dagger}^{(1)}}$ are defined as
$\underline{f}^{{\dagger}^{(1)}}_{ij}=\underline{f}^{*^{(1)}}_{ji},\underline{g}^{{\dagger}^{(1)}}_{ij}=\underline{g}^{*^{(1)}}_{ji}$The
differential equations (5) and (5) are two frames of the complicated nonlinear
coupled differential equations. The solution of these equations exactly is
impossible unless an iteration method is used. One may apply the first order
approximation and keep only the first term in the right hand of the equations
(5) and (5) and write the solutions of the differential equations (5)and (5)
as

$\displaystyle\underline{X}_{\omega i}({\bf k},t)=\underline{X}_{N\omega
i}({\bf
k},t)+\int_{-\infty}^{t}dt^{\prime}\frac{\sin\omega(t-t^{\prime})}{\omega}\int
d^{3}k_{1}\underline{f}^{{\dagger}^{(1)}}_{ij}(\omega,{\bf k}_{1},{\bf
k})\underline{E}^{j}({\bf k}_{1},t^{\prime})$
$\displaystyle\underline{Y}_{\omega i}({\bf k},t)=\underline{Y}_{N\omega
i}({\bf
k},t)+\int_{-\infty}^{t}dt^{\prime}\frac{\sin\omega(t-t^{\prime})}{\omega}\int
d^{3}k_{1}\underline{g}^{{\dagger}^{(1)}}_{ij}(\omega,{\bf k}_{1},{\bf
k})\underline{B}^{j}({\bf k}_{1},t^{\prime})$

where $\underline{{\bf X}}_{N\omega}({\bf k},t)$ and $\underline{{\bf
Y}}_{N\omega}({\bf k},t)$ are the solutions of the homogenous equations

$\displaystyle\underline{\ddot{{\bf X}}}_{N\omega}({\bf
k},t)+\omega^{2}\underline{{\bf X}}_{N\omega}({\bf k},t)=0$
$\displaystyle\underline{\ddot{{\bf Y}}}_{N\omega}({\bf
k},t)+\omega^{2}\underline{{\bf Y}}_{N\omega}({\bf k},t)=0$ (26)

respectively which can be written as

$\displaystyle\underline{{\bf X}}_{N\omega}({\bf
k},t)=\sqrt{\frac{\hbar}{2\omega}}\sum_{\lambda=1}^{3}[b_{\lambda}(\omega,{\bf
k})e^{-i\omega t}+b^{\dagger}_{\lambda}(\omega,{\bf k})e^{i\omega t}]{\bf
e}_{{\bf k}\lambda}$ $\displaystyle\underline{{\bf Y}}_{N\omega}({\bf
k},t)=\sqrt{\frac{\hbar}{2\omega}}\sum_{\lambda=1}^{3}[d_{\lambda}(\omega,{\bf
k})e^{-i\omega t}+d^{\dagger}_{\lambda}(\omega,{\bf k})e^{i\omega t}]{\bf
e}_{{\bf k}\lambda}$ (27)

where the operators $b_{\lambda}(\omega,{\bf
k}),b^{\dagger}_{\lambda}(\omega,{\bf k})$ and $d_{\lambda}(\omega,{\bf
k}),d^{\dagger}_{\lambda}(\omega,{\bf k})$ satisfy the commutation relations

$\displaystyle[b_{\lambda}(\omega,{\bf
k}),b^{\dagger}_{\lambda^{\prime}}(\omega^{\prime},{\bf
k^{\prime}})]=\delta_{\lambda\lambda^{\prime}}\delta(\omega-\omega^{\prime})\delta({\bf
k}-{\bf k^{\prime}})$ $\displaystyle[d_{\lambda}(\omega,{\bf
k}),d^{\dagger}_{\lambda^{\prime}}(\omega^{\prime},{\bf
k^{\prime}})]=\delta_{\lambda\lambda^{\prime}}\delta(\omega-\omega^{\prime})\delta({\bf
k}-{\bf k^{\prime}})$ (28)

in compatible to the commutation relations (4).
Now substituting $\underline{{\bf X}}_{\omega}({\bf k},t)$ from (5) in (3) and
$\underline{{\bf Y}}_{\omega}({\bf k},t)$ from (5) in (3) give us the
constitutive relations of the nonlinear anisotropic magnetodielectric medium
as

$\displaystyle\underline{P}_{i}({\bf k},t)=\underline{P}_{Ni}({\bf
k},t)+\int_{-\infty}^{+\infty}dt_{1}\int
d^{3}k_{1}\chi^{{(1)}}_{ij}(t-t_{1},{\bf k},{\bf k}_{1})\underline{E}^{j}({\bf
k}_{1},t_{1})$
$\displaystyle+\int_{-\infty}^{+\infty}dt_{1}\int_{-\infty}^{+\infty}dt_{2}\int
d^{3}k_{1}\int d^{3}k_{2}\chi^{{(2)}}_{ijk}(t-t_{1},t-t_{2},{\bf k},{\bf
k}_{1},{\bf k}_{2})\underline{E}^{j}({\bf k}_{1},t_{1})\underline{E}^{k}({\bf
k_{2}},t_{2})$ $\displaystyle+......$ (29)
$\displaystyle\underline{M}_{i}({\bf k},t)=\underline{M}_{Ni}({\bf
k},t)+\int_{-\infty}^{+\infty}dt_{1}\int
d^{3}k_{1}\zeta^{{(1)}}_{ij}(t-t_{1},{\bf k},{\bf
k}_{1})\underline{B}^{j}({\bf k}_{1},t_{1})$
$\displaystyle+\int_{-\infty}^{+\infty}dt_{1}\int_{-\infty}^{+\infty}dt_{2}\int
d^{3}k_{1}\int d^{3}k_{2}\zeta^{{(2)}}_{ijk}(t-t_{1},t-t_{2},{\bf k},{\bf
k}_{1},{\bf k}_{2})\underline{B}^{j}({\bf k}_{1},t_{1})\underline{B}^{k}({\bf
k_{2}},t_{2})$ $\displaystyle+......$ (30)

where the tensors $\chi^{n}$ and $\zeta^{n}$ for $n=1,2,...$ are respectively
the electric and magnetic susceptibility tensors and define in terms of the
coupling tensors $\underline{{\bf f}}^{(n)}$ and $\underline{{\bf g}}^{(n)}$
as

$\displaystyle\chi^{(n)}_{ii_{1}i_{2}...i_{n}}(t_{1},t_{2},...,t_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{n})=$
$\displaystyle\Theta(t_{1})\Theta(t_{2})....\Theta(t_{n})\int_{0}^{\infty}d\omega_{1}\int_{0}^{d}\omega_{2}....\int_{0}^{\infty}d\omega_{n}\frac{\sin\omega_{1}t_{1}}{\omega_{1}}\frac{\sin\omega_{2}t_{2}}{\omega_{2}}....\frac{\sin\omega_{n}t_{n}}{\omega_{n}}$
$\displaystyle\times\int d^{3}p_{1}\int d^{3}p_{2}....\int
d^{3}p_{n}[\underline{f}^{(n)}_{ij_{1}j_{2}...j_{n}}(\omega_{1},\omega_{2},...,\omega_{n},{\bf
k},{\bf p}_{1},{\bf p}_{2},...{\bf p}_{n})$
$\displaystyle\times\underline{f}^{{\dagger}^{(1)}}_{j_{1}i_{1}}(\omega_{1},{\bf
k}_{1},{\bf
p}_{1})\underline{f}^{{\dagger}^{(1)}}_{j_{2}i_{2}}(\omega_{2},{\bf
k}_{2},{\bf
p}_{2})....\underline{f}^{{\dagger}^{(1)}}_{j_{n}i_{n}}(\omega_{n},{\bf
k}_{n},{\bf p}_{n})]$ (31)
$\displaystyle\zeta^{(n)}_{ii_{1}i_{2}...i_{n}}(t_{1},t_{2},...,t_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{n})=$ (32)
$\displaystyle\Theta(t_{1})\Theta(t_{2})....\Theta(t_{n})\int_{0}^{\infty}d\omega_{1}\int_{0}^{d}\omega_{2}....\int_{0}^{\infty}d\omega_{n}\frac{\sin\omega_{1}t_{1}}{\omega_{1}}\frac{\sin\omega_{2}t_{2}}{\omega_{2}}....\frac{\sin\omega_{n}t_{n}}{\omega_{n}}$
$\displaystyle\times\int d^{3}p_{1}\int d^{3}p_{2}....\int
d^{3}p_{n}[\underline{g}^{(n)}_{ij_{1}j_{2}...j_{n}}(\omega_{1},\omega_{2},...,\omega_{n},{\bf
k},{\bf p}_{1},{\bf p}_{2},...{\bf p}_{n})$
$\displaystyle\times\underline{g}^{{\dagger}^{(1)}}_{j_{1}i_{1}}(\omega_{1},{\bf
k}_{1},{\bf
p}_{1})\underline{g}^{{\dagger}^{(1)}}_{j_{2}i_{2}}(\omega_{2},{\bf
k}_{2},{\bf
p}_{2})....\underline{g}^{{\dagger}^{(1)}}_{j_{n}i_{n}}(\omega_{n},{\bf
k}_{n},{\bf p}_{n})]$

where $\Theta(t)$ is the step function and over the repeated indices
$j_{1},j_{2},...j_{n}$ should be summed. There are some symmetry conditions
that the electric and magnetic susceptibility tensors of various rank should
satisfy.The tensors $\chi^{(n)}$ and $\zeta^{(n)}$ may satisfy the symmetry
conditions [21]

$\displaystyle\chi^{(n)}_{ii_{1}i_{2}...i_{\alpha}...i_{\beta}...i_{n}}(t_{1},t_{2},...,t_{\alpha},...,t_{\beta},...,t_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\alpha},...,{\bf k}_{\beta},...,{\bf
k}_{n})$
$\displaystyle=\chi^{(n)}_{ii_{1}i_{2}...i_{\beta}...i_{\alpha}...i_{n}}(t_{1},t_{2},...,t_{\beta},...,t_{\alpha},...,t_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\beta},...,{\bf k}_{\alpha},...,{\bf
k}_{n})$
$\displaystyle\zeta^{(n)}_{ii_{1}i_{2}...i_{\alpha}...i_{\beta}...i_{n}}(t_{1},t_{2},...,t_{\alpha},...,t_{\beta},...,t_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\alpha},...,{\bf k}_{\beta},...,{\bf
k}_{n})$
$\displaystyle=\zeta^{(n)}_{ii_{1}i_{2}...i_{\beta}...i_{\alpha}...i_{n}}(t_{1},t_{2},...,t_{\beta},...,t_{\alpha},...,t_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\beta},...,{\bf k}_{\alpha},...,{\bf
k}_{n})$

From the definitions (5) and (32) it is clear that the symmetry relations (5)
and (5) are satisfied provided that we impose the conditions

$\displaystyle\underline{f}^{(n)}_{ij_{1}j_{2}...j_{\alpha}...j_{\beta}...j_{n}}(\omega_{1},\omega_{2},...,\omega_{\alpha},...,\omega_{\beta},...,\omega_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\alpha},...,{\bf k}_{\beta},...,{\bf
k}_{n})$
$\displaystyle=\underline{f}^{(n)}_{ij_{1}j_{2}...j_{\beta}...j_{\alpha}...j_{n}}(\omega_{1},\omega_{2},...,\omega_{\beta},...,\omega_{\alpha},...,\omega_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\beta},...,{\bf k}_{\alpha},...,{\bf
k}_{n})$
$\displaystyle\underline{g}^{(n)}_{ij_{1}j_{2}...j_{\alpha}...j_{\beta}...j_{n}}(\omega_{1},\omega_{2},...,\omega_{\alpha},...,\omega_{\beta},...,\omega_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\alpha},...,{\bf k}_{\beta},...,{\bf
k}_{n})$
$\displaystyle=\underline{g}^{(n)}_{ij_{1}j_{2}...j_{\beta}...j_{\alpha}...j_{n}}(\omega_{1},\omega_{2},...,\omega_{\beta},...,\omega_{\alpha},...,\omega_{n},{\bf
k},{\bf k}_{1},{\bf k}_{2},...,{\bf k}_{\beta},...,{\bf k}_{\alpha},...,{\bf
k}_{n})$

on the coupling tensors $\underline{{\bf f}}^{(n)}$ and $\underline{{\bf
g}}^{(n)}$ for $n=2,3,...$.
In constitutive relations (5) and (5) $\underline{{\bf P}}_{N}({\bf k},t)$ and
$\underline{{\bf M}}_{N}({\bf k},t)$ are the noise polarization and
magnetization fields, respectively, where their components are written in
terms of the $\underline{{\bf X}}_{N\omega}({\bf k},t)$ and $\underline{{\bf
Y}}_{N\omega}({\bf k},t)$ as

$\displaystyle\underline{P}_{Ni}({\bf k},t)=\int_{0}^{\infty}d\omega_{1}\int
d^{3}k_{1}\underline{f}^{(1)}_{ij}(\omega_{1},{\bf k},{\bf
k}_{1})\underline{X}_{N\omega_{1}}^{j}({\bf k}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}k_{1}\int d^{3}k_{2}[\underline{f}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},{\bf k}_{2})\underline{X}_{N\omega_{1}}^{j}({\bf
k}_{1},t)\underline{X}_{N\omega_{2}}^{k}({\bf k}_{2},t)]$
$\displaystyle\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}k_{1}\int
d^{3}k_{2}\left[\underline{f}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf k},{\bf
k}_{1},{\bf
k}_{2})\int_{0}^{t}dt^{\prime}\frac{\sin\omega_{2}(t-t^{\prime})}{\omega_{2}}\right.$
$\displaystyle\times\left.\int
d^{3}p_{2}\underline{f}^{{\dagger}^{(1)}}_{kl}(\omega_{2},{\bf p}_{2},{\bf
k}_{2})\left(\underline{X}_{N\omega_{1}}^{j}({\bf
k}_{1},t)\underline{E}^{l}({\bf p}_{2},t^{\prime})+\underline{E}^{l}({\bf
p}_{2},t^{\prime})\underline{X}^{j}_{N\omega_{1}}({\bf
k}_{1},t)\right)\right]+...$ $\displaystyle\underline{M}_{Ni}({\bf
k},t)=\int_{0}^{\infty}d\omega_{1}\int
d^{3}k_{1}\underline{g}^{(1)}_{ij}(\omega_{1},{\bf k},{\bf
k}_{1})\underline{Y}_{N\omega_{1}}^{j}({\bf k}_{1},t)$
$\displaystyle+\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}k_{1}\int d^{3}k_{2}[\underline{g}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf
k},{\bf k}_{1},{\bf k}_{2})\underline{Y}_{N\omega_{1}}^{j}({\bf
k}_{1},t)\underline{Y}_{N\omega_{2}}^{k}({\bf k}_{2},t)]$
$\displaystyle\int_{0}^{\infty}d\omega_{1}\int_{0}^{\infty}d\omega_{2}\int
d^{3}k_{1}\int
d^{3}k_{2}\left[\underline{g}^{(2)}_{ijk}(\omega_{1},\omega_{2},{\bf k},{\bf
k}_{1},{\bf
k}_{2})\int_{0}^{t}dt^{\prime}\frac{\sin\omega_{2}(t-t^{\prime})}{\omega_{2}}\right.$
$\displaystyle\times\left.\int
d^{3}p_{2}\underline{g}^{{\dagger}^{(1)}}_{kl}(\omega_{2},{\bf p}_{2},{\bf
k}_{2})\left(\underline{Y}_{N\omega_{1}}^{j}({\bf
k}_{1},t)\underline{B}^{l}({\bf p}_{2},t^{\prime})+\underline{B}^{l}({\bf
p}_{2},t^{\prime})\underline{Y}^{j}_{N\omega_{1}}({\bf
k}_{1},t)\right)\right]+...$

where the symmetry conditions (5) and (5) have been used , the summation
should be done on the repeated indices and and the three points … indicate the
terms containing the coupling tensors more than the third rank.

##  6 The time dependence of the electric and magnetic field

Let we write the the electric field $\underline{{\bf E}}({\bf k},t)$ in terms
of its temporal Fourier transform as

$\underline{{\bf E}}({\bf
k},t)=\frac{1}{(2\pi)^{\frac{3}{2}}}\int_{-\infty}^{+\infty}d\omega\underline{\tilde{{\bf
E}}}({\bf k},\omega)e^{i\omega t}$ (39)

then, using this relation and similar transformations for the other fields and
combination of the Maxwell equations $i{\bf k}\times{\bf E}({\bf
k},t)=\dot{\underline{{\bf B}}}({\bf k},t),\dot{\underline{{\bf D}}}({\bf
k},t)=i{\bf k}\times\underline{{\bf H}}({\bf k},t)$ give us

${\bf k}\times{\bf k}\times\underline{\tilde{{E}}}({\bf
k},\omega)+\frac{\omega^{2}}{c^{2}}\underline{\tilde{{\bf E}}}({\bf
k},\omega)=-\mu_{0}\omega^{2}\underline{\tilde{{\bf P}}}({\bf
k},\omega)-\mu_{0}\omega{\bf k}\times\underline{\tilde{{\bf M}}}({\bf
k},\omega)$ (40)

for the electric field in frequency domain where $\underline{\tilde{{\bf
P}}}({\bf k},\omega)$ and $\underline{\tilde{{\bf M}}}({\bf k},\omega)$ are
the temporal Fourier transforms of $\underline{{\bf P}}({\bf k},t)$ and
$\underline{{\bf M}}({\bf k},t)$, respectively. Now substituting
$\underline{\tilde{{\bf P}}}({\bf k},\omega)$ and $\underline{\tilde{{\bf
M}}}({\bf k},\omega)$, using the constitutive relations (5) and (5), in (40)
leads to the integral equation

$\displaystyle\Lambda_{ii_{1}}({\bf
k},\omega)\underline{\tilde{E}}_{i_{1}}({\bf
k},\omega)+(2\pi)^{\frac{3}{2}}\mu_{0}\omega^{2}\int
d^{3}k_{1}\tilde{\chi}^{(1)}_{ii_{1}}(\omega,{\bf k},{\bf
k}_{1})\tilde{E}_{i_{1}}({\bf k}_{1},\omega)$
$\displaystyle+(2\pi)^{\frac{3}{2}}\mu_{0}\omega^{2}\int d^{3}k_{1}\int
d^{3}k_{2}\int_{0}^{\infty}d\omega_{1}$
$\displaystyle\times[\tilde{\chi}^{(2)}_{ii_{1}i_{2}}(\omega_{1},\omega-\omega_{1},{\bf
k},{\bf k}_{1},{\bf k}_{2})\tilde{E}_{i_{1}}({\bf
k}_{1},\omega_{1})\tilde{E}_{i_{2}}({\bf k}_{2},\omega-\omega_{1})]$
$\displaystyle+......$ $\displaystyle+(2\pi)^{\frac{3}{2}}\mu_{0}\omega\int
d^{3}k_{1}\epsilon_{ij_{1}j_{2}}k_{j_{1}}\tilde{\zeta}^{(1)}_{j_{2}i_{1}}(\omega,{\bf
k},{\bf k}_{1})\frac{({\bf k}_{1}\times\tilde{E}({\bf
k}_{1},\omega))_{i_{1}}}{-\omega}$
$\displaystyle+(2\pi)^{\frac{3}{2}}\mu_{0}\omega\int d^{3}k_{1}\int
d^{3}k_{2}\int_{0}^{\infty}d\omega_{1}$
$\displaystyle\times\left[\epsilon_{ij_{1}j_{2}}k_{j_{1}}\tilde{\zeta}^{(2)}_{j_{2}i_{1}i_{2}}(\omega_{1},\omega-\omega_{1},{\bf
k},{\bf k}_{1},{\bf k}_{2})\right.$ $\displaystyle\times\left.\frac{({\bf
k}_{1}\times\tilde{E}({\bf
k}_{1},\omega_{1}))_{i_{1}}}{-\omega_{1}}\frac{({\bf
k}_{2}\times\tilde{E}({\bf
k}_{2},\omega-\omega_{1}))_{i_{2}}}{-(\omega-\omega_{1})}\right]+...=\underline{\tilde{J}}_{i}({\bf
k},\omega)$ (41)

for the electric field where the Maxwell equation $i{\bf
k}\times\underline{{\bf E}}({\bf k},t)=\underline{\dot{{\bf B}}}({\bf k},t)$
has been used, $\epsilon_{ijk}$ is three dimensional levi-Civita symbol and

$\displaystyle\Lambda_{ii_{1}}({\bf
k},\omega)=\epsilon_{ijm}\epsilon_{mni_{1}}k_{j}k_{n}+\frac{\omega^{2}}{c^{2}}$
$\displaystyle\underline{\tilde{J}}_{i}({\bf
k},\omega)=-\mu_{0}\omega^{2}\underline{\tilde{P}}_{Ni}({\bf
k},\omega)-\mu_{0}\omega\epsilon_{ijl}k_{j}\underline{\tilde{M}}_{Nl}({\bf
k},\omega)$

In integral equation (6) $\tilde{\chi}^{(n)},\tilde{\zeta}^{(n)},n=1,2,...$
are the temporal Fourier transforms of the Susceptibility tensors
$\chi^{(n)},\zeta^{(n)}$ which is defined as

$\displaystyle\chi^{(n)}_{ii_{1}i_{2}...i_{n}}(t_{1},t_{2},...,t_{n},{\bf
k},{\bf k}_{1},...,{\bf
k}_{n})=\frac{1}{(2\pi)^{\frac{3(n+1)}{2}}}\int_{-\infty}^{+\infty}d\omega_{1}\int_{-\infty}^{+\infty}d\omega_{2}...\int_{-\infty}^{+\infty}d\omega_{n}$
$\displaystyle\times[\chi^{(n)}_{ii_{1}...i_{n}}(\omega_{1},\omega_{2},...,\omega_{n},{\bf
k}_{1},{\bf k}_{2},...,{\bf
k}_{n})e^{i\omega_{1}t_{1}+i\omega_{2}t_{2}+...i\omega_{n}t_{n}}]$

and $\tilde{\zeta}^{(n)}$ is defined in a similar way.
The equation (6) is an integral equation with kerenels $\tilde{\chi}^{(n)}$
and $\tilde{\zeta}^{(n)}$ and the source term $\underline{\tilde{{\bf
J}}}({\bf k},\omega)$. As well, from the definitions of the noise polarization
and magnetization densities given by (5) and(5), it is seen that the source
term $\underline{\tilde{{\bf J}}}({\bf k},\omega)$ is dependent on the
electric and magnetic fields. The integral equations such as (6) can be solved
by an iteration method[22]. In the zero order approximation we neglect the
integrals in the left hand of the equation(6) and write the components of the
electric field $\underline{\tilde{{\bf E}}}({\bf k},\omega)$ as follows

$\underline{\tilde{E}}^{(0)}_{i}({\bf k},\omega)=\Lambda^{-1}_{ij}({\bf
k},\omega)\underline{\tilde{J}}^{(0)}_{j}({\bf k},\omega)$ (44)

where $\underline{\tilde{{\bf J}}}^{(0)}({\bf k},\omega)$ is the source term
$\underline{\tilde{{\bf J}}}({\bf k},\omega)$ given by (6) from which the
terms containing the electric and magnetic fields have been eliminated.
Therefore in the zero order approximation the time dependence of the electric
and magnetic fields can be written as

$\displaystyle\underline{{\bf E}}^{(0)}({\bf
k},t)=\frac{1}{(2\pi)^{\frac{3}{2}}}\int_{-\infty}^{+\infty}d\omega\underline{\tilde{{\bf
E}}}^{(0)}({\bf k},\omega)e^{i\omega t}$ $\displaystyle\underline{{\bf
B}}^{(0)}({\bf
k},t)=\frac{1}{(2\pi)^{\frac{3}{2}}}\int_{-\infty}^{+\infty}d\omega\frac{{\bf
k}\times\underline{\tilde{{\bf E}}}^{(0)}({\bf k},\omega)}{-\omega}e^{i\omega
t}$ (45)

In $n^{,}$th order approximation the components of the temporal Fourier
transform of the electric field for $n=1,2,3,...$ can be computed using the
following recurrence relation

$\displaystyle\underline{\tilde{E}}^{(n)}_{i}({\bf
k},\omega)=-(2\pi)^{\frac{3}{2}}\mu_{0}\omega^{2}\int
d^{3}k_{1}\left[\Lambda^{(-1)}_{ij}({\bf
k},\omega)\tilde{\chi}^{(1)}_{ji_{1}}(\omega,{\bf k},{\bf
k}_{1})\tilde{E}^{(n-1)}_{i_{1}}({\bf k}_{1},\omega)\right]$
$\displaystyle-(2\pi)^{\frac{3}{2}}\mu_{0}\omega^{2}\int d^{3}k_{1}\int
d^{3}k_{2}\int_{0}^{\infty}d\omega_{1}$
$\displaystyle\times\left[\Lambda^{(-1)}_{ij}({\bf
k},\omega)\tilde{\chi}^{(2)}_{ji_{1}i_{2}}(\omega_{1},\omega-\omega_{1},{\bf
k},{\bf k}_{1},{\bf k}_{2})\tilde{E}^{(n-1)}_{i_{1}}({\bf
k}_{1},\omega_{1})\tilde{E}^{(n-1)}_{i_{2}}({\bf
k}_{2},\omega-\omega_{1})\right]+...$
$\displaystyle-(2\pi)^{\frac{3}{2}}\mu_{0}\omega\int
d^{3}k_{1}\left[\Lambda^{(-1)}_{ij}({\bf
k},\omega)\epsilon_{jj_{1}j_{2}}k_{j_{1}}\tilde{\zeta}^{(1)}_{j_{2}i_{1}}(\omega,{\bf
k},{\bf k}_{1})\frac{[{\bf k}_{1}\times\tilde{E}^{(n-1)}({\bf
k}_{1},\omega)]_{i_{1}}}{-\omega}\right]$
$\displaystyle-(2\pi)^{\frac{3}{2}}\mu_{0}\omega\int d^{3}k_{1}\int
d^{3}k_{2}\int_{0}^{\infty}d\omega_{1}$
$\displaystyle\times\left[\Lambda^{(-1)}_{ij}({\bf
k},\omega)\epsilon_{jj_{1}j_{2}}k_{j_{1}}\tilde{\zeta}^{(2)}_{j_{2}i_{1}i_{2}}(\omega_{1},\omega-\omega_{1},{\bf
k},{\bf k}_{1},{\bf k}_{2})\right.$ $\displaystyle\times\left.\frac{[{\bf
k}_{1}\times\tilde{E}^{(n-1)}({\bf
k}_{1},\omega_{1})]_{i_{1}}}{-\omega_{1}}\frac{[{\bf
k}_{2}\times\tilde{E}^{(n-1)}({\bf
k}_{2},\omega-\omega_{1})]_{i_{2}}}{-(\omega-\omega_{1})}\right]+...$
$\displaystyle+\Lambda^{(-1)}_{ij}({\bf
k},\omega)\underline{\tilde{J}}^{(n-1)}_{j}({\bf k},\omega)$ (46)

where the three points …indicate the terms containing the electric and
magnetic susceptibility tensors more than the third rank and the components of
$\underline{\tilde{{\bf J}}}^{(n-1)}({\bf k},\omega)$ are defined as

$\displaystyle\underline{\tilde{J}}^{(n-1)}_{i}({\bf
k},\omega)=-\mu_{0}\omega^{2}\underline{\tilde{P}}^{(n-1)}_{Ni}({\bf
k},\omega)-\mu_{0}\omega\epsilon_{ijl}k_{j}\underline{\tilde{M}}^{(n-1)}_{Nl}({\bf
k},\omega)$ (47)

and $\underline{\tilde{{\bf P}}}_{N}^{(n-1)}({\bf k},\omega)$ and
$\underline{{\tilde{\bf M}}}_{N}^{(n-1)}({\bf k},\omega)$ are the temporal
Fourier transforms of the noise polarization and magnetization fields given by
(5) and (5) in which the electric and magnetic fields $\underline{{\bf
E}}^{(n-1)}({\bf k},t),\underline{{\bf B}}^{(n-1)}({\bf k},t)$ are inserted in
those terms of $\underline{{\bf P}}_{N}({\bf k},t),\underline{{\bf
M}}_{N}({\bf k},t)$ that are dependent on the electric and magnetic fields.

##  7 Conclusion

By introducing a Lagrangian for the total system , that is the nonlinear
anisotropic magnetodielectric medium plus the electromagnetic field, a
canonical quantization was introduced for both the medium and the
electromagnetic field. The coupling tensors appeared in the Lagrangian of the
total system had a crucial role in this quantization scheme. So that the
electric and magnetic polarization fields could be written in terms of the
coupling tensors and the space dependent three dimensional harmonic
oscillators modeling the medium. As well the electric and magnetic
susceptibility tensors of the medium were obtained in terms of the coupling
tensors. Also the constitutive relations of the magnetodielectric medium
successfully could be obtained using the Heisenberg equations of the harmonic
oscillators modeling the medium.

## References

  * [1] H. B. G. Casimir, Proc. K. Ned. Akad. Wet.51,793(1948).
  * [2] P. W. Milonni, _The quantum vacuum: An introduction to quantum electrodynamics_ (Academic Press, San Diego,1994).
  * [3] H. Casimir, D. Polder, Phys. Rev. A 73, 360(1948).
  * [4] E. Lifshitz, _The theory of molecular attractive forces betweensolids_(1956).
  * [5] I. E. e. Dzyaloshinskii, E. Lifshitz, L. P. Pitaevskii,_General theory of Van der waals forces_ , Physics-Uspekhi 4,153-176(1961).
  * [6] C. Raabe, D. G. Welsch, Phys. Rev. A 71, 013814 (2005).
  * [7] R. Matloob, Phys. Rev. A 70, 062110(2004).
  * [8] R. Matloob, Phys. Rev. A 60, 3421(1999).
  * [9] R. Podgornik, P. L. Hansen, V. A. Parsegian, J. Chem. Phys. Vol. 119, No. 2(2003).
  * [10] R. Golestanian, Phys. Rev. Lett. 95, 230601 (2005).
  * [11] E. M. purcell, Phys. Rev. 69, 681 (1946).
  * [12] E. Yablonovitch, Phys. Rev. Lett. 58,2059(1987).
  * [13] S. M. Barnett, B. Huttner, R. Loudon, Phys. Rev. Lett.68,3698(1992).
  * [14] S. M. Barnett, B. Huttner, R. Loudon, R. Matloob, J. Phys. B 29, 3763(1996).
  * [15] B. Huttner, S. M. Barnett, Phys. Rev. A 46, 4306(1992)
  * [16] H. T. Dung, L. Knöll, D. G. Welsch, Phys. Rev. A 57, 3931(1998).
  * [17] F. Kheirandish, M. Amooshahi, M. Soltani, J. Phys. B: At. Mol. Opt. Phys. 42, 075504(2009).
  * [18] M. Amooshahi, J. Math. Phys. 50, 062301(2009)
  * [19] M. Amooshahi, Eur. Phys. J. D 69 (2015)
  * [20] M. Amooshahi, E. Amooghorban, Annals of physics, 325, 1976-1986 (2010)
  * [21] Guang S. He, Song H. Liu, _Physics of nonlinear optics_ , World Scientific (1999).
  * [22] Sadri Hassani, _Foundations of Mathematical Physics_ , McGraw-Hill )1991)