Quantum dynamics of plasmons in nanomaterials

The origin of the appearance of the resonance in collective oscillations of a system of itinerant electrons with their Coulomb interaction is the existence of a quasiparticle called plasmon. In recent years there has been significant progress in the study of plasmon resonances in the physical and chemical processes based on the light–matter interaction in nanomaterials. Many interesting and valuable results were obtained in the studies of plasmonic–molecular resonance coupling [1–13], which is also called exciton–plasmon coupling [9,10] or coherent coupling of molecular excitons to electronic polarizations of metal nanoparticles [2], and takes place also in the coupled quantum dot-metal nanoparticle systems [14–16]. Plasmon-enhanced fluorescence [17–23] and plasmon resonance energy transfer [24, 25] were applied for improving the fluorescence labeling and sensing of biomolecules. Numerous applications of plasmon resonance in science were presented in a recent comprehensive review [26] (with more than a thousand references). In the theoretical studies of above-mentioned interaction processes of plasmons, two methods were used: either one introduced the plasmon destruction and creation operators and assume some simple expression of the effective Hamiltonian containing these operators [10], or one wrote the effective Hamiltonian in the form containing the polarization of the metallic nanoparticles and used the corresponding term in the Hamiltonian at the plasmon resonance to describe the interaction of the plasmon with other quantum particles [14–16].

In this work, we present a new approach to the study of quantum dynamics of plasmons in metallic nanostructures. It is based on the functional integral technique [27–32]. Denote by ψ (x,t) the Grassmann (anticommuting) space-time dependent variable playing the role of the wave function of electron and $\bar \psi ({\bf x},t)$ its hermitian conjugate, where x is the radius vector determining the position of the relevant point in the nanostructure. For simplifying formulae, we omit the spin index of the electron wave function and therefore consider electron as a spinless fermion (without the loss of the generality of the obtained results). For elaborating mathematical tools and presenting basics of the theory in section 2, we consider the system of non-interacting itinerant electrons in a confining potential. We define the contraction of two field operators and the Green functions in terms of the functional integrals, and re-derive the well-known Wick theorem in quantum field theory. In section 3 we start from the generating functional of the system of interacting itinerant electrons in a confining potential, the electron–electron interaction being the Coulomb repulsion between electrons. By means of the Hubbard–Stratonovich transformation [33, 34] we rewrite the expression containing the electron–electron interaction Hamiltonian in the generating functional in the new form with the presence of a scalar field φ (x,t), which plays the role of the order parameter characterizing the collective excitation of the electron gas and is related to the quantized scalar field whose quanta are the plasmons. The effective action and the dynamical equations for this scalar field are explicitly presented. As a test of the validity of the derived formulae and expressions, in section 4 we apply them to a special case—that of the electron gas in the three-dimentional free space. We check that in this special case they lead to the same physical consequences as those of the conventional well-known method. In section 5 we establish the quantum field with plasmons as its quanta—plasmonic quantum field. Conclusion and discussion are presented in section 6.

Consider first the electron gas in a three-dimensional nanostructure with some confining potential U(x), the electron–electron Coulomb repulsion being neglected. For simplifying formulae, we introduce four-dimensional space-time with the coordinate x =  (x,t) = (x,x₀), where t = x₀ is the time variable, and denote

$$\begin{eqnarray*} \psi (x) = \psi ({\bf x},t), & \bar \psi (x) = \bar \psi ({\bf x},t), & \\ \varphi (x) = \varphi ({\bf x},t), & \int {{\rm{d}}x} = \int\limits_V {{\rm{d}}{\bf x}} \int {{\rm{d}}t}, \end{eqnarray*}$$

where V is the volume of the nanostructure. The boundary surface of V is denoted by S =  ∂V. The fermionic functional integration variables ψ (x) and $\bar \psi (x)$ are the Grassman variables. They anticommute one with another

$$\begin{eqnarray} &&\psi (x)\psi (y) + \psi (y)\psi (x) = 0, \nonumber\\ &&\psi (x)\bar \psi (y) + \bar \psi (y)\psi (x) = 0, \\ &&\bar \psi (x)\bar \psi (y) + \bar \psi (y)\bar \psi (x) = 0,\nonumber \end{eqnarray} \tag{ 1 }$$

while bosonic functional integration variables φ (x) are commuting

$$\begin{equation} \varphi (x)\varphi (y) = \varphi (y)\varphi (x).\end{equation} \tag{ 2 }$$

In the functional integration formalism one often uses Grassmann parameter η (x) and its conjugate $\bar \eta (x)$. They anticommute one with another

$$\begin{eqnarray} &&\eta (x)\eta (y) + \eta (y)\eta (x) = 0, \nonumber\\ &&\eta (x)\bar \eta (y) + \bar \eta (y)\eta (x) = 0, \\ &&\bar \eta (x)\bar \eta (y) + \bar \eta (y)\bar \eta (x) = 0 \nonumber \end{eqnarray} \tag{ 3 }$$

and also anticommute with fermionic fields ψ (x) and $\bar \psi (x)$:

$$\begin{equation}\begin{array}{l} \eta (x)\psi (y) + \psi (y)\eta (x) = 0, \\ \eta (x)\bar \psi (y) + \bar \psi (y)\eta (x) = 0, \\ \bar \eta (x)\psi (y) + \psi (y)\bar \eta (x) = 0, \\ \bar \eta (x)\bar \psi (y) + \bar \psi (y)\bar \eta (x) = 0. \end{array}\end{equation} \tag{ 4 }$$

The functional derivation operators $\frac{\delta }{{\delta \eta (x)}}$ and $\frac{\delta }{{\delta \bar \eta (x)}}$ anticommute one with another

$$\begin{eqnarray} \frac{\delta }{{\delta \eta (x)}}\frac{\delta }{{\delta \eta (y)}} + \frac{\delta }{{\delta \eta (y)}}\frac{\delta }{{\delta \eta (x)}} &=& 0, \nonumber\\ \frac{\delta }{{\delta \eta (x)}}\frac{\delta }{{\delta \bar \eta (y)}} + \frac{\delta }{{\delta \bar \eta (y)}}\frac{\delta }{{\delta \eta (x)}} &= &0, \\ \frac{\delta }{{\delta \bar \eta (x)}}\frac{\delta }{{\delta \bar \eta (y)}} + \frac{\delta }{{\delta \bar \eta (y)}}\frac{\delta }{{\delta \bar \eta (x)}} &= &0 \nonumber\end{eqnarray} \tag{ 5 }$$

and also anticommute with fermionic fields ψ (x) and $\bar \psi (x)$:

$$\begin{equation}\begin{array}{l} \displaystyle\frac{\delta }{{\delta \eta (x)}}\psi (y) + \psi (y)\frac{\delta }{{\delta \eta (x)}} = 0, \\ \displaystyle \frac{\delta }{{\delta \eta (x)}}\bar \psi (y) + \bar \psi (y)\frac{\delta }{{\delta \eta (x)}} = 0, \\ \displaystyle\frac{\delta }{{\delta \bar \eta (x)}}\psi (y) + \psi (y)\frac{\delta }{{\delta \bar \eta (x)}} = 0, \\ \displaystyle\frac{\delta }{{\delta \bar \eta (x)}}\bar \psi (y) + \bar \psi (y)\frac{\delta }{{\delta \bar \eta (x)}} = 0. \end{array} \end{equation} \tag{ 6 }$$

From above-mentioned anticommutativity properties of fermionic fields ψ (x), $\bar \psi (x)$, Grassman variables η (x), $\bar \eta (x)$ and functional derivation operators $\frac{\delta }{{\delta \eta (x)}}$, $\frac{\delta }{{\delta \bar \eta (x)}}$, it is straightforward to derive following formulae which will be often used in this work:

$$\begin{eqnarray} &&\frac{\delta }{{\delta \bar \eta (y)}}\,{\rm exp}\!\left\{\! {{\rm{i}}\int \!\!{{\rm{d}}x\bar \eta (x)\psi (x)} } \!\right\} = {\rm{i}}\psi (y)\,{\rm exp}\!\left\{\! {{\rm{i}}\int {{\rm{d}}x\bar \eta (x)\psi (x)} } \!\right\}\!\!,\nonumber\\ \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} &&\frac{\delta }{{\delta \eta (z)}}\,{\rm exp}\!\left\{\! {{\rm{i}}\int \!\!{{\rm{d}}x\bar \psi (x)\eta (x)} } \!\right\} = - {\rm{i}}\bar \psi (z)\,{\rm exp}\!\left\{\! {{\rm{i}}\int \!\!{{\rm{d}}x\bar \psi (x)\eta (x)} } \!\right\}\!\!,\nonumber\\ \end{eqnarray}\noindent \tag{ 8 }$$

$$\begin{eqnarray} &&\frac{{\delta ^2 }}{{\delta \eta (z)\delta \bar \eta (y)}}\,{\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x\left[ {\bar \psi (x)\eta (x) + \bar \eta (x)\psi (x)} \right]} } \right\} \nonumber\\ &&= \bar \psi (z)\psi (y)\,{\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x\left[ {\bar \psi (x)\eta (x) + \bar \eta (x)\psi (x)} \right]} } \right\}, \end{eqnarray}\noindent \tag{ 9 }$$

$$\begin{eqnarray} &&\displaystyle\frac{{\delta ^2 }}{{\delta \bar \eta (z)\delta \bar \eta (y)}}\,{\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x\bar \eta (x)\psi (x)\,} } \right\} \nonumber\\ &&= - \psi (z)\psi (y)\,{\rm exp}\,\left\{ {{\rm{i}}\bar \eta (x)\psi (x)} \right\}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &&\frac{{\delta ^2 }}{{\delta \eta (z)\delta \eta (y)}}\,{\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x\bar \psi (x)\eta (x)} } \right\} \nonumber\\ && = - \bar \psi (z)\bar \psi (y)\,{\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x\bar \psi (x)\eta (x)} } \right\}. \end{eqnarray} \tag{ 11 }$$

We denote $H\left( { - {\rm i}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)$ the Hamiltonian of electron with the mass m in a confining potential U(x)

$$\begin{equation}H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right) = \frac{1}{{2m}}\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}}} \right)^2 + U({\bf x}),\end{equation} \tag{ 12 }$$

introduce the generating functional of the electron system

$$\begin{eqnarray} &&Z_0^\psi [ {\eta,\bar \eta } ] = \int \!\!{[ {D\psi } ]} [ {D\bar \psi } ]{\rm exp}\!\left\{\! {{\rm{i}}\int \!\!{{\rm{d}}x[ {\bar \psi (x)\eta (x) {+} \bar \eta (x)\psi (x)} ]} } \!\right\} \nonumber\\ && \times {\rm exp}\,\left\{ {{\rm{i}}\int {\rm{d}} x\bar \psi (x)\left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm i}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]\psi (x)} \right\}, \end{eqnarray} \tag{ 13 }$$

and define the contraction between two fermionic fields at two different points x,y of the nanostructure, and at two different time momenta x₀, y₀ as follows:

$$\begin{equation}\begin{array}{l} \displaystyle\psi (y)_c \bar \psi (z)_c = \frac{1}{{Z_0^\psi }}\int {[ {D\psi } ]} [ {D\bar \psi } ]\psi (y)\bar \psi (z) \\ \displaystyle\times{\rm exp}\,\left\{ {{\rm{i}}\int {\rm{d}} x\bar \psi (x)\left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]\psi (x)} \right\}, \end{array} \end{equation} \tag{ 14 }$$

$$\begin{equation}\begin{array}{l} \displaystyle \psi (y)_c \psi (z)_c = \frac{1}{{Z_0^\psi }}\int {[ {D\psi } ]} [ {D\bar \psi } ]\psi (y)\psi (z) \\ \displaystyle\times {\rm exp}\,\left\{ {{\rm{i}}\int {\rm{d}} x\bar \psi (x)\left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]\psi (x)} \right\}, \end{array} \end{equation} \tag{ 15 }$$

$$\begin{equation}\begin{array}{l} \displaystyle\bar \psi (y)_c \bar \psi (z)_c = \frac{1}{{Z_0^\psi }}\int {[ {D\psi } ]} [ {D\bar \psi } ]\bar \psi (y)\bar \psi (z) \\ \displaystyle\times {\rm exp}\,\left\{ {{\rm{i}}\int {\rm{d}} x\bar \psi (x)\left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]\psi (x)} \right\}, \end{array} \end{equation} \tag{ 16 }$$

where

$$\begin{eqnarray} Z_0^\psi &=& Z_0^\psi [ {0,0} ] = \int {[ {D\psi } ]} [ {D\bar \psi } ] \, {\rm exp}\,\left\{ {\rm{i}}\int {\rm{d}} x\bar \psi (x)\right.\nonumber\\ &&\left.\times \left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]\psi (x) \right\}. \end{eqnarray} \tag{ 17 }$$

Using formulae (9)–(11), from the definitions (14)–(16) we obtain

$$\begin{equation}\psi (y)_c \bar \psi (z)_c = - \frac{1}{{Z_0^\psi }}\left. {\frac{{\delta ^2 Z_0^\psi \left[ {\eta,\bar \eta } \right]}}{{\delta \eta (z)\delta \bar \eta (y)}}} \right|_{\eta = \bar \eta = 0},\end{equation} \tag{ 18 }$$

$$\begin{equation}\psi (y)_c \psi (z)_c = \frac{1}{{Z_0^\psi }}\left. {\frac{{\delta ^2 Z_0^\psi \left[ {\eta,\bar \eta } \right]}}{{\delta \bar \eta (z)\delta \bar \eta (y)}}} \right|_{\eta = \bar \eta = 0},\end{equation} \tag{ 19 }$$

$$\begin{equation}\bar \psi _c (y)\bar \psi _c (z) = \frac{1}{{Z_0^\psi }}\left. {\frac{{\delta ^2 Z_0^\psi \left[ {\eta,\bar \eta } \right]}}{{\delta \eta (z)\delta \eta (y)}}} \right|_{\eta = \bar \eta = 0}.\end{equation} \tag{ 20 }$$

Now we establish the explicit form of the generating functional (13). Denote $\{ {u_\alpha ({\bf x})} \}$ the complete set of eigenstates of Hamiltonian (12) with corresponding eigenvalues E_α

$$\begin{equation}Hu_\alpha ({\bf x}) = E_\alpha u_\alpha ({\bf x}),\end{equation} \tag{ 21 }$$

and S(x,y) =S(x,y;x₀ − y₀ ) the solution of the inhomogeneous differential equation

$$\begin{eqnarray} \left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left({ - \frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]S(x,y) &=& \delta (x - y)\nonumber\\ & =& \delta ({\bf x} {-} {\bf y})\delta (x_0 {-} y_0 ).\quad\quad\end{eqnarray} \tag{ 22 }$$

S(x,y) has the general form

$$\begin{eqnarray} S\left( {{\bf x},{\bf y};x_0 - y_0 } \right) &=& \frac{1}{{2\pi }}\int {\rm{d}} \omega \,{\rm{e}}^{ - {\rm{i}}\omega (x_0 - y_0 )} \sum\limits_\alpha {\frac{{u_\alpha ({\bf x})u_\alpha ({\bf y})^ * }}{{\omega - E_\alpha + {\rm i}o}}} \nonumber\\ &&+ {\rm{i}}\sum\limits_\alpha {C_\alpha \,{\rm{e}}^{ - {\rm{i}}E_\alpha (x_0 - y_0 )} u_\alpha ({\bf x})u_\alpha ({\bf y})^ * } \nonumber\\ &=& \frac{1}{{2\pi }}\int {\rm{d}} \omega \,{\rm{e}}^{ - {\rm{i}}\omega (x_0 - y_0 )} \sum\limits_\alpha u_\alpha ({\bf x})u_\alpha ({\bf y})^ * \nonumber\\ &&\times\left\{ {\frac{{1 - C_\alpha }}{{\omega - E_\alpha + {\rm i}o}} + \frac{{C_\alpha }}{{\omega - E_\alpha - {\rm i}o}}} \right\}, \end{eqnarray} \tag{ 23 }$$

where the constants C_α are related with the physical characterizations of the electron gas. Their physical meaning will be determined later.

Let us now perform the shift of the functional integration variables

$$\begin{equation} \begin{array}{l} \displaystyle\psi (x) \to \psi (x) + \int {{\rm{d}}y\,S(x,y)\eta (y)}, \\ \displaystyle\bar \psi (x) \to \bar \psi (x) + \int {{\rm{d}}y\bar \eta (y)S(y,x)} \end{array} \end{equation} \tag{ 24 }$$

in the rhs of formula (17). The functional integral in this formula must be invariant under the shift (24) of the functional integration variables. After lengthy but standard transformations and calculations related to the rhs of formula (17) we obtain following relation:

$$\begin{eqnarray} Z_0^\psi &=& \int {[ {D\psi } ]} \left[ {D\bar \psi } \right]{\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x\left[ {\bar \psi (x)\eta (x) + \bar \eta (x)\psi (x)} \right]} } \right\} \nonumber\\ &&\times \,{\rm exp}\,\left\{ {{\rm{i}}\int {\rm{d}} x\bar \psi (x)\left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm{i}}\frac{\partial }{{\partial \vec x}},\vec x} \right)} \right]\psi (x)} \right\} \nonumber\\ &&\times{\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x} \int {{\rm{d}}y\bar \eta (x)S(x,y)\eta (y)} } \right\}\; \cdot \end{eqnarray} \tag{ 25 }$$

Comparing the rhs of relation (25) with the definition (13) of $Z_0^\psi \left[ {\eta,\bar \eta } \right]$, we derive the explicit expression of this generating functional

$$\begin{equation}Z_0^\psi \left[ {\eta,\bar \eta } \right] = Z_0^\psi \,{\rm exp}\,\left\{ { - {\rm{i}}\int {{\rm{d}}x\int {dy} \bar \eta (x)S(x,y)\eta (y)} } \right\}.\end{equation} \tag{ 26 }$$

From expression (26) and formulae (18)–(20) it follows that

$$\begin{equation}\psi (y)_c \bar \psi (z)_c = - \bar \psi (z)_c \psi (y)_c = {\rm{i}}S(y,z),\end{equation} \tag{ 27 }$$

$$\begin{equation}\psi (y)_c \psi (z)_c = \bar \psi (y)_c \bar \psi (z)_c = 0.\end{equation} \tag{ 28 }$$

Let us introduce the notion of the average of an expression of fermionic fields ψ (x) and $\bar \psi (x)$ over all configurations of these fields. We consider some expression containing these fields at different points x₁, x₂,...,x_n and y₁, y₂,...,y_m of the four-dimensional space-time $A[ {\psi (x_1 ),\psi (x_2 ), \ldots,\psi (x_n );} \bar \psi (y_1 ),\bar \psi (y_2 ), \ldots,{\bar \psi (y_m )} ]$ and define its average over all configurations of the fermionic fields of electrons as follows:

$$\begin{eqnarray} &&\left\langle\kern-1pt {A\left[\kern-1pt{\psi (x_1 ),\psi (x_2 ), \ldots,\psi (x_n );\bar \psi (y_1 ),\bar \psi (y_2 ), \ldots,\bar \psi (y_m )}\kern-1pt\right]} \right\rangle _{\psi,\bar \psi } \nonumber\\ && = \frac{1}{{Z_0^\psi }}\int {[ {D\psi } ][ {D\bar \psi } ]} \;A[ \psi (x_1 ),\psi (x_2 ), \ldots,\psi (x_n );\bar \psi (y_1 ),\nonumber\\ &&\quad \bar \psi (y_2 ), \ldots,\bar \psi (y_m ) ] {\rm exp}\,\left\{ {\rm{i}}\int {\rm{d}} x\bar \psi (x)\left[ {\rm{i}}\frac{\partial }{{\partial x_0 }} \right.\right.\nonumber\\ &&\left.\left.- H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right) \right]\psi (x) \right\} \cdot \end{eqnarray} \tag{ 29 }$$

The average of the product $\psi (y)\bar \psi (z)$ over all field configurations is called two-point (one-particle) Green function G(y,z). We have

$$\begin{equation}G(y,z) = \left\langle {\psi (y)\bar \psi (z)} \right\rangle _{\psi,\bar \psi } = \psi (y)_c \bar \psi (z)_c = {\rm{i}}S(y,z).\end{equation} \tag{ 30 }$$

The average of the product $\psi (y_1)\psi (y_2)\bar\psi (z_1)\bar\psi(z_2)$ over all field configurations is called four-point (two-particle) Green function G(y₁,y₂ ;z₁,z₂ ). It has following explicit expression:

$$\begin{eqnarray} && G\left( {y_1,y_2 ;z_1,z_2 } \right) = \left\langle {\psi (y_1 )\psi (y_2 )\bar \psi (z_1 )\bar \psi (z_2 )} \right\rangle _{\psi,\bar \psi } \nonumber\\ & &= \frac{1}{{Z_0^\psi }}\int {[ {D\psi } ]\,[ {D\bar \psi } ]} \;\psi (y_1 )\psi (y_2 )\bar \psi (z_1 )\bar \psi (z_2 ) \nonumber\\ &&\times {\rm exp}\,\left\{ {{\rm{i}}\int {\rm{d}} x\bar \psi (x)\left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]\psi (x)} \right\},\quad\quad\quad \end{eqnarray} \tag{ 31 }$$

and can be expressed in terms of a functional derivative of the generating functional (13):

$$\begin{eqnarray} &&G(y_1,y_2 ;z_1,z_2 ) = \frac{1}{{Z_0^\psi }}\left. {\frac{{\delta ^4 Z_0^\psi \left[ {\eta,\bar \eta } \right]}}{{\delta \eta (z_2 )\delta \eta (z_1 )\delta \bar \eta (y_2 )\delta \bar \eta (y_1 )}}} \right|_{\eta = \bar \eta = 0}.\nonumber\end{eqnarray} \tag{ 32 }$$

Using explicit expression (26) of $Z_0 \left[ {\eta,\bar \eta } \right]$, after lengthy but standard calculations we obtain

$$\begin{eqnarray} G(y_1,y_2 ;z_1,z_2 ) &=& G(y_2,z_1 )G(y_1,z_2 )\nonumber\\ && - G(y_1,z_1 )G(y_2,z_2 ). \end{eqnarray} \tag{ 33 }$$

In terms of the contractions of two fermionic fields, formula (33) becomes

$$\begin{eqnarray} \left\langle\kern-1pt {\psi (y_1 )\psi (y_2 )\bar \psi (z_1 )\bar \psi (z_2 )} \right\rangle _{\psi,\bar \psi }\kern-1pt&=&\kern-1pt\psi (y_2 )_c \bar \psi (z_1 )_c\kern-0.5pt \psi (y_1 )_{c'}\kern-0.5pt \bar\psi (z_2 )_{c'} \nonumber\\ & -&\kern-0.5pt \psi (y_1 )_c\kern-0.5pt \bar \psi (z_1 )_c\kern-0.5pt \psi (y_2 )_{c'}\kern-0.5pt \bar\psi (z_2 )_{c'}.\nonumber\\ \end{eqnarray} \tag{ 34 }$$

Note that ψ (y₁) and $\bar \psi (z_1 )$ or ψ (y₂) and $\bar \psi (z_2 )$ are not two adjacent fields in the starting product of four fields in the lhs of formula (34). Let us rewrite formula (34) in another form

$$\begin{eqnarray} && \left\langle {\psi (y_1 )\psi (y_2 )\bar \psi (z_1 )\bar \psi (z_2 )} \right\rangle _{\psi,\bar \psi } \nonumber\\ &&= \left\langle {\psi (y_1 )_c \psi (y_2 )_{c'} \bar \psi (z_1 )_c \bar \psi (z_2 )_{c'} } \right\rangle _{\psi,\bar \psi } \nonumber\\ &&\quad+ \left\langle {\psi (y_1 )_c \psi (y_2 )_{c'} \bar \psi (z_1 )_{c'} \bar \psi (z_2 )_c } \right\rangle _{\psi,\bar \psi }, \end{eqnarray} \tag{ 35 }$$

where we have introduced the definition of the contractions between two non-adjacent fields: a contraction between two non-adjacent fields coincide with their contraction if the number of fermion fields between them is even, equal in magnitude but opposite in sign to their contraction if the number of fermionic fields between them is odd. Then formula (35) means that in order to find the average of the product in its lhs, one should form all possible contractions between the field ψ (y₁ ), ψ (y₂ ), and $\bar \psi (z_1 )$ and $\bar \psi (z_2 )$ and then sum up all obtained products of contractions. For any 2n-point Green function

$$\begin{eqnarray} &&G\left( {y_1,y_2, \ldots,y_n ;z_1,z_2, \ldots,z_n } \right) \nonumber\\ &&= \left\langle {\psi (y_1 )\psi (y_2 ) \cdots \psi (y_n )\bar \psi (z_1 )\bar \psi (z_2 ) \cdots \bar \psi (z_n )} \right\rangle _{\psi,\bar \psi },\nonumber\\ \end{eqnarray} \tag{ 36 }$$

we have the same rule: in order to derive the expression of the 2n-point Green function (36) one should form all possible contractions between each of the fields ψ (y₁),ψ (y₂), ...,ψ (y_n ) and each of the fields $\bar \psi (z_1 ),\bar \psi (z_2 ), \ldots,\bar \psi (z_n )$, and then sum up all obtained products of n! contractions. This is the well-known Wick theorem in quantum field theory.

To conclude this section we discuss a physical consequence of formula (27) and explain the physical meaning of the constants C_α in the rhs of formula (23). Consider the limit $\varepsilon \to + 0$ of the average of the product $\bar \psi ({\bf x},t + \varepsilon )\psi ({\bf x},t)$. It is the particle density n(x) :

$$\begin{eqnarray} n({\bf x}) &=& \mathop {\lim }\limits_{\varepsilon \to + 0} \left\langle {\bar \psi ({\bf x},t + \varepsilon )\psi ({\bf x},t)} \right\rangle _{\psi,\bar \psi } \nonumber\\ &=& - {\rm{i}}\mathop {\lim }\limits_{\varepsilon \to + 0} S({\bf x},{\bf x};\varepsilon ). \end{eqnarray} \tag{ 37 }$$

From expression (23) of S(x,y;ε ) it follows that

$$\begin{equation}n({\bf x}) = \sum\limits_\alpha {C_\alpha \left| {u_\alpha ({\bf x})} \right|^2 }\end{equation} \tag{ 38 }$$

Thus, C_α is the electron occupation number n_α at the quantum state with wave function u_α (x) : C_α =  n_α.

Consider now the system of interacting electrons in the confining potential U(x) of some nanostructure, the electron–electron interaction being the Coulomb repulsion. Denote u(x− x') the potential energy of the interaction between two electrons locating at two points x,x', and write

$$\begin{equation}u(x - x') = u({\bf x} - {\bf x'})\delta (x_0 - x'_0 ).\end{equation} \tag{ 39 }$$

The interacting electron system has following generating functional

$$\begin{eqnarray} Z &{=}&\!\! \int \!\!{[D\psi ]} [D\bar \psi ]\,{\rm exp}\!\left\{\! {{\rm{i}}\int \!\!{\rm{d}} x\bar \psi (x)\left[\! {{\rm{i}}\frac{\partial }{{\partial x_0 }} {-} H\!\left(\! { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \!\right)} \!\right]\!\psi (x)} \!\right\}\! \nonumber\\ &&\times {\rm exp}\!\left\{\! { {-} \frac{{\rm{i}}}{2}\int \!\!{{\rm{d}}x} \int \!\!{{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\bar \psi (x')\psi (x')} } \right\} \cdot \nonumber\\ \end{eqnarray} \tag{ 40 }$$

The Coulomb interaction term in the generating functional (40) is bilinear with respect to the electron density $\bar \psi (x)\psi (x)$. We linearize the interaction Hamiltonian with respect to the electron density $\bar \psi (x)\psi (x)$ by introducing a scalar field φ (x) playing the role of the order parameter of the collective oscillations in the electron gas and using following functional integral

$$\begin{equation}Z_0^\varphi = \int {[D\varphi ]} {\rm exp}\,\left\{ {\frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}x'\varphi (x)u(x - x')\varphi (x')} } \right\} \cdot\end{equation} \tag{ 41 }$$

Performing a shift of the functional integration variable

$$\begin{equation}\varphi (x) \to \varphi (x) - \bar \psi (x)\psi (x),\end{equation} \tag{ 42 }$$

we rewrite Z₀^φ in a new form

$$\begin{eqnarray} Z_0^\varphi & = & \int {[D\varphi ]} {\rm exp}\,\left\{ {\frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}x'\varphi (x)u(x - x')\varphi (x')} } \right\} \nonumber\\ & & \times {\rm exp}\left\{ { - {\rm{i}}\int {{\rm{d}}x} \int {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\varphi (x')} } \right\} \nonumber\\ && \times {\rm exp}\left\{\! {\frac{{\rm{i}}}{2}\kern-2pt\int {{\rm{d}}x}\kern-2pt\int\kern-2pt {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\bar \psi (x')\psi (x')} } \!\right\}\!.\nonumber\end{eqnarray} \tag{ 43 }$$

From this relation we derive the celebrated Hubbard–Stratonovich transformation

$$\begin{equation}\begin{array}{l} \displaystyle{\rm exp}\,\left\{ { - \frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\bar \psi (x')\psi (x')} } \right\} \\ \displaystyle= \frac{1}{{Z_0^\varphi }}\int {[D\varphi ]} {\rm exp}\,\left\{ {\frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}x'\varphi (x)u(x - x')\varphi (x')} } \right\} \\ \displaystyle \times {\rm exp}\,\left\{ { - {\rm{i}}\int {{\rm{d}}x} \int {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\varphi (x')} } \right\} \cdot \end{array} \end{equation} \tag{ 44 }$$

Substituting the expression in the rhs of relation (44) instead of its lhs into the rhs of formula (40), we transform the generating functional Z of the system into the form

$$\begin{eqnarray} && Z = \frac{1}{{Z_0^\varphi }}\int {[D\varphi ]} {\rm exp}\,\left\{ {\frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}x'\varphi (x)u(x - x')\varphi (x')} } \right\} \nonumber\\ &&\int\!\!\kern-1pt {[D\psi ][D\bar \psi ]} \,{\rm exp}\left\{\kern-2pt{ - {\rm{i}}\kern-2pt\int\kern-2pt {{\rm{d}}x} \kern-2pt\int\kern-2pt {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\varphi (x')}}\kern-1pt \right\} \nonumber\\ &&\times {\rm exp}\,\left\{ {{\rm{i}}\int {{\rm{d}}x\bar \psi (x)\left[ {{\rm{i}}\frac{\partial }{{\partial x_0 }} - H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)} \right]\psi (x)} } \right\} \cdot \end{eqnarray} \tag{ 45 }$$

According to the definition (29) of the average of an expression over all field configurations, formula (45) can be rewritten as follows

$$\begin{equation}\begin{array}{l} \displaystyle Z = \frac{{Z_0^\psi }}{{Z_0^\varphi }}\int {[D\varphi ]} {\rm exp}\,\left\{ {\frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}x'\varphi (x)u(x - x')\varphi (x')} } \right\} \\ \displaystyle \times \left\langle {{\rm exp}\,\left\{ { - {\rm{i}}\int {{\rm{d}}x} \int {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\varphi (x')} } \right\}} \right\rangle _{\psi,\bar \psi } \cdot \end{array} \end{equation} \tag{ 46 }$$

Expanding the exponential function

$$\begin{eqnarray*}&&{\rm exp}\,\left\{ { - {\rm{i}}\int {{\rm{d}}x} \int {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\varphi (x')} } \right\}\end{eqnarray*}$$

into a functional power series of the scalar field φ (x) and denoting F⁽ⁿ⁾ [φ ] the term of the nth order,

$$\begin{eqnarray} &&\left\langle {{\rm exp}\,\left\{ { - {\rm{i}}\int {{\rm{d}}x} \int {{\rm{d}}x'\bar \psi (x)\psi (x)u(x - x')\varphi (x')\;} } \right\}} \right\rangle _{\psi,\bar \psi } \nonumber\\ & =& F[\varphi ] = \sum\limits_{n = 0}^\infty {F^{(n)} [\varphi ]},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} &&F^{(n)} [\varphi ] {=} \frac{1}{{n!}}\left\langle {\left[\! { {-} {\rm{i}}\!\int \!{{\rm{d}}x} \!\int\! {{\rm{d}}x'\bar \psi (x)\psi (x)u(x {-} x')\varphi (x')} } \!\right]^n } \right\rangle _{\psi,\bar \psi }\!,\nonumber\\ \end{eqnarray} \tag{ 48 }$$

we obtain

$$\begin{equation}F^{(0)} [\varphi ] = 1,\end{equation} \tag{ 49 }$$

$$\begin{equation}F^{(1)} [\varphi ] = {\rm{i}}W^{(1)} [\varphi ],\end{equation} \tag{ 50 }$$

$$\begin{equation}W^{(1)} [\varphi ] = - \int {{\rm{d}}x} \int {{\rm{d}}x'n({\bf x})u(x - x')\varphi (x')},\end{equation} \tag{ 51 }$$

$$\begin{equation}F^{(2)} [\varphi ] = \frac{1}{{2!}}\left( {{\rm{i}}W^{(1)} [\varphi ]} \right)^2 + {\rm{i}}W^{(2)} [\varphi ],\end{equation} \tag{ 52 }$$

$$\begin{eqnarray} W^{(2)} [\varphi ] &=& \frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}x'\int {{\rm{d}}y} \int {{\rm{d}}y'S(x,y)S(y,x)} } \nonumber\\ &&\times u(x - x')u(y - y')\varphi (x')\varphi (y'),\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray} F^{(3)} [\varphi ] &=& \frac{1}{{3!}}\left( {{\rm{i}}W^{(1)} [\varphi ]} \right)^3 + \left( {{\rm{i}}W^{(1)} [\varphi ]} \right)\nonumber\\ &&\times\left( {{\rm{i}}W^{(2)} [\varphi ]} \right) + {\rm{i}}W^{(3)} [\varphi ], \end{eqnarray} \tag{ 54 }$$

$$\begin{equation}\begin{array}{l} \displaystyle W^{(3)} [\varphi ] = \frac{{\rm{i}}}{3}\!\!\int \!\!{{\rm{d}}x}\! \int \!\!{{\rm{d}}x'\int {{\rm{d}}y} \!\int \!\!{{\rm{d}}y'\!\int \!\!{{\rm{d}}z}\! \int \!\!{{\rm{d}}z'} S(x,y)S(y,z)}} \\ \displaystyle \times S(z,x) u(x - x')u(y - y')u(z - z')\varphi (x')\varphi (y')\varphi (z'),\end{array} \end{equation} \tag{ 55 }$$

and so on. Continuing similar calculations in higher orders, finally we have

$$\begin{equation}F[\varphi ] = {\rm exp}\,\left\{ {{\rm{i}}W[\varphi ]} \right\}\end{equation} \tag{ 56 }$$

with

$$\begin{equation}\;W[\varphi ] = \sum\limits_{n = 1}^\infty {W^{(n)} [\varphi ]}.\end{equation} \tag{ 57 }$$

Using expression (46) and relations (47)–(57), we obtain the generating functional Z in the form of a functional integral of the scalar field φ (x) only

$$\begin{equation}Z = \frac{{Z_0^\psi }}{{Z_0^\varphi }}\int {[D\varphi ]{\rm exp}\,\left\{ {{\rm{i}}I[\varphi ]} \right\}}\end{equation} \tag{ 58 }$$

with the effective action

$$\begin{eqnarray} I[\varphi ] &=& \frac{1}{2}\int {{\rm{d}}x} \int {{\rm{d}}y\varphi (x)u(x - y)\varphi (y)} \nonumber\\ &&- \int {{\rm{d}}x} \int {{\rm{d}}yn({\bf x})u(x - y)\varphi (y)} \nonumber\\ && + \frac{{\rm{i}}}{2}\int {{\rm{d}}x} \int {{\rm{d}}y\int {{\rm{d}}x'} \int {{\rm{dy'}}} S(x',y')S(y',x')} \nonumber\\ && \times u(x' - x)u(y' - y)\varphi (x)\varphi (y) \nonumber\\ && + \frac{{\rm{i}}}{3}\int {{\rm{d}}x} \int {{\rm{d}}y} \int {{\rm{d}}z} \int {{\rm{d}}x'} \int {{\rm{d}}y'} \int {{\rm{d}}z'} \nonumber\\ &&\times S(x',y')S(y',z')S(z',x') u(x' - x)\nonumber\\ &&\times u(y' - y)u(z' - z)\varphi (x)\varphi (y)\varphi (z) + \cdots. \end{eqnarray} \tag{ 59 }$$

In the approximation of the second order with respect to the scalar field φ (x), the effective action can be written as follows:

$$\begin{eqnarray} I_0 [\varphi ] &=& - \int {{\rm{d}}xv(x)\varphi (x)} \nonumber\\ &&+ \frac{1}{2}\int {{\rm{d}}x} \int {{\rm{d}}y} \varphi (x)A(x,y)\varphi (y),\end{eqnarray} \tag{ 60 }$$

where

$$\begin{equation}v(x) = \int {{\rm{d}}yu(x - y)n({\bf y})}\end{equation} \tag{ 61 }$$

and

$$\begin{eqnarray} A(x,y) &=& u(x - y) \nonumber\\ && + {\rm{i}}\int {{\rm{d}}x'} \int {\rm{d}}y'S(x',y')\nonumber\\ &&\times S(y',x')u(x' - x)u(y' - y).\end{eqnarray} \tag{ 62 }$$

From the principle of extreme action

$$\begin{equation}\frac{{\delta I_0 [\varphi ]}}{{\delta \varphi (x)}} = 0,\end{equation} \tag{ 63 }$$

we derive the dynamical equation for the scalar field φ₀ (x) corresponding to the extreme action:

$$\begin{equation}\int {{\rm{d}}yA(x,y)} \varphi _0 (y) = v(x).\end{equation} \tag{ 64 }$$

It has the solution

$$\begin{equation}\varphi _0 (y) = \int {{\rm{d}}yA^{ - 1} (y,x)} v(x),\end{equation} \tag{ 65 }$$

where A⁻¹ (y,x) is the kernel of the integral operator inverse to that with the kernel A(x,y):

$$\begin{equation} \int {{\rm{d}}zA(x,z)A^{ - 1} (z,y)} = \int {{\rm{d}}zA^{ - 1} (x,z)} A(z,y) = \delta (x - y).\end{equation} \tag{ 66 }$$

The extreme value of the action equals

$$\begin{equation}I_0 \left[ {\varphi _0 } \right] = - \frac{1}{2}\int {{\rm{d}}x} \int {{\rm{d}}y\varphi _0 (x)A(x,y)} \varphi _0 (y).\end{equation} \tag{ 67 }$$

The fluctuations of the scalar field φ (x) around the extreme field φ₀ (x) are described by the difference

$$\begin{equation}\zeta (x) = \varphi (x) - \varphi _0 (x).\end{equation} \tag{ 68 }$$

In terms of this new field ζ (x) the effective action $I_0 \left[ \varphi \right]$ has the quadratic form

$$\begin{equation}I_0 \left[ {\varphi _0 + \zeta } \right] = I_0 \left[ {\varphi _0 } \right] + \frac{1}{2}\int {{\rm{d}}x} \int {{\rm{d}}y\zeta (x)A(x,y)} \zeta (y).\end{equation} \tag{ 69 }$$

This new field ζ (x) is related directly to (but not identical to) quantum field whose quanta are the plasmons. Dynamical equation for the free plasmons is

$$\begin{equation}\int {{\rm{d}}yA(x,y)} \zeta (y) = 0.\end{equation} \tag{ 70 }$$

The terms of the orders higher than the second one give rise to the part of the effective action describing the plasmon–plasmon interaction. Note that due to the special shape of the nanostructure and/or the presence of the coordinate-dependent potential energy U(x), the plasma resonances are localized.

For testing the validity of the above-presented functional integral formalism, in this section we apply it to the simple case of a homogeneous electron gas in the three-dimensional space and compare its result with that of conventional theories. In this special case the energy potential U(x) in Hamiltonian (12) vanishes and function S(x, y) determined by inhomogeneous differential equation (22) depends only on the difference x–y of the variables x and y. This function S(x–y) has the Fourier transformation

$$\begin{eqnarray} S(x - y) &=& S({\bf x} - {\bf y};x_0 - y_0 ) \nonumber\\ &=& \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm{d}}\omega \;{\rm{e}}^{{\rm{i}}[{\bf k}({\bf x} - {\bf y}) - k_0 (x_0 - y_0 )]} {\tilde{S}}({\bf k},\omega )} \nonumber\\ \end{eqnarray} \tag{ 71 }$$

with the Fourier transform

$$\begin{equation}{\tilde{S}}({\bf k},\omega ) = \frac{{1 - n({\bf k})}}{{\omega - E({\bf k}) + {\rm{i}}o}} + \frac{{n({\bf k})}}{{\omega - E({\bf k}) - {\rm{i}}o}},\end{equation} \tag{ 72 }$$

where E(k) is the kinetic energy E(k) = k² /2m of free electron with momentum k and mass m, n(k) is the electron occupation number at the state with momentum k, 0 ⩽n(k) ⩽ 1. Function A(x, y) determined by formula (62) also depends only on the difference x–y and has the following Fourier transformation:

$$\begin{eqnarray} A(x - y) &=& A({\bf x} - {\bf y};x_0 - y_0 )\nonumber\\ & =& \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm{d}}\omega \;{\rm{e}}^{{\rm{i}}[{\bf k}({\bf x} - {\bf y}) - k_0 (x_0 - y_0 )]} \tilde A({\bf k},\omega )}\nonumber\\ \end{eqnarray} \tag{ 73 }$$

with the Fourier transform

$$\begin{equation}\tilde A({\bf k},\omega ) = \tilde u({\bf k}) + {\tilde{\Pi}} ({\bf k},\omega )\tilde u({\bf k})^2,\end{equation} \tag{ 74 }$$

where

$$\begin{eqnarray} {\tilde{\Pi}} ({\bf k},\omega ) &=& \frac{{\rm{i}}}{{(2\pi )^4 }}\int {{\rm{d}}{\bf p}} \int {\rm{d}} p_0 {\tilde{S}}\left( {{\bf p} + \frac{{\bf k}}{2},p_0 + \frac{\omega }{2}} \right)\nonumber\\ &&\times{\tilde{S}}\left( {{\bf p} - \frac{{\bf k}}{2},p_0 - \frac{\omega }{2}} \right)\end{eqnarray} \tag{ 75 }$$

and $\tilde u({\bf k})$ is the Fourier transform of $u({\bf x})$. From the dynamical equation (70) for the free plasmon there follows the algebraic equation determining the plasmon dispersion

$$\begin{equation}1 + u({\bf k}){\tilde{\Pi}} ({\bf k},\omega ) = 0.\end{equation} \tag{ 76 }$$

In the conventional quantum theory of many-body systems [35–37] the plasmon dispersion equation (76) can be derived by studying the wave vector and frequency dependence of the longitudinal dielectrical constant of a homogeneous system of itinerant electrons. In this medium the external charge with the density ρ^ext (x,t) creates electrical field with potential u^ext (x,t) determined by the Poisson equation

$$\begin{equation}\nabla ^2 u^{{\rm{ext}}} ({\bf x},t) = - 4\pi \rho ^{{\rm{ext}}} ({\bf x},t).\end{equation} \tag{ 77 }$$

This external field with potential u^ext (x,t) induces the polarization of the medium and creates a charge distribution with the density

$$\begin{equation}\rho ^{{\rm{in d}}} ({\bf x},t) = \int {{\rm d}{\bf x'}} \int {{\rm d}t'} \chi ({\bf x} - {\bf x'},t - t')u^{{\rm{ext}}} ({\bf x'},t'),\end{equation} \tag{ 78 }$$

and the total electrical field E(x,t) in the medium is determined by the equation

$$\begin{equation}\nabla {\bf E}({\bf x},t) = 4\pi \left[ {\rho ^{{\rm{ext}}} ({\bf x},t) + \rho ^{{\rm{ind}}} ({\bf x},t)} \right].\end{equation} \tag{ 79 }$$

By definition, the electrical induction vector D(x,t) is determined by the equation

$$\begin{equation}\nabla {\bf D}^{{\rm{ext}}} ({\bf x},t) = 4\pi \rho ^{{\rm{ext}}} ({\bf x},t).\end{equation} \tag{ 80 }$$

Its components are expressed in terms of those of the electrical field E(x,t) as follows:

$$\begin{equation}D_i ({\bf x},t) = \sum\limits_{j = 1}^3 {\int {{\rm{d}}{\bf x'}} } \int {{\rm{d}}t'} \varepsilon _{ij} ({\bf x} - {\bf x'},t - t')E_j ({\bf x'},t'),\end{equation} \tag{ 81 }$$

where ε_ij (x− x',t −t') are the components of the dielectrical constant tensor of the electron gas. Introducing the Fourier transformations

$$\begin{equation}\rho ^{{\rm{ext}},{\rm{ind}}} ({\bf x},t) = \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm{d}}\omega \;{\rm{e}}^{{\rm{i}}({\bf kx} - \omega t)} \tilde \rho ^{{\rm{ext}},{\rm{ind}}} ({\bf k},\omega )},\end{equation} \tag{ 82 }$$

$$\begin{equation}u^{{\rm{ext}}} ({\bf x},t) = \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm{d}}\omega \;{\rm{e}}^{{\rm{i}}({\bf kx} - \omega t)} \tilde u^{{\rm{ext}}} ({\bf k},\omega )},\end{equation} \tag{ 83 }$$

$$\begin{equation}E_i ({\bf x},t) = \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm{d}}\omega \;{\rm{e}}^{{\rm{i}}({\bf kx} - \omega t)} \tilde E_i ({\bf k},\omega )},\end{equation} \tag{ 84 }$$

$$\begin{equation}D_i ({\bf x},t) = \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm{d}}\omega \;{\rm{e}}^{{\rm{i}}({\bf kx} - \omega t)} \tilde D_i ({\bf k},\omega )},\end{equation} \tag{ 85 }$$

$$\begin{equation}\chi ({\bf x},t) = \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm d}\omega \;{\rm{e}}^{{\rm{i}}({\bf k}{\bf x} - \omega t)} \tilde \chi ({\bf k},\omega )},\end{equation} \tag{ 86 }$$

$$\begin{equation}\varepsilon _{i,j} ({\bf x},t) = \frac{1}{{(2\pi )^4 }}\int {{\rm{d}}{\bf k}} \int {{\rm{d}}\omega \;{\rm{e}}^{{\rm{i}}({\bf kx} - \omega t)} \tilde \varepsilon _{ij} ({\bf k},\omega )},\end{equation} \tag{ 87 }$$

we rewrite equations (77)–(81) in the form of algebraic equations for the Fourier transforms of physical quantities

$$\begin{equation}{\bf k}^2 \tilde v^{{\rm ext}} ({\bf k},\omega ) = 4\pi \tilde \rho ^{{\rm ext}} ({\bf k},\omega ),\end{equation} \tag{ 88 }$$

$$\begin{equation}\tilde \rho ^{{\rm ind}} ({\bf k},\omega ) = \tilde \chi ({\bf k},\omega )\tilde u^{{\rm ext}} ({\bf k},\omega ),\end{equation} \tag{ 89 }$$

$$\begin{equation} {\rm i}{\bf k}{\tilde{\bf E}}({\bf k},\omega ) = 4\pi \left[ {\tilde \rho ^{{\rm ext}} ({\bf k},\omega ) + \tilde \rho ^{{\rm ind}} ({\bf k},\omega )} \right]\;,\end{equation} \tag{ 90 }$$

$$\begin{equation}{\rm{i}}{\bf k}{\tilde{\bf D}}({\bf k},\omega ) = 4\pi \tilde \rho ^{{\rm{ext}}} ({\bf k},\omega ),\end{equation} \tag{ 91 }$$

$$\begin{equation}\tilde D_i ({\bf k},\omega ) = \sum\limits_{j = 1}^3 {\tilde \varepsilon _{ij} ({\bf k},\omega )\tilde E_j ({\bf k},\omega )} \;.\end{equation} \tag{ 92 }$$

In a homogeneous and isotropic medium tensor $\tilde \varepsilon _{ij} ({\bf k},\omega )$ has the general form

$$\begin{equation}\tilde \varepsilon _{ij} ({\bf k},\omega ) = \left( {\delta _{ij} - \frac{{k_i k_j }}{{{\bf k}^2 }}} \right)\tilde \varepsilon _{\rm{T}} ({\bf k},\omega ) + \frac{{k_i k_j }}{{{\bf k}^2 }}\tilde \varepsilon _{\rm{L}} ({\bf k},\omega ),\end{equation} \tag{ 93 }$$

where $\tilde \varepsilon _{\rm{T}} ({\bf k},\omega )$ and $\tilde \varepsilon _{\rm{L}} ({\bf k},\omega )$ are called (wave vector and frequency dependent) transverse and longitudinal dielectrical constants. From equations (88)–(93) it follows that

$$\begin{equation}\frac{1}{{\tilde \varepsilon _{\rm{L}} ({\bf k},\omega )}} = \frac{{{\rm{i}}{\bf k}{\tilde{\bf E}}({\bf k},\omega )}}{{{\rm{i}}{\bf k}{\tilde{\bf D}}({\bf k},\omega )}} = 1 + \frac{1}{{{\bf k}^2 }}\tilde \chi ({\bf k},\omega ).\end{equation} \tag{ 94 }$$

Using perturbation theory with respect to the Coulomb electron–electron interaction, in the random phase approximation [38–41] we obtained

$$\begin{equation}\tilde \chi ({\bf k},\omega ) = - \frac{{{\bf k}^2 4\pi e^2 \tilde \Pi ({\bf k},\omega )}}{{{\bf k}^2 + 4\pi e^2 \tilde \Pi ({\bf k},\omega )}}\;,\end{equation} \tag{ 95 }$$

where $\tilde \Pi ({\bf k},\omega )$ was determined by formula (75) and e is the electron charge. Substituting expression (95) into the rhs of equation (94), we have

$$\begin{equation}\tilde \varepsilon _{\rm{L}} ({\bf k},\omega ) = 1 + \frac{{4\pi e^2 }}{{{\bf k}^2 }}{\tilde{\Pi}} ({\bf k},\omega ).\end{equation} \tag{ 96 }$$

Because the plasmon resonance appears when the longitudinal dielectrical constant $\tilde \varepsilon _{\rm{L}} ({\bf k},\omega )$ vanishes, the wave vector dependence of the plasmon frequency ω is determined by the dispersion equation

$$\begin{equation}1 + \frac{{4\pi e^2 }}{{{\bf k}^2 }}{\tilde{\Pi}} ({\bf k},\omega ) = 0\;.\end{equation} \tag{ 97 }$$

Since

$$\begin{equation}\tilde u({\bf k}) = \frac{{4\pi e^2 }}{{{\bf k}^2 }}\;,\end{equation} \tag{ 98 }$$

dispersion equation (97) derived in the dielectrical formalism exactly coincides with formula (76) derived in functional integral formalism.

In earlier works on dielectrical formalism [38–41], it was shown that for a homogeneous and isotropic electron gas with average electron density n, at very small values of the wave vector, ${\bf k} \to 0$, and zero temperature (T = 0) we have

$$\begin{equation}{\rm Re}{\tilde{\Pi}} ({\bf k},\omega ) \approx - \frac{{k^2 }}{{\omega ^2 }}\frac{n}{m}\end{equation} \tag{ 99 }$$

and therefore

$$\begin{equation}\varepsilon _{\rm{L}} ({\bf k},\omega ) \approx 1 - \frac{{\omega _p^2 }}{{\omega ^2 }},\end{equation} \tag{ 100 }$$

where

$$\begin{equation}\omega _{\rm{p}} = \sqrt {\frac{{4\pi ne^2 }}{m}}\end{equation} \tag{ 101 }$$

is called the plasma frequency of the electron gas.

Formula (101) for the plasma frequency can be also easily derived from the classical equations of mechanics and electrodynamics by means of the canonical quantization rules in quantum mechanics. We start from the classical equations and denote R(x,t) the displacement of the electron located at the point x. Due to this displacement, the electron density changes from the mean value n to a varied value n(x,t). As a result, in the electron gas there arises some charge distribution with the density

$$\begin{equation}\rho ({\bf x},t) = - e\left[ {n({\bf x},t) - n} \right].\end{equation} \tag{ 102 }$$

This charge distribution induces electrical field with the potential

$$\begin{equation}u({\bf x},t) = \int {{\rm{d}}{\bf x'}\frac{{\rho ({\bf x'},t)}}{{\left| {{\bf x} - {\bf x'}} \right|}}}\end{equation} \tag{ 103 }$$

satisfying Poisson equation

$$\begin{equation}\nabla ^2 u({\bf x},t) = - 4\pi \rho ({\bf x},t).\end{equation} \tag{ 104 }$$

Due to the Coulomb repulsion between induced charges, the electron gas gains the potential energy

$$\begin{eqnarray} U & = & \frac{1}{2}\int {{\rm{d}}{\bf x}} \int {{\rm{d}}{\bf x'}\rho ({\bf x},t)\frac{1}{{\left| {{\bf x} - {\bf x'}} \right|}}} \rho ({\bf x'},t) \nonumber\\ & = & \frac{1}{2}\int {{\rm{d}}{\bf x}\rho ({\bf x},t)u({\bf x},t)}. \end{eqnarray} \tag{ 105 }$$

Due to the motion of each electron with the velocity ${\bf \dot R}({\bf x},t)$, the electron gas gains also the kinetic energy

$$\begin{equation}T = \frac{1}{2}\int {{\rm{d}}{\bf x}\,n({\bf x},t){\bf \dot R}({\bf x},t)^2 } \approx \frac{1}{2}n\int {{\rm{d}}{\bf x\,\dot R}({\bf x},t)^2 }.\end{equation} \tag{ 106 }$$

Denote N_S the number of the particles going through the boundary S of some volume V from the inside of V to its outside, and δN_V the decrease of the number of particles inside V. We have

$$\begin{equation} N_S {=} \oint\nolimits_S \!{n({\bf x},t){\bf R}({\bf x},t)\,{\rm{d}}{\bf S}} \approx n\oint\nolimits_S\! {{\bf R}({\bf x},t)\,{\rm{d}}{\bf S}} = n\int\nolimits_V \!{{\rm{d}}{\bf x}\nabla {\bf R}({\bf x},t)}\end{equation} \tag{ 107 }$$

and

$$\begin{equation}\delta N_V = \int {{\rm{d}}{\bf x}} [ {n - n({\bf x},t)} ].\end{equation} \tag{ 108 }$$

Because of the conservation of the total number of particles, there must be the equality

$$\begin{equation}N_S = \delta N_V.\end{equation} \tag{ 109 }$$

From equations (102) and (106)–(108) it follows that

$$\begin{equation}\rho ({\bf x},t) = en\nabla {\bf R}({\bf x},t).\end{equation} \tag{ 110 }$$

In order to perform the canonical quantization procedure it is convenient to decompose the functions ρ (x,t), u(x,t) and R(x,t) into the Fourier series with the use of a basis consisting of plane waves with discrete wave vectors. For this purpose as the basis we choose the complete set of plane waves satisfying the periodic boundary conditions at the opposite surfaces of a cube with volume V and normalize them in this cube. In this case we have following decompositions:

$$\begin{equation}\rho ({\bf x},t) = \frac{1}{{\sqrt V }}\sum\limits_{\bf k} {{\rm e}^{{\rm{i}}{\bf kx}} \rho _{\bf k} (t)},\end{equation} \tag{ 111 }$$

$$\begin{equation}u({\bf x},t) = \frac{1}{{\sqrt V }}\sum\limits_{\bf k} {{\rm e}^{{\rm{i}}{\bf kx}} u_{\bf k} (t)},\end{equation} \tag{ 112 }$$

$$\begin{equation}{\bf R}({\bf x},t) = \frac{1}{{\sqrt V }}\sum\limits_{\bf k} {{\rm{e}}^{{\rm{i}}{\bf kx}} } \sum\limits_{i = 1}^3 {{\bf \xi }_{\bf k}^{(i)} q_{\bf k}^{(i)} (t)},\end{equation} \tag{ 113 }$$

where ξ_k⁽ⁱ⁾ for each k are three unit vectors satisfying the conditions

$$\begin{equation}{\bf k\xi }_{\bf k}^{(i)} = \left\{ \begin{array}{@{}ll@{}} 0, & i = 1,2 \\ k, & i = 3, \\ \end{array} \right.,\end{equation} \tag{ 114 }$$

$k = | {\bf k}|$. Terms with i = 1, 2 are the transverse displacements, while that with i = 3 is the longitudinal displacement along the direction of the wave vector k. In terms of the Fourier transforms ρ_k (t), u_k (t) and q_k⁽ⁱ⁾ (t), equations (104) and (110) become

$$\begin{equation}u_{\bf k} (t) = \frac{{4\pi }}{{{\bf k}^2 }}\rho _{\bf k} (t)\end{equation} \tag{ 115 }$$

and

$$\begin{equation}\rho _{\bf k} (t) = {\rm{i}}enkq_{\bf k}^{(3)} (t).\end{equation} \tag{ 116 }$$

Therefore the potential energy (105) and the kinetic energy (106) are expressed in terms of the generalized coordinates q_k⁽ⁱ⁾ (t) of the system as follows:

$$\begin{equation}U = 2\pi (en)^2 \sum\limits_{\bf k} {q_{\bf k}^{(3)} (t)^ * q_{\bf k}^{(3)} (t)},\end{equation} \tag{ 117 }$$

$$\begin{equation}T = \frac{1}{2}nm\sum\limits_{\bf k} {\sum\limits_{i = 1}^3 {\dot q_{\bf k}^{(i)} (t)^ * \dot q_{\bf k}^{(i)} (t)} }.\end{equation} \tag{ 118 }$$

The Lagrangian of electron gas as a classical mechanical system has the form

$$\begin{eqnarray} L &=& T - U \nonumber\\ &=& \frac{1}{2}\sum\limits_{\bf k} {\left\{ {nm\sum\limits_{i = 1}^3 {\dot q_{\bf k}^{(i)} (t)^ * \dot q_{\bf k}^{(i)} (t) - 4\pi e^2 n^2 q_{\bf k}^{(3)} (t)^ * q_{\bf k}^{(3)} (t)} } \right\}}.\nonumber\\ \end{eqnarray} \tag{ 119 }$$

The equations of motion are

$$\begin{equation}\frac{{{\rm{d}}^2 q_{\bf k}^{(i)} (t)}}{{{\rm{d}}t^2 }} + 4\pi e^2 \frac{n}{m}\delta _{i3} q_{\bf k}^{(i)} (t) = 0.\end{equation} \tag{ 120 }$$

It follows that the transverse coordinates q_k⁽ⁱ⁾ (t), i = 1, 2 do not oscillate, and longitudinal one q_k⁽³⁾ (t) periodically oscillates with the angular frequency $\omega _{\rm p} = \sqrt {{{4\pi e^2 n}/m}}$ which is the plasma frequency of the electron gas as a classical mechanical system. Because the transverse coordinates q_k⁽¹⁾ (t) and q_k⁽²⁾ (t) do not oscillate, we discard them and replace the notation q_k⁽³⁾ (t) by a simpler one q_k (t). The Lagrangian of the electron gas becomes

$$\begin{equation}L = \frac{1}{2}\sum\limits_{\bf k} {\left\{\!{\;nm\dot q_{\bf k} (t)^ * \dot q_{\bf k} (t) - 4\pi e^2 n^2 q_{\bf k} (t)^ * q_{\bf k} (t)} \right\}}.\end{equation} \tag{ 121 }$$

The canonical momentum p_k (t) corresponding to the generalized coordinate q_k (t) is

$$\begin{equation}p_{\bf k} (t) = \frac{{\partial L}}{{\partial \dot q_{\bf k} (t)}} = \frac{1}{2}nm\dot q_{\bf k} (t)^ *.\end{equation} \tag{ 122 }$$

The Hamiltonian of the system equals

$$\begin{eqnarray} H &=& \sum\limits_{\bf k} p_{\bf k} (t)\dot q_{\bf k} (t) - L \nonumber\\ &=& \frac{1}{2}\sum\limits_{\bf k} {\left\{ {\frac{1}{{mn}}p_{\bf k} (t)^ * p_{\bf k} (t) + 4\pi e^2 n^2 q_{\bf k} (t)^ * q_{\bf k} (t)} \right\}}.\nonumber\\ \end{eqnarray} \tag{ 123 }$$

Because ρ (x,t), u(x,t) and R(x,t) are real functions, the canonical coordinates and momenta q_k (t) and p_k (t) have the following property:

$$\begin{equation}q_{\bf k} (t)^ * = q_{ - {\bf k}} (t), \quad p_{\bf k} (t)^ * = p_{ - {\bf k}} (t).\end{equation} \tag{ 124 }$$

The canonical quantization procedure consists of replacing canonical coordinates and momenta q_k (t) and p_k (t) by corresponding operators $\hat q_{\bf k}$ and $\hat p_{\bf k}$, satisfying canonical commutation relations

$$\begin{equation} \left[ {\hat q_{\bf k},\hat q_{\bf l} } \right] = \left[ {\hat p_{\bf k},\hat p_{\bf l} } \right] = 0,\quad \left[ {\hat p_{\bf k},\hat q_{\bf l} } \right] = - {\rm{i}}\delta _{{\bf kl}}. \end{equation} \tag{ 125 }$$

The complex conjugate canonical variables q_k (t) * and p_k (t) * are then replaced by corresponding hermitian conjugate operators

$$\begin{equation}\hat q_{\bf k}^ + = \hat q_{ - {\bf k}}, \quad \hat p_{\bf k}^ + = \hat p_{ - {\bf k}}.\end{equation} \tag{ 126 }$$

Setting

$$\begin{equation}\begin{array}{l} \displaystyle\hat a_{\bf k} = \frac{1}{{\sqrt {2mn\omega _{\rm p} } }}\hat p_{\bf k} - {\rm{i}}\sqrt {\frac{{mn\omega _{\rm{p}} }}{2}} \hat q_{\bf k}^ +, \\ [12pt] \displaystyle\hat a_{\bf k}^ + = \frac{1}{{\sqrt {2mn\omega _{\rm p} } }}\hat p_{\bf k}^ + + {\rm{i}}\sqrt {\frac{{mn\omega _{\rm{p}} }}{2}} \hat q_{\bf k}, \end{array} \end{equation} \tag{ 127 }$$

we obtain plasmon destruction and creation operators $\hat a_{\bf k}$ and $\hat a_{\bf k}^ +$ satisfying commutation relations

$$\begin{equation} \left[ {\hat a_{\bf k},\hat a_{\bf l} } \right] = \left[ {\hat a_{\bf k}^ +,\hat a_{\bf l}^ + } \right] = 0, \quad \left[ {\hat a_{\bf k},\hat a_{\bf l}^ + } \right] = \delta _{{\bf kl}}. \end{equation} \tag{ 128 }$$

Operators $\hat q_{\bf k}$ and $\hat p_{\bf k}$ are expressed in terms of destruction and creation operators as follows:

$$\begin{equation} \begin{array}{l} \displaystyle\hat q_{\bf k} = \frac{{\rm{i}}}{{\sqrt {2mn\omega _{\rm{p}} } }}\left( {\hat a_{ - {\bf k}} - \hat a_{\bf k}^ + } \right), \\ \displaystyle\hat p_{\bf k} = \sqrt {\frac{{mn\omega _{\rm{p}} }}{2}} \left( {\hat a_{\bf k} - \hat a_{ - {\bf k}}^ + } \right). \end{array} \end{equation} \tag{ 129 }$$

Substituting these expressions into the rhs of formula (123), we obtain the Hamiltonian

$$\begin{equation}H = \sum\limits_{\bf k} {\omega _{\rm{p}} \left( {\hat a_{\bf k}^ + \hat a_{\bf k} + \frac{1}{2}} \right)} \cdot\end{equation} \tag{ 130 }$$

This means that all plasmons have the same energy ω_p.

After checking the validity of the dynamical equation (70) for the free plasmons by comparing the plasmon dispersion equation (76) derived from this general dynamical equation for the special case of the homogeneous and isotropic three-dimensional electron gas with that previously obtained from conventional theories, equation (97), in this section, we demonstrate how to establish the explicit form of the quantum field, whose quanta are the plasmons, from the bosonic field ζ (x) in the expression (69) of the effective action. ζ (x) is related but does not coincide with the canonical quantum field of plasmons which we call the canonical quantum plasmonic field. This canonical quantum field must obey the canonical quantization rules, while ζ (x) does not.

For sounding out the way to reach our aim, we start from a simple case: the homogeneous isotropic three-dimensional electron gas. In this case the function A(x,y) in expression (69) of the effective action depends only on the coordinate difference x–y. Its Fourier component $\tilde A({\bf k},\omega )$ was given in equations (74) and (75). For small values of k and k_F such that k² ≪k_F² ⩽m²,k² ≪ ω_p² we have the approximate expression

$$\begin{equation}\tilde A({\bf k},\omega) = \frac{{4\pi e^2 }}{k^2\omega ^2}\{ {\omega ^2 - \omega _{\rm p}^2 - \gamma ^2 k^2 }\},\end{equation} \tag{ 131 }$$

where

$$\begin{equation}\gamma ^2 = \frac{3}{5}\frac{{k_{\rm{F}}^2 }}{{m^2 }} \cdot\end{equation} \tag{ 132 }$$

Introducing the Fourier transform $\tilde \zeta ({\bf k},\omega )$ of ζ (x),

$$\begin{equation}\zeta (x) = \frac{1}{{2\pi }}\int {{\rm{d}}\omega } \frac{1}{{\left( {2\pi } \right)^3 }}\int {{\rm{d}}{\bf k}\,{\rm{e}}^{ - {\rm{i}}(\omega t - {\bf kx})} \tilde \zeta ({\bf k},\omega )},\end{equation} \tag{ 133 }$$

and setting

$$\begin{equation}\tilde \sigma ({\bf k},\omega ) = \frac{{\sqrt {4\pi e^2 } }}{{\omega k}}\tilde \zeta ({\bf k},\omega ),\end{equation} \tag{ 134 }$$

we rewrite the effective action in the form

$$\begin{eqnarray} I_0^{{\rm{eff}}} [\sigma ] &=& \frac{1}{{2\pi }}\int {{\rm{d}}\omega } \frac{1}{{\left( {2\pi } \right)^3 }}\int {\rm{d}}{\bf k}\frac{1}{2}\tilde \sigma ({\bf k},\omega )[ {\omega ^2 - \omega _{\rm p}^2 - \gamma ^2 k^2 } ]\nonumber\\ &&\times\tilde \sigma ({\bf k},\omega ).\end{eqnarray} \tag{ 135 }$$

Considering $\tilde \sigma ({\bf k},\omega )$ as the Fourier transform of a scalar quantum field σ (x,t), we represent the effective action (135) in terms of this scalar field

$$\begin{eqnarray} I_0^{{\rm{eff}}} \left[ \sigma \right] &=& \frac{1}{2}\int {{\rm{d}}t} \int {{\rm{d}}{\bf x}} \left\{ \left[ {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right]^2 \right.\nonumber\\ &&\left.- \gamma ^2 \left[ {\nabla \sigma ({\bf x},t)} \right]^2 - \omega _{\rm p}^2 \sigma ({\bf x},t)^2 \vphantom{{\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}}}\vphantom{\left[ {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right]^2}\right\}.\end{eqnarray} \tag{ 136 }$$

It is similar to the expression of the action for a free relativistic scalar Klein–Gordon quantum field [42–44], except for a scaling factor γ of the spatial coordinates. The Lagrangian of the field σ (x,t) can be written in the standard form

$$\begin{eqnarray} &&L = \frac{1}{2}\int {{\rm{d}}{\bf x}} \left\{ {\left[ {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right]^2 - \gamma ^2 \left[ {\nabla \sigma ({\bf x},t)} \right]^2 - \omega _{\rm{p}}^2 \sigma ({\bf x},t)^2 } \right\}.\nonumber\\ \end{eqnarray} \tag{ 137 }$$

The canonical momentum π (x,t) corresponding to the canonical coordinate σ (x,t), by definition, is the following operator:

$$\begin{equation}\pi ({\bf x},t) = \frac{{\delta L}}{{\delta \left[ {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right]}} = \frac{{\partial \sigma ({\bf x},t)}}{{\partial t}} \cdot\end{equation} \tag{ 138 }$$

Canonical quantum operators σ (x,t) and π (x,t) must satisfy following equal-time canonical quantization rules:

$$\begin{equation}\begin{array}{l} \left[ {\sigma ({\bf x},t),\sigma ({\bf y},t)} \right] = \left[ {\pi ({\bf x},t),\pi ({\bf y},t)} \right] = 0, \\ \left[ {\pi ({\bf x},t),\sigma ({\bf y},t)} \right] = - {\rm{i}}\delta ({\bf x} - {\bf y}). \end{array} \end{equation} \tag{ 139 }$$

It is straightforward to extend the above-presented quantization procedure to the general case of the electron gas in a nanostructure with arbitrary shape and nanosize. In this case, function A(x,y) separately depends on both spatial coordinates x, y and has the form

$$\begin{eqnarray} A(x,y) &=& A({\bf x},{\bf y};x_0 - y_0 ) = \frac{1}{{2\pi }}\int {{\rm{d}}\omega \;{\rm{e}}^{ - {\rm{i}}\omega (x_0 - y_0 )} } \nonumber\\ &&\times\left[ u({\bf x} - {\bf y}) + \int {{\rm{d}}{\bf x'}} \int {\rm{d}}{\bf y'}\right.\nonumber\\ &&\left.\times u({\bf x} - {\bf x'})u({\bf y} - {\bf y'})\tilde \Pi ({\bf x'},{\bf y'};\omega ) \vphantom{\int }\right]. \end{eqnarray} \tag{ 140 }$$

Denote by u_α(x) the eigenfunctions of Hamiltonian $H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)$ with corresponding eigenvalues E_α:

$$\begin{equation}H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)u_\alpha ({\bf x}) = E_\alpha u_\alpha ({\bf x}).\end{equation} \tag{ 141 }$$

It is straightforward to show that

$$\begin{eqnarray} \tilde \Pi ({\bf x},{\bf y};\omega ) &=& \sum\limits_\alpha {u_\alpha ({\bf x})u_\alpha ({\bf y})^ * u_\beta ({\bf y})u_\beta ({\bf x})^ * } \nonumber\\ &&\times \left\{ {\frac{{(1 - n_\alpha )n_\beta }}{{E_\alpha - E_\beta - \omega - {\rm i}o}} + \frac{{(1 - n_\beta )n_\alpha }}{{E_\beta - E_\alpha + \omega - {\rm{i}}o}}} \right\}, \nonumber\\ \end{eqnarray} \tag{ 142 }$$

where n_α and n_β are the occupation numbers at the states with energies E_α and E_β, respectively. Consider the case when Hamiltonian $H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)$ is real in the sense that

$$\begin{equation}H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right)^ * = H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right).\end{equation} \tag{ 143 }$$

For example

$$\begin{eqnarray*}&&H\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}},{\bf x}} \right) = \frac{1}{{2m}}\left( { - {\rm{i}}\frac{\partial }{{\partial {\bf x}}}} \right)^2 + U({\bf x}).\end{eqnarray*}$$

Then both u_α (x) and u_α (x) *are two eigenfunctions of $H( - {\rm{i}}{\textstyle{\partial \over {\partial {\bf x}}}},{\bf x})$ with one and the same eigenvalue E_α. In this case the expression in the rhs of equation (142) can be transformed into a new form

$$\begin{eqnarray} \tilde \Pi ({\bf x},{\bf y};\omega ) &=& - 2\sum\limits_{\alpha,\beta } {u_\alpha ({\bf x})u_\alpha ({\bf y})^ * u_\beta ({\bf y})u_\beta ({\bf x})^ * (1 - n_\alpha )n_\beta } \nonumber\\ &&\times\frac{{E_\alpha - E_\beta }}{{\omega ^2 - (E_\alpha - E_\beta - io)^2 }} \cdot \end{eqnarray} \tag{ 144 }$$

Suppose that for all pairs of states with energies E_α and E_β we always have

$$\begin{eqnarray*} &&(E_\alpha - E_\beta )^2 \ll \omega ^2. \end{eqnarray*}$$

Then $\tilde \Pi ({\bf x},{\bf y};\omega )$ has approximate expression

$$\begin{eqnarray}\tilde \Pi ({\bf x},{\bf y};\omega ) &=& - \frac{2}{{\omega ^2 }}\sum\limits_{\alpha,\beta } u_\alpha ({\bf x})u_\alpha ({\bf y})^ * u_\beta ({\bf y})u_\beta ({\bf x})^ * (1 - n_\alpha )n_\beta \nonumber\\ &&\times (E_\alpha - E_\beta ) \cdot\end{eqnarray} \tag{ 145 }$$

Denote $\tilde \Pi ({\bf k},{\bf l};\omega )$ the Fourier transform of $\tilde \Pi ({\bf x},{\bf y};\omega )$ with respect to the spatial coordinates

$$\begin{equation}\tilde \Pi ({\bf x},{\bf y};\omega ) = \frac{1}{{(2\pi )^3 }}\int {{\rm{d}}{\bf k}} \frac{1}{{(2\pi )^3 }}\int {{\rm{d}}{\bf l}} \,{\rm{e}}^{{\rm{i}}({\bf kx} - {\bf ly})} \tilde \Pi ({\bf k},{\bf l};\omega ).\end{equation} \tag{ 146 }$$

In terms of $\tilde \Pi ({\bf k},{\bf l};\omega )$ the effective action has the form

$$\begin{eqnarray} I_0 [\varphi ] &=& I_0 [\varphi _0 ] {+} \frac{1}{{2\pi }}\int \!\!{{\rm{d}}\omega } \!\left\{\! {\frac{1}{{(2\pi )^3 }}\!\int\!\! {{\rm{d}}{\bf k}} \frac{1}{2}\tilde \zeta ({\bf k},\omega )^+ \tilde u({\bf k})\tilde \zeta ({\bf k},\omega )} \right. \nonumber\\ &&+ \frac{1}{{(2\pi )^3 }}\!\int\!\! {{\rm{d}}{\bf k}} \frac{1}{{(2\pi )^3 }}\!\int \!\! {{\rm{d}}{\bf l}} \frac{1}{2}\tilde \zeta ({\bf k},\omega )^ + \tilde u({\bf k})\nonumber\\ &&\left.\times\tilde \Pi ({\bf k},{\bf l};\omega )\tilde u({\bf l})\tilde \zeta ({\bf l},\omega ) \vphantom{\frac{1}{2}}\right\} \cdot \end{eqnarray} \tag{ 147 }$$

Setting

$$\begin{equation}\sqrt {\tilde u({\bf k})} \;\tilde \zeta ({\bf k},\omega ) = \omega \tilde \sigma ({\bf k},\omega ),\end{equation} \tag{ 148 }$$

$$\begin{equation}\omega ^2 \sqrt {\tilde u({\bf k})} \;\tilde \Pi ({\bf k},{\bf l};\omega )\sqrt {\tilde u({\bf l})} = - \tilde K({\bf k},{\bf l}),\end{equation} \tag{ 149 }$$

we obtain the effective action in a new form

$$\begin{eqnarray} I_0 [\varphi ] &= & I_0 [\varphi _0 ] + I_{\rm{0}}^{{\rm{eff}}} [\sigma ], \nonumber\\ I_0^{{\rm{eff}}} [\sigma ] &= & \frac{1}{{2\pi }}\int {{\rm{d}}\omega } \left\{ {\frac{1}{{(2\pi )^3 }}\int {{\rm{d}}{\bf k}} \frac{1}{2}\omega ^2 \tilde \sigma ({\bf k},\omega )^ + \tilde \sigma ({\bf k},\omega )} \right. \nonumber\\ && \left. { {-} \frac{1}{{(2\pi )^3 }}\int \!{{\rm{d}}{\bf k}} \frac{1}{{(2\pi )^3 }}\int \!{{\rm{d}}{\bf l}} \tilde \sigma ({\bf k},\omega )^ + \tilde K({\bf k},{\bf l})\tilde \sigma ({\bf l},\omega )} \right\}\!. \nonumber\\ \end{eqnarray} \tag{ 150 }$$

Considering $\tilde K({\bf k},{\bf l})$as the Fourier transform of a function K(x,y) and $\tilde \sigma ({\bf k},\omega )$ as that of a scalar field σ (x,t),

$$\begin{equation}K({\bf x},{\bf y}) = \frac{1}{{(2\pi )^3 }}\int {{\rm{d}}{\bf k}} \frac{1}{{(2\pi )^3 }}\int {{\rm{d}}{\bf l}\,} {\rm{e}}^{{\rm{i}}({\bf kx} - {\bf ly})} \tilde K({\bf k},{\bf l}),\end{equation} \tag{ 151 }$$

$$\begin{equation}\sigma ({\bf x},t) = \frac{1}{{2\pi }}\int {{\rm{d}}\omega } \frac{1}{{(2\pi )^3 }}\int {{\rm{d}}{\bf k}} \,{\rm{e}}^{ - {\rm{i}}(\omega t - {\bf kr})} \tilde \sigma ({\bf k},\omega ),\end{equation} \tag{ 152 }$$

we can rewrite the effective action (150) in the form

$$\begin{eqnarray} I_0^{{\rm{eff}}} [\sigma ] &=& \frac{1}{2}\int {{\rm{d}}t} \left\{\int {\mathop{\rm d}\nolimits} {\bf x}\left( {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right)^2 \right.\nonumber\\ &&\left.- \int {{\rm{d}}{\bf x}} \int {{\rm{d}}{\bf y}\,\sigma ({\bf x},t)K({\bf x},{\bf y})\sigma ({\bf y},t)} \vphantom{\left( {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right)} \right\} \cdot\end{eqnarray} \tag{ 153 }$$

In the rhs of equation (153), the last term contains the linear integral transformation of the field σ (x,t):

$$\begin{eqnarray*}&&\sigma ({\bf x},t) \to \int {{\rm{d}}{\bf y}} K({\bf x},{\bf y})\sigma ({\bf y},t). \end{eqnarray*}$$

It is a non-local linear transformation. This non-locality property is a consequence as well as an exhibition of the essential structure of localized plasmon as the extended quasiparticles resulted in the collective oscillations of the electron gas in some potential field U(x).

Lagrangian of the scalar field σ (x,t) has the form

$$\begin{eqnarray} &&L {=} \frac{1}{2}\left\{\! {\int \!\!{{\rm{d}}{\bf x}\left(\! {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \!\right)^2 {-} \!\int\!\! {{\rm{d}}{\bf x}} \!\int\!\! {{\rm{d}}{\bf y}\,\sigma ({\bf x},t)K({\bf x},{\bf y})\sigma ({\bf y},t)} } } \!\right\} \cdot\nonumber\\ \end{eqnarray}\noindent \tag{ 154 }$$

By definition, the canonical momentum corresponding to the canonical coordinate σ (x,t) is the operator

$$\begin{equation}\pi ({\bf x},t) = \frac{{\delta L}}{{\delta \left( {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right)}} = \frac{{\partial \sigma ({\bf x},t)}}{{\partial t}} \cdot\end{equation} \tag{ 155 }$$

Hamiltonian of the scalar field σ (x,t) equals

$$\begin{eqnarray} H &=& \int {{\rm{d}}{\bf x}\,\pi ({\bf x},t)} \frac{{\partial \sigma ({\bf x},t)}}{{\partial t}} - L = \frac{1}{2}\left\{\int {\rm{d}}{\bf x}\left( {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right)^2 \right.\nonumber\\ &&\left.+ \int {{\rm{d}}{\bf x}} \int {{\rm{d}}{\bf y}\,\sigma ({\bf x},t)K({\bf x},{\bf y})\sigma ({\bf y},t)} \vphantom{\left( {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right)}\vphantom{\left( {\frac{{\partial \sigma ({\bf x},t)}}{{\partial t}}} \right)^2}\right\} \cdot \end{eqnarray}\noindent \tag{ 156 }$$

Canonical coordinate and momentum σ (x,t) and π (x,t) must satisfy the canonical quantization rules (139). Linear integral transformation with the kernel K(x,y) is a definitely positive operator with eigenvectors f_i (x) and eigenvalues ε_i² :

$$\begin{eqnarray} \int {{\rm d}{\bf y}} \;K({\bf x},{\bf y})f_i ({\bf y}) &=& \varepsilon _i^2 f_i ({\bf x}), \nonumber\\ \int {{\rm d}{\bf y}} \;K({\bf x},{\bf y})f_i ({\bf y})^ * &=& \varepsilon _i^2 f_i ({\bf x})^ *, \\ \int {{\rm d}{\bf x}\;} f_i ({\bf x})^ * f_j ({\bf x}) &=& \delta _{ij}. \nonumber\end{eqnarray}\noindent \tag{ 157 }$$

Let us expand the scalar field σ (x,t) over the complete set of eigenvectors f_i (x) and f_i (x) *,

$$\begin{equation}\sigma ({\bf x},t) = \sum {\frac{1}{{\sqrt {2\varepsilon _i } }}} \left\{ {a_i \,{\rm{e}}^{ - {\rm{i}}\varepsilon _i t} f_i ({\bf x}) + a_i^ + \,{\rm{e}}^{{\rm{i}}\varepsilon _i t} f_i ({\bf x})^ * } \right\},\end{equation} \tag{ 158 }$$

where a_i and a_i ⁺ are the destruction and creation operators of the plasmon with energy ε_i,

$$\begin{equation}\begin{array}{l} {[}a_i, a_j] = [a_i^+, a_j^+] = 0, \\ [6pt] { }[a_i, a_j^+] = \delta _{ij}. \end{array} \end{equation} \tag{ 159 }$$

Then the Hamiltonian (156) of the scalar field σ (x,t) becomes

$$\begin{equation}H = \sum\limits_i {\varepsilon _i \left( {a_i^ + a_i + \frac{1}{2}} \right)}.\end{equation} \tag{ 160 }$$

A new approach to the quantum dynamics of plasmons in nanomaterials based on the functional integral technique was proposed. Starting from the generating functional of the fermionic field ψ (x) and its hermitian conjugate $\bar \psi (x)$ as the Grassman variables and introducing a new scalar field φ (x) playing the role of the order parameter characterizing the collective oscillations of the electron gas in the nanomaterials, we transform the generating functional of the electron gas into that of scalar field φ (x). The latter is then rewritten in another form convenient for the derivation of the dynamical equation of the plasmonic field. In the special case of a homogeneous and isotropic electron gas in the three-dimensional space the plasmon energy calculated from the above-mentioned dynamical equation exactly coincides with that previously obtained from conventional theories. However, the conventional theories are not able to be applied to the study of plasmons in nanostructures with different sizes and shapes, while, in principle, the functional integral approach can be effectively applied to the study of plasmons in nanostructures with complicated shapes.

However, in this beginning article of a series of theoretical works on the functional integral formalism of the quantum theory of plasmons, which can also be called briefly quantum plasmodynamics (QPD), we just presented the basics of the theory and showed how to derive the dynamical equation determining the plasmon dispersion law as well as elaborated the mathematical method for establishing the canonical quantum plasmonic field. It is then straightforward to find the plasmon destruction and creation operators and effective Hamiltonians of the interactions of plasmons with other quasiparticles in nanostructures from first principles as well as to investigate the physical processes with the participation of plasmons.

The authors express their gratitude to Vietnam Academy of Science and Technology for the support.