Quantum theory of plasmon

The existence of the resonance in collective excitations of the electron gas in metals was demonstrated in early works on plasma oscillations in solids (see, for example, [1–4]). It behaves like a quasiparticle of a special type called plasmon. Recently the research on the fundamental processes with the participation of plasmons has led to the emergence of a new scientific discipline —plasmonics. The variety of plasmonic processes and phenomena is quite broad: the formation of hybrid systems consisting of semiconductor quantum dot and metallic nanoparticles [5, 6], interaction between a metal nanoparticle and a dipole emitter [7], exiton–plasmon coupling (plexciton) [8–10], plasmon resonance energy transfer (PRET) [11, 12], plasmon-enhanced light absorption [13, 14] and fluorescene [15–20], plasmonic-molecular resonance [21–28] etc. The results of the research on plasmonic processes have led to the creation of plasmonic nanoantenae for various efficient applications [29].

In many of the above-mentioned fundamental research works on plasmonic processes there was the need to use the Hamiltonians of the interactions between plasmon and other elementary excitations in matter. All those interaction Hamiltonians were introduced in a phenomenological manner. In order to exactly derive the Hamiltonians of the interactions of plasmon from first principles it is necessary to elaborate a rigorous procedure for quantizing the collective excitations of electron gas. This task will be performed in the present work. The main ideas were outlined in our previous publication [30].

Consider a simple model of metal consisting of a gas of itinerant electrons freely moving inside metal over a background of ions with homogeneous distribution of the positive charge. Denote $n(\mathbf{r},t)$ the electron density (number of electrons per unit volume) and n₀ its mean value (averaged over space and time). The average charge density −en₀ of electron, −e being the electron charge, compensates the average positive charge of the ions in the background, and the fluctuating charge density in the metal is

$$\begin{eqnarray}&&\rho (\mathbf{r},t)=-e[n(\mathbf{r},t)-{{n}_{0}}].\end{eqnarray} \tag{ 1 }$$

According to the Coulomb law the effective charge distribution (1) creates a time-dependent electrical field with the potential

$$\begin{eqnarray}&&\phi (\mathbf{r},t)=\mathop{\int }^{}d\mathbf{r}\prime \frac{\rho (\mathbf{r}\prime ,t)}{\left| \mathbf{r}-\mathbf{r}\prime \right|}\cdot \end{eqnarray} \tag{ 2 }$$

From formula (2) there follows the Poisson equation

$$\begin{eqnarray}&&{{\nabla }^{2}}\phi (\mathbf{r},t)=-4\pi \rho (\mathbf{r},t).\end{eqnarray} \tag{ 3 }$$

The mutual interaction between effective charge densities at two different regions in the space gives rise to the potential energy of the electron gas

$$\begin{eqnarray}&&U(t)=\frac{1}{2}\mathop{\int }^{}d\mathbf{r}\mathop{\int }^{}d\mathbf{r}\prime \rho (\mathbf{r},t)\frac{1}{\left| \mathbf{r}-\mathbf{r}\prime \right|}\rho (\mathbf{r}\prime ,t),\end{eqnarray} \tag{ 4 }$$

which can be also written in the form

$$\begin{eqnarray}&&U(t)=\frac{1}{2}\mathop{\int }^{}d\mathbf{r}\rho (\mathbf{r},t)\phi (\mathbf{r},t).\end{eqnarray} \tag{ 5 }$$

As a consequence of the oscillating displacements of electrons, the fluctuation of the electron density $n(\mathbf{r},t)$ generates the total kinetic energy of the electron gas. Denote $\delta \mathbf{r}(\mathbf{r},t)$ the displacement of the electron having the coordinate r at the time moment t, and m the electron mass. Since the electron has the velocity

$$\begin{eqnarray}&&\delta \mathbf{\dot{r}}(\mathbf{r},t)=\frac{\partial \delta \mathbf{r}(\mathbf{r},t)}{\partial t},\end{eqnarray} \tag{ 6 }$$

the whole electron gas has following total kinetic energy

$$\begin{eqnarray}&&T(t)=\frac{1}{2}m\mathop{\int }^{}d\mathbf{r}n(\mathbf{r},t)\delta \mathbf{\dot{r}}{{(\mathbf{r},t)}^{2}}.\end{eqnarray} \tag{ 7 }$$

Thus we have derived following expression of the Lagrangian of the electron gas

$$\begin{eqnarray}\begin{array}{*{35}{l}} L(t) & = & \frac{1}{2}m\mathop{\int }^{}d\mathbf{r}n(\mathbf{r},t)\delta \mathbf{\dot{r}}{{(\mathbf{r},t)}^{2}} \\ {} & {} & -\frac{1}{2}\mathop{\int }^{}d\mathbf{r}\mathop{\int }^{}d\mathbf{r}\prime \rho (\mathbf{r},t)\frac{1}{\left| \mathbf{r}-\mathbf{r}\prime \right|}\rho (\mathbf{r}\prime ,t), \\ \end{array}\end{eqnarray} \tag{ 8 }$$

where the spatial integrations are performed over the whole volume of the electron gas. Replacing the expression

$$\begin{eqnarray*}&&n(\mathbf{r},t)={{n}_{0}}-\frac{1}{e}\rho (\mathbf{r},t)\end{eqnarray*}$$

into the r.h.s. of formula (8), we rewrite this formula in the form containing only two types of dynamical variables—the electron displacement vectors $\delta \mathbf{r}(\mathbf{r},t)$ and the charge density $\rho (\mathbf{r},t)$:

$$\begin{eqnarray}\begin{array}{*{35}{l}} L(t) & = & \frac{1}{2}m{{n}_{0}}\mathop{\int }^{}d\mathbf{r}\delta \mathbf{\dot{r}}{{(\mathbf{r},t)}^{2}} \\ {} & {} & -\frac{1}{2}\mathop{\int }^{}d\mathbf{r}\mathop{\int }^{}d\mathbf{r}\prime \rho (\mathbf{r},t)\frac{1}{\left| \mathbf{r}-\mathbf{r}\prime \right|}\rho (\mathbf{r}\prime ,t) \\ {} & {} & -\frac{1}{2}\frac{m}{e}\mathop{\int }^{}d\mathbf{r}\rho (\mathbf{r},t)\delta \mathbf{\dot{r}}{{(\mathbf{r},t)}^{2}}. \\ \end{array}.\end{eqnarray} \tag{ 9 }$$

The dynamical variables $\delta \mathbf{r}(\mathbf{r},t)$ and $\rho (\mathbf{r},t)$ cannot be completely independent, and we must establish the relationship between these physical quantities. We note that since oscillating displacements of electrons cause the fluctuations of the electron number $n(\mathbf{r},t)$, there must exist some direct relationship between $n(\mathbf{r},t)$ and $\delta \mathbf{r}(\mathbf{r},t)$. For establishing this relationship we consider any finite volume Ω bounded by a closed surface Σ in the spatial region of the electron gas. In comparison with the average electron number density n₀ needed for compensating the positive charge of the ions in the lattice of the metal, the number of excess electrons in the volume Ω is

$$\begin{eqnarray}&&{{N}_{\Omega }}(t)=\underset{\Omega }{\mathop \int }\,d\mathbf{r}[n(\mathbf{r},t)-{{n}_{0}}].\end{eqnarray} \tag{ 10 }$$

Due to the fluctuation of $n(\mathbf{r},t)$, this number changes with the time and its increment during a very short time interval (t, t + δt) is

$$\begin{eqnarray}&&\delta {{N}_{\Omega }}(t)=\underset{\Omega }{\mathop \int }\,d\mathbf{r}\dot{n}(\mathbf{r},t)\delta t,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\dot{n}(\mathbf{r},t)=\frac{\partial n(\mathbf{r},t)}{\partial t}\cdot \end{eqnarray} \tag{ 12 }$$

On the other side, the fluctuation of $n(\mathbf{r},t)$ is caused by the oscillating displacements of electrons. Denote $\delta {{N}_{\Sigma }}(t)$ the number of electrons displacing across the boundary Σ of the volume Ω and leaving this volume, i.e. moving from the inside of the closed surface Σ to its outside during the same time interval (t, t + δt). We have

$$\begin{eqnarray}&&\delta {{N}_{\Sigma }}(t)=\underset{\Omega }{\mathop \oint }\,d\mathbf{S}\delta \mathbf{\dot{r}}(\mathbf{r},t)n(\mathbf{r},t)\delta t.\end{eqnarray} \tag{ 13 }$$

By means of the Ostrogradski–Gauss theorem we transform the surface integral over Σ in the r.h.s. of formula (13) into a volume integral over Ω and obtain

$$\begin{eqnarray}&&\delta {{N}_{\Sigma }}(t)=\underset{\Omega }{\mathop \int }\,d\mathbf{r}\nabla \left[ n(\mathbf{r},t)\delta \mathbf{\dot{r}}(\mathbf{r},t) \right]\delta t.\end{eqnarray} \tag{ 14 }$$

Because the number $\delta {{N}_{\Sigma }}(t)$ of electrons leaving the volume Ω across its surface Σ must be equal to the decrement $-\delta {{N}_{\Omega }}(t)$ of the number of electrons contained in the volume Ω

$$\begin{eqnarray*}&&\delta {{N}_{\Sigma }}(t)=-\delta {{N}_{\Omega }}(t),\end{eqnarray*}$$

we obtain following formula

$$\begin{eqnarray}\begin{array}{*{35}{l}} \dot{n}(\mathbf{r},t) & = & \frac{\partial n(\mathbf{r},t)}{\partial t} \\ {} & = & -\nabla \left[ n(\mathbf{r},t)\delta \mathbf{\dot{r}}(\mathbf{r},t) \right] \\ \end{array}\end{eqnarray} \tag{ 15 }$$

called the continuity equation. In terms of the charge density (1), the continuity equation has the form

$$\begin{eqnarray}\begin{array}{*{35}{l}} \dot{\rho }(\mathbf{r},t) & = & \frac{\partial \rho (\mathbf{r},t)}{\partial t} \\ {} & = & e{{n}_{0}}\nabla \delta \mathbf{\dot{r}}(\mathbf{r},t) \\ {} & {} & -\nabla \left[ \rho (\mathbf{r},t)\delta \mathbf{\dot{r}}(\mathbf{r},t) \right]. \\ \end{array}\end{eqnarray} \tag{ 16 }$$

Thus the dynamical variables $\delta \mathbf{r}(\mathbf{r},t)$ and $\rho (\mathbf{r},t)$ in the Lagrangian (9) are not independent. They must satisfy the subsidiary condition (16).

The appearance of the small effective charge density $\rho (\mathbf{r},t)$ in comparison with the average charge density (having the absolute value en₀) of the electron gas is the consequence of the very small oscillating displacements $\delta \mathbf{r}(\mathbf{r},t)$ of the electrons. In the Lagrangian (9) and the subsidiary condition (16) they are two very small quantities of the same order. Consider Lagrangian (9) and the subsidiary condition (16) in the lowest order with respect to these very small quantities. Then the Lagrangian has the following approximate expression

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{0}}(t) & = & \frac{1}{2}m{{n}_{0}}\mathop{\int }^{}d\mathbf{r}\delta \mathbf{\dot{r}}{{(\mathbf{r},t)}^{2}} \\ {} & {} & -\frac{1}{2}\mathop{\int }^{}d\mathbf{r}\mathop{\int }^{}d\mathbf{r}\prime \rho (\mathbf{r},t)\frac{1}{\mathbf{r}-\mathbf{r}\prime }\rho (\mathbf{r}\prime ,t), \\ \end{array}\end{eqnarray} \tag{ 17 }$$

and the subsidiary condition becomes

$$\begin{eqnarray}&&\dot{\rho }(\mathbf{r},t)=e{{n}_{0}}\nabla \delta \mathbf{\dot{r}}(\mathbf{r},t).\end{eqnarray} \tag{ 18 }$$

Now we establish the canonical formalism of classical mechanics and then apply the quantization procedure to the study of the electron gas with the charge density $\rho (\mathbf{r},t)$ and the electron displacement vectors $\delta \mathbf{r}(\mathbf{r},t)$ being the solution of the system of differential equations consisting of Lagrange equations with the approximate Lagrangian (17) and the subsidiary condition (18). For this purpose we decompose functions $\rho (\mathbf{r},t)$ and $\delta \mathbf{r}(\mathbf{r},t)$ into the Fourier series of plane waves normalized in a cube with the volume V and orthogonalized by means of the periodic boundary conditions:

$$\begin{eqnarray}&&\rho (\mathbf{r},t)=\frac{1}{\sqrt{V}}\underset{\mathbf{k}}{\mathop \sum }\,{{\rho }_{\mathbf{k}}}(t){{e}^{i\mathbf{kr}}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\delta \mathbf{r}(\mathbf{r},t)=\frac{1}{\sqrt{V}}\underset{\mathbf{k}}{\mathop \sum }\,\underset{\alpha =1}{\overset{3}{\mathop \sum }}\,\mathbf{e}_{\mathbf{k}}^{(\alpha )}f_{\mathbf{k}}^{(\alpha )}(t){{e}^{i\mathbf{kr}}}\end{eqnarray} \tag{ 20 }$$

where $\mathbf{e}_{\mathbf{k}}^{(\alpha )}$ for each pair of indices k and α, α = 1, 2, 3, are three real unit vectors satisfying the orthogonalization and normalization conditions

$$\begin{eqnarray}&&\mathbf{e}_{\mathbf{k}}^{(\alpha )}\mathbf{e}_{\mathbf{k}}^{(\alpha )}={{\delta }_{\alpha \beta }}\end{eqnarray} \tag{ 21 }$$

and having properties

$$\begin{eqnarray}\begin{matrix} \mathbf{ke}_{\mathbf{k}}^{(1,2)} & = & 0, \\ \mathbf{e}_{\mathbf{k}}^{(3)} & = & \frac{\mathbf{k}}{k}, \\ &&\end{matrix}\;\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\mathbf{e}_{-k}^{(\alpha )}=-\mathbf{e}_{\mathbf{k}}^{(\alpha )}.\end{eqnarray} \tag{ 23 }$$

Since $\rho (\mathbf{r},t)$ and $\delta \mathbf{r}(\mathbf{r},t)$ are the real quantities, their Fourier transforms must satisfy conditions

$$\begin{eqnarray}\begin{matrix} {{\rho }_{-\mathbf{k}}}(t) & = & {{\rho }_{\mathbf{k}}}(t)*, \\ f_{-\mathbf{k}}^{(\alpha )}(t) & = & f_{\mathbf{k}}^{(\alpha )}(t)*. \\ &&\end{matrix}\end{eqnarray} \tag{ 24 }$$

In terms of the Fourier components ${{\rho }_{\mathbf{k}}}(t)$ and $f_{\mathbf{k}}^{(\alpha )}(t)$ the subsidiary condition (18) becomes

$$\begin{eqnarray}&&{{\dot{\rho }}_{\mathbf{k}}}(t)=\frac{d{{\rho }_{\mathbf{k}}}(t)}{dt}=ie{{n}_{0}}k\dot{f}_{\mathbf{k}}^{(3)}(t).\end{eqnarray} \tag{ 25 }$$

Substituting the Fourier series (19) and (20) into the r.h.s. of formula (17), we obtain

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{0}}(t) & = & \frac{1}{2}\underset{\mathbf{k}}{\mathop \sum }\,\left[ m{{n}_{0}}\underset{\alpha =1}{\overset{3}{\mathop \sum }}\,\dot{f}_{\mathbf{k}}^{(\alpha )}(t)*\dot{f}_{\mathbf{k}}^{(\alpha )}(t) \right. \\ {} & {} & \left. -\frac{4\pi }{{{k}^{2}}}{{\rho }_{\mathbf{k}}}(t)*{{\rho }_{\mathbf{k}}}(t) \right]. \\ \end{array}\end{eqnarray} \tag{ 26 }$$

Using formula (25) to express $\dot{f}_{\mathbf{k}}^{(3)}(t)$ in the r.h.s. of relation (26) in terms of ${{\dot{\rho }}_{\mathbf{k}}}(t)$, we rewrite the Lagrangian L₀(t) in the new form containing $\dot{f}_{\mathbf{k}}^{(1,2)}(t)$, ${{\rho }_{\mathbf{k}}}(t)$ and ${{\dot{\rho }}_{\mathbf{k}}}(t)$:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{0}}(t) & = & \frac{1}{2}\underset{\mathbf{k}}{\mathop \sum }\,\left\{ m{{n}_{0}} \right. \\ {} & {} & \times \underset{\alpha =1,2}{\mathop \sum }\,\dot{f}_{\mathbf{k}}^{(\alpha )}(t)*\dot{f}_{\mathbf{k}}^{(\alpha )}(t) \\ {} & {} & +\frac{1}{{{e}^{2}}}\frac{m}{{{n}_{0}}}\frac{1}{{{k}^{2}}}\left[ {{{\dot{\rho }}}_{\mathbf{k}}}(t)*{{{\dot{\rho }}}_{\mathbf{k}}}(t) \right. \\ {} & {} & \left. \left. -\omega _{pl}^{2}{{\rho }_{\mathbf{k}}}(t)*{{\rho }_{\mathbf{k}}}(t) \right] \right\} \\ \end{array}\end{eqnarray} \tag{ 27 }$$

with

$$\begin{eqnarray}&&{{\omega }_{pl}}=\sqrt{\frac{4\pi {{e}^{2}}m}{{{n}_{0}}}}\end{eqnarray} \tag{ 28 }$$

called the plasma frequency.

From the Lagrangian equations

$$\begin{eqnarray*}&&\frac{d}{dt}\left( \frac{\partial {{L}_{0}}(t)}{\partial \dot{f}_{\mathbf{k}}^{(1,2)}(t)} \right)=\frac{\partial {{L}_{0}}(t)}{\partial f_{\mathbf{k}}^{(1,2)}(t)},\end{eqnarray*}$$

$$\begin{eqnarray*}&&\frac{d}{dt}\left( \frac{\partial {{L}_{0}}(t)}{\partial \dot{f}_{\mathbf{k}}^{(1,2)}(t)*} \right)=\frac{\partial {{L}_{0}}(t)}{\partial f_{\mathbf{k}}^{(1,2)}(t)*}\end{eqnarray*}$$

and expression (27) of L₀(t) it follows that the generalized velocities $\dot{f}_{\mathbf{k}}^{(1,2)}$ and $\dot{f}_{\mathbf{k}}^{{{(1,2)}^{*}}}$ must be time-independent. Because the addition of arbitrary constants to the Lagrangian L₀(t) does not affect the equation of motion of the system, the terms containing the constants $\dot{f}_{\mathbf{k}}^{(1,2)}$ and $\dot{f}_{\mathbf{k}}^{{{(1,2)}^{*}}}$ in the r.h.s. of formula (27) can be discarded. The absence of these constants in the Lagrangian has the following clear physical meaning: only the longitudinal displacements of electrons can cause the wave of propagating fluctuations of electron number density $n(\mathbf{r},t)$, i.e. of charge density, in the electron gas. Therefore instead of the expansion formula (20) we shall use following expression

$$\begin{eqnarray}&&\delta \mathbf{r}(\mathbf{r},t)=\frac{1}{\sqrt{V}}\underset{\mathbf{k}}{\mathop \sum }\,\frac{\mathbf{k}}{k}{{f}_{\mathbf{k}}}(t){{e}^{i\mathbf{kr}}}.\end{eqnarray} \tag{ 29 }$$

In this case the Lagrangian (26) becomes

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{0}}(t) & = & \frac{1}{2}\underset{\mathbf{k}}{\mathop \sum }\,\left[ m{{n}_{0}}{{{\dot{f}}}_{\mathbf{k}}}(t)*{{{\dot{f}}}_{\mathbf{k}}}(t) \right. \\ {} & {} & \left. -\frac{4\pi }{{{k}^{2}}}{{\rho }_{\mathbf{k}}}(t)*{{\rho }_{\mathbf{k}}}(t) \right] \\ \end{array}\end{eqnarray} \tag{ 30 }$$

and the subsidiary condition gives

$$\begin{eqnarray}\begin{matrix} {{{\dot{\rho }}}_{\mathbf{k}}}(t) & = & \frac{d{{\rho }_{\mathbf{k}}}(t)}{dt} \\ {} & = & ie{{n}_{0}}k{{{\dot{f}}}_{\mathbf{k}}}(t). \\ &&\end{matrix}\end{eqnarray} \tag{ 31 }$$

Using this condition we can rewrite the expression of the Lagrangian in the form containing only ${{\rho }_{\mathbf{k}}}(t),{{\rho }_{\mathbf{k}}}(t)*$ and ${{\dot{\rho }}_{\mathbf{k}}}(t),{{\dot{\rho }}_{\mathbf{k}}}(t)*$:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{0}}(t) & = & \frac{1}{2{{e}^{2}}}\frac{m}{{{n}_{0}}} \\ {} & {} & \times \underset{\mathbf{k}}{\mathop \sum }\,\frac{1}{{{k}^{2}}}\left[ {{{\dot{\rho }}}_{\mathbf{k}}}(t)*{{{\dot{\rho }}}_{\mathbf{k}}}(t)-\omega _{pl}^{2}{{\rho }_{\mathbf{k}}}(t)*{{\rho }_{\mathbf{k}}}(t) \right]. \\ \end{array}\end{eqnarray} \tag{ 32 }$$

Setting

$$\begin{eqnarray}\begin{matrix} \frac{1}{e}\sqrt{\frac{m}{{{n}_{0}}}}\frac{1}{k}{{\rho }_{\mathbf{k}}}(t) & = & q_{\mathbf{k}}^{(1)}(t)+iq_{\mathbf{k}}^{(2)}(t), \\ \frac{1}{e}\sqrt{\frac{m}{{{n}_{0}}}}\frac{1}{k}{{\rho }_{\mathbf{k}}}(t)* & = & q_{\mathbf{k}}^{(1)}(t)-iq_{\mathbf{k}}^{(2)}(t), \\ &&\end{matrix}\;\end{eqnarray} \tag{ 33 }$$

$q_{\mathbf{k}}^{(i)}(t)$ with i = 1, 2 being the real functions of t, we obtain the formula of the Lagrangian of a system of harmonic oscillators with the generalized real coordinates $q_{\mathbf{k}}^{(i)}(t)$ and the frequency ω_pl :

$$\begin{eqnarray}&&{{L}_{0}}(t)=\frac{1}{2}\underset{i=1,2}{\mathop \sum }\,\underset{\mathbf{k}}{\mathop \sum }\,\left[ \dot{q}_{\mathbf{k}}^{(i)}{{(t)}^{2}}-\omega _{pl}^{2}q_{\mathbf{k}}^{(i)}{{(t)}^{2}} \right].\end{eqnarray} \tag{ 34 }$$

In order to quantize the system with the Lagrangian (34) we must calculate the generalized momenta $p_{\mathbf{k}}^{(i)}(t)$ as well as the Hamiltonian H. We obtain

$$\begin{eqnarray*}&&p_{\mathbf{k}}^{(i)}(t)=\frac{\partial {{L}_{0}}(t)}{\partial \dot{q}_{\mathbf{k}}^{(i)}(t)}=\dot{q}_{\mathbf{k}}^{(i)}(t)\end{eqnarray*}$$

and

$$\begin{eqnarray}&&H=\frac{1}{2}\underset{i=1,2}{\mathop \sum }\,\underset{\mathbf{k}}{\mathop \sum }\,\left[ p_{\mathbf{k}}^{(i)}{{(t)}^{2}}+\omega _{pl}^{2}q_{\mathbf{k}}^{(i)}{{(t)}^{2}} \right].\end{eqnarray} \tag{ 35 }$$

The quantization procedure consists of the replacement of dynamical canonical variables $q_{\mathbf{k}}^{(i)}(t)$ and $p_{\mathbf{k}}^{(i)}(t)$ by the corresponding hermitian operators ${\mathop{\hat{q}}}\,_{\mathbf{k}}^{(i)}$ and ${\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)}$, i = 1, 2, satisfying the following canonical commutation relations

$$\begin{eqnarray}&&\begin{matrix} \left[ {\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)},{\mathop{\hat{p}}}\,_{\mathbf{l}}^{(j)} \right]-\left[ {\mathop{\hat{q}}}\,_{\mathbf{k}}^{(i)},{\mathop{\hat{q}}}\,_{\mathbf{l}}^{(j)} \right]=0, \\ &&\quad \left[ {\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)},{\mathop{\hat{q}}}\,_{\mathbf{l}}^{(j)} \right]=-i\hbar {{\delta }_{ij}}{{\delta }_{\mathbf{kl}}}. \\ &&\end{matrix}\end{eqnarray} \tag{ 36 }$$

Hamiltonian of the quantized system is the operator

$$\begin{eqnarray}&&H=\frac{1}{2}\underset{i=1,2}{\mathop \sum }\,\underset{\mathbf{k}}{\mathop \sum }\,\left[ {\mathop{\hat{p}}}\,_{\mathbf{k}}^{{{(i)}^{2}}}+\omega _{pl}^{2}{\mathop{\hat{q}}}\,_{\mathbf{k}}^{{{(i)}^{2}}} \right].\end{eqnarray} \tag{ 37 }$$

Instead of two hermitian operators ${\mathop{\hat{q}}}\,_{\mathbf{k}}^{(i)}$ and ${\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)}$ for each set of two indices i and k we introduce two operators

$$\begin{eqnarray}\begin{array}{*{35}{l}} {\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)} & = & \frac{1}{\sqrt{2\hbar {{\omega }_{pl}}}}{\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)}-i\sqrt{\frac{{{\omega }_{pl}}}{2\hbar }}{\mathop{\hat{q}}}\,_{\mathbf{k}}^{(i)}, \\ {\mathop{\hat{a}}}\,_{\mathbf{k}}^{{{(i)}^{+}}} & = & \frac{1}{\sqrt{2\hbar {{\omega }_{pl}}}}{\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)}+i\sqrt{\frac{{{\omega }_{pl}}}{2\hbar }}{\mathop{\hat{q}}}\,_{\mathbf{k}}^{(i)}. \\ \end{array}\end{eqnarray} \tag{ 38 }$$

They satisfy commutation relations

$$\begin{eqnarray}\begin{array}{*{35}{l}} \left[ {\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)},{\mathop{\hat{a}}}\,_{\mathbf{l}}^{(j)} \right] & = & \left[ {\mathop{\hat{a}}}\,_{\mathbf{k}}^{{{(i)}^{+}}},{\mathop{\hat{a}}}\,_{\mathbf{l}}^{{{(j)}^{+}}} \right]=0, \\ \left[ {\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)},{\mathop{\hat{a}}}\,_{\mathbf{l}}^{{{(j)}^{+}}} \right] & = & {{\delta }_{ij}}{{\delta }_{\mathbf{kl}}}. \\ \end{array}\end{eqnarray} \tag{ 39 }$$

Inversely we can express ${\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)}$ and ${\mathop{\hat{q}}}\,_{\mathbf{k}}^{(i)}$ in terms of ${\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)}$ and ${\mathop{\hat{a}}}\,_{\mathbf{k}}^{{{(i)}^{+}}}$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {\mathop{\hat{p}}}\,_{\mathbf{k}}^{(i)} & = & \sqrt{\frac{\hbar {{\omega }_{pl}}}{2}}\left[ {\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)}+{\mathop{\hat{a}}}\,_{\mathbf{k}}^{{{(i)}^{+}}} \right], \\ {\mathop{\hat{q}}}\,_{\mathbf{k}}^{(i)} & = & i\sqrt{\frac{\hbar }{2{{\omega }_{pl}}}}\left[ {\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)}-{\mathop{\hat{a}}}\,_{\mathbf{k}}^{{{(i)}^{+}}} \right], \\ \end{array}\end{eqnarray} \tag{ 40 }$$

and obtain the Hamiltonian operator in the form

$$\begin{eqnarray}&&H=\hbar {{\omega }_{pl}}\underset{i=1,2}{\mathop \sum }\,\underset{\mathbf{k}}{\mathop \sum }\,\left( {\mathop{\hat{a}}}\,_{\mathbf{k}}^{{{(i)}^{+}}}{\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)}+\frac{1}{2} \right).\end{eqnarray} \tag{ 41 }$$

Operators ${\mathop{\hat{a}}}\,_{\mathbf{k}}^{(i)}$ and ${\mathop{\hat{a}}}\,_{\mathbf{k}}^{{{(i)}^{+}}}$ satisfying the commutation relations (39) play the role of the destruction operator and creation operator, respectively, for a quasiparticle of the type i and with the momentum $\hbar \mathbf{k}$. The quasiparticles of this kind are the quanta of the collective excitations, i.e. of the fluctuations of the electron number density $n(\mathbf{r},t)$ in the electron gas. They are called plasmons.

Since the electron gas behaves like a system of harmonic oscillators in the above presented lowest order approximation with respect to the small oscillating displacements of electrons, this approximation in often called the harmonic approximation. Note that in the harmonic approximation the plasma energy does not depend on its momentum and equals $\hbar {{\omega }_{pl}}$. The plasmons in the system with Lagrangian L₀(t) do not participate in any interaction process and are called the free plasmons.

Now we return to the total Lagrangian (9) and the subsidiary condition (16) and study the fluctuations of the charge density $\rho (\mathbf{r},t)$ beyond the harmonic approximation. Substituting the Fourier series (19) and (29) for $\rho (\mathbf{r},t)$ and $\delta \mathbf{r}(\mathbf{r},t)$ into both sides of the continuity equation (16), we obtain following system of equations between their Fourier components

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\dot{\rho }}}_{\mathbf{k}}}(t) & = & ie{{n}_{0}}k{{{\dot{f}}}_{\mathbf{k}}}(t)-\frac{i}{\sqrt{V}} \\ {} & {} & \times \underset{\mathbf{l}}{\mathop \sum }\,{{\rho }_{\mathbf{k}-\mathbf{l}}}(t)\frac{\mathbf{kl}}{l}{{{\dot{f}}}_{\mathbf{l}}}(t). \\ \end{array}\end{eqnarray} \tag{ 42 }$$

In order to express ${{\dot{f}}_{\mathbf{k}}}(t)$ in term of ${{\rho }_{\mathbf{k}}}(t)$ and ${{\dot{\rho }}_{\mathbf{k}}}(t)$ we rewrite equation (42) in another form

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\dot{f}}}_{\mathbf{k}}}(t) & = & -\frac{i}{e{{n}_{0}}}\frac{1}{k}{{{\dot{\rho }}}_{\mathbf{k}}}(t)+\frac{1}{e{{n}_{0}}}\frac{1}{\sqrt{V}} \\ {} & {} & \times \underset{\mathbf{l}}{\mathop \sum }\,\frac{1}{k}{{\rho }_{k-l}}(t)\mathbf{kl}\frac{1}{l}{{{\dot{f}}}_{\mathbf{l}}}(t) \\ \end{array}\end{eqnarray} \tag{ 43 }$$

and solve the system of equations (43) with respect to ${{\dot{f}}_{\mathbf{k}}}(t)$ by means of the iteration procedure. Then we obtain the expression of ${{\dot{f}}_{\mathbf{k}}}(t)$ in the form of a functional power series of ${{\rho }_{\mathbf{k}}}(t)$ and ${{\dot{\rho }}_{\mathbf{k}}}(t)$. Up to the third order we have

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\dot{f}}}_{\mathbf{k}}}(t) & = & -\frac{i}{e{{n}_{0}}}\frac{1}{k}{{{\dot{\rho }}}_{\mathbf{k}}}(t)-\frac{i}{{{\left( e{{n}_{0}} \right)}^{2}}}\frac{1}{k}\frac{1}{\sqrt{V}} \\ {} & {} & \times \underset{\mathbf{l}}{\mathop \sum }\,\frac{1}{k}{{\rho }_{\mathbf{k}-\mathbf{l}}}(t)\frac{\mathbf{kl}}{{{l}^{2}}}{{{\dot{\rho }}}_{\mathbf{l}}}(t)-\frac{i}{{{\left( e{{n}_{0}} \right)}^{3}}}\frac{1}{k}\frac{1}{\sqrt{V}} \\ {} & {} & \times \underset{\mathbf{l}}{\mathop \sum }\,\frac{1}{k}{{\rho }_{\mathbf{k}-\mathbf{l}}}(t)\frac{\mathbf{kl}}{{{l}^{2}}}\frac{1}{\sqrt{V}} \\ {} & {} & \times \underset{\mathbf{l}\prime }{\mathop \sum }\,{{\rho }_{\mathbf{l}-{{\mathbf{l}}^{\prime }}}}(t)\frac{\mathbf{ll}\prime }{{{l}^{\prime 2}}}\frac{1}{\sqrt{V}}{{{\dot{\rho }}}_{{{\mathbf{l}}^{\prime }}}}(t) \\ {} & {} & +\text{higher}\;\text{order}\;\text{terms}. \\ \end{array}\end{eqnarray} \tag{ 44 }$$

Consider the total Lagrangian (9). Using the decompositions (19) and (29) for $\rho (\mathbf{r},t)$ and $\delta \mathbf{r}(\mathbf{r},t)$, we express L(t) in terms of their Fourier components ${{\rho }_{\mathbf{k}}}(t)$ and ${{f}_{\mathbf{k}}}(t)$ as follows

$$\begin{eqnarray}\begin{array}{*{35}{l}} L(t) & = & \frac{1}{2}\underset{\mathbf{k}}{\mathop \sum }\,\left[ m{{n}_{0}}{{{\dot{f}}}_{\mathbf{k}}}(t)*{{{\dot{f}}}_{\mathbf{k}}}(t) \right. \\ {} & {} & \left. -\frac{4\pi }{{{k}^{2}}}{{{\dot{\rho }}}_{\mathbf{k}}}(t)*{{{\dot{\rho }}}_{\mathbf{k}}}(t) \right] \\ {} & {} & -\frac{1}{4}\frac{m}{e}\frac{1}{\sqrt{V}} \\ {} & {} & \times \underset{\mathbf{k}}{\mathop \sum }\,\underset{\mathbf{l}}{\mathop \sum }\,\frac{\mathbf{kl}}{kl}\left[ {{\rho }_{\mathbf{k}+\mathbf{l}}}(t)*{{{\dot{f}}}_{\mathbf{k}}}(t){{{\dot{f}}}_{\mathbf{l}}}(t) \right. \\ {} & {} & \left. +{{{\dot{f}}}_{\mathbf{k}}}(t)*{{{\dot{f}}}_{\mathbf{l}}}(t)*{{\rho }_{\mathbf{k}+\mathbf{l}}}(t) \right]. \\ \end{array}\end{eqnarray} \tag{ 45 }$$

Substituting the expansion (44) of ${{\dot{f}}_{\mathbf{k}}}(t)$ into the r.h.s. of formula (45) and dividing the derived expression of the total Lagrangian L(t) into two parts

$$\begin{eqnarray}&&L(t)={{L}_{0}}(t)+{{L}_{\operatorname{int}}}(t),\end{eqnarray} \tag{ 46 }$$

where L₀(t) is the expression (32) of the Lagrangian in the harmonic approximation and L_int(t) is called the interaction Lagrangian. In the fourth order with respect to the functions ${{\rho }_{\mathbf{k}}}(t)$, ${{\rho }_{\mathbf{k}}}(t)*$ and ${{\dot{\rho }}_{\mathbf{k}}}(t)$, ${{\dot{\rho }}_{\mathbf{k}}}(t)*$ we obtain

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{\operatorname{int}}}(t) & = & \frac{1}{2}\frac{m{{n}_{0}}}{{{(e{{n}_{0}})}^{3}}} \\ {} & {} & \left[ \frac{1}{\sqrt{V}}\underset{\mathbf{k}}{\mathop \sum }\,\underset{\mathbf{l}}{\mathop \sum }\,\frac{\mathbf{kl}}{{{k}^{2}}{{l}^{2}}}{{{\dot{\rho }}}_{\mathbf{k}}}(t)*{{\rho }_{\mathbf{k}-\mathbf{l}}}(t){{{\dot{\rho }}}_{\mathbf{l}}}(t) \right. \\ {} & {} & +\frac{1}{e{{n}_{0}}}\frac{1}{V}\underset{\mathbf{k}}{\mathop \sum }\,\underset{\mathbf{l}}{\mathop \sum }\,\underset{\mathbf{l}\prime }{\mathop \sum }\,\frac{(\mathbf{kl})(\mathbf{ll}\prime )}{{{k}^{2}}{{l}^{2}}{{l}^{\prime 2}}} \\ {} & {} & \left. \times {{\rho }_{\mathbf{l}-{{\mathbf{l}}^{\prime }}}}(t)*{{{\dot{\rho }}}_{{{\mathbf{l}}^{\prime }}}}(t)*{{{\dot{\rho }}}_{\mathbf{k}}}(t){{\rho }_{\mathbf{l}-\mathbf{k}}}(t) \right]. \\ \end{array}\end{eqnarray} \tag{ 47 }$$

Note that due to the property (24) of the functions ${{\rho }_{\mathbf{k}}}(t)$ and ${{\dot{\rho }}_{\mathbf{k}}}(t)$, interaction Lagrangian (47) is a real function of t.

In this work we have elaborated the quantum theory of collective excitations of the isotropic and homogeneous electron gas in the framework of the canonical formalism. The quanta of these collective excitations are the quasiparticles of some special type called plasmons. In this formalism there naturally appear the destruction and creation operators of plasmons during the standard quantization procedure. By applying this method we can establish the Lagrangian of the interaction between plasmons and other particles or quasiparticles in condensed matters from the first principle.

In the harmonic approximation the Lagrangian of the system of electrons in the electron gas can be represented as that of the free plasmons with a common momentum-independent energy $\hbar {{\omega }_{pl}}$. In higher order approximation with respect to the collective excitations of the electron gas, the total Lagrangian consists of two parts, one is that of free plasmons, another is the interaction Lagrangian describing the plasmon–plasmon interaction. The dependence of the energy of a plasmon on its momentum is the consequence of the renormalization of the one-plasmon state vector (or two-point Green function of the plasmon) due to the plasmon–plasmon interaction.

As a simple example of the interaction of plasmons with other particles let us consider the photon–plasmon interaction. The photons are described by a quantized elecromagnetic field with the vector potential ${\mathop{\hat{\mathbf{A}}}}\,(\mathbf{r},t)$ and the scalar potential ${\hat{\mathop{V}}}\,(\mathbf{r},t)$. The Lagrangian of the interaction between photons and all electrons of the electron gas is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{\gamma -e}}(t) & = & m\frac{e}{c}\times \mathop{\int }^{}d\mathbf{r}n(\mathbf{r},t)\delta \mathbf{\dot{r}}(\mathbf{r},t){\mathop{\hat{\mathbf{A}}}}\,(\mathbf{r},t) \\ {} & {} & +\frac{1}{2}m\frac{{{e}^{2}}}{{{c}^{2}}}\mathop{\int }^{}d\mathbf{r}n(\mathbf{r},t)\mathbf{A}{{(\mathbf{r},t)}^{2}} \\ {} & {} & +e\mathop{\int }^{}d\mathbf{r}\left[ n(\mathbf{r},t)-{{n}_{0}} \right]{\hat{\mathop{V}}}\,(\mathbf{r},t). \\ \end{array}\end{eqnarray} \tag{ 48 }$$

where c is light velocity in the vacuum.

In terms of the displacements $\delta \mathbf{r}(\mathbf{r},t)$ and the charge density $\rho (\mathbf{r},t)$ we have

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{L}_{\gamma -e}}(t) & = & m{{n}_{0}}\frac{e}{c}\mathop{\int }^{}d\mathbf{r}\delta \mathbf{\dot{r}}(\mathbf{r},t){\mathop{\hat{\mathbf{A}}}}\,(\mathbf{r},t) \\ {} & {} & -\frac{m}{c}\mathop{\int }^{}d\mathbf{r}\;\rho (\mathbf{r},t)\delta \mathbf{\dot{r}}(\mathbf{r},t){\mathop{\hat{\mathbf{A}}}}\,(\mathbf{r},t) \\ {} & {} & +\frac{1}{2}m{{n}_{0}}\frac{{{e}^{2}}}{{{c}^{2}}}\mathop{\int }^{}d\mathbf{r}\mathbf{A}{{(\mathbf{r},t)}^{2}} \\ {} & {} & -\frac{1}{2}m\frac{e}{{{c}^{2}}}\mathop{\int }^{}d\mathbf{r}\;\rho (\mathbf{r},t){\mathop{\hat{\mathbf{A}}}}\,{{(\mathbf{r},t)}^{2}} \\ {} & {} & -\mathop{\int }^{}d\mathbf{r}\;\rho (\mathbf{r},t){\hat{\mathop{V}}}\,(\mathbf{r},t). \\ \end{array}\end{eqnarray} \tag{ 49 }$$

By canonically quantizing the dynamical variable $\delta \mathbf{r}(\mathbf{r},t)$ and $\rho (\mathbf{r},t)$, we obtain the plasmon–photon interaction Lagrangian.

The authors would like to express their gratitude to the Vietnam Academy of Science and technology for the support to this work.