Comprehensive Microscopic Theory for Coupling of Longitudinal–Transverse Fields and Individual–Collective Excitations

Elucidation of the light–matter interaction is one of the core research topics of contemporary physics. Recently, the behaviors of light-induced plasmons in metals have attracted considerable interest Tame13 . A plasmon is a quantum of collective electron motion due to the electron–electron interaction, and has been clearly observed in experiment by Powell and Swan Powell59 . In bulk metals, a plasmon has a longitudinal (L) component only and, hence, cannot be excited by transverse (T) electromagnetic (EM) fields. To induce a plasmon using light, the L component of the EM field, e.g., the evanescent mode on a surface, is required book:NovotnyHecht . Cardinal methods to excite surface plasmon-polaritons (SPPs) on a 2D metal surface have been developed by Otto Otto68 and Kretschmann Kretschmann68 . The induced SPPs cause a charge–density spatial distribution on the surface, which in turn generates L fields. Therefore, the EM fields and induced plasmons must be determined self-consistently Feibelman82 ; Gerhardts84 . In metal nanostructures, the SPPs are localized and enhance the L fields significantly. Thus, SPPs are sensitive to surface states and are, therefore, applied in sensors for gases, molecules, and biological matter Garcia11 ; Zeng11 ; Gandhi19 .

One of the fascinating applications of light-induced plasmons is hot-carrier generation  Govorov14 ; Brongersma15 and injection to combined materials White12 ; Govorov13 ; Kumar19 , which can be utilized for photocatalysis Ueno16 , photodetection Li17 , photocarrier injection Tatsuma17 , and photovoltaics Clavero14 . The hot-carrier generation processes involving localized SPPs are categorized into two types: those achieved through plasmon relaxation or through coherent coupling between plasmons and individual single-particle excitations (SPEs). Some theoretical studies have elucidated the relevant mechanisms of the plasmon relaxation approach. The first-principles calculation has been examined for hot-carrier generation on 2D metal surfaces Sundararaman14 ; Bernardi15 , where the interband and intraband excitations make dominant and tunable contributions to the generation, respectively. Hot-carrier generation for metallic nanostructures has also been investigated based on a phenomenological model with several relaxation times Zhang14 ; Ross15 ; Besteiro17 and using density functional theory Manjavacas14 . It was found that the nanostructure sizes and shapes influence the plasmon via the electronic wavefunction and the enhanced electric field due to the hot spots. In those studies, a unidirectional energy transfer from the plasmon(-like) state to the hot carriers by the electron SPEs was discussed Zhang14 ; Ross15 ; Besteiro17 ; Manjavacas14 ; Chapkin18 . As noted above, the second type of hot-carrier generation is due to coherent coupling between the plasmons and individual SPEs. Recently, Ma et al. explicitly discussed the interplay between the plasmons and SPEs in a nanocluster of Ag atoms using a real-time simulation based on time-dependent density functional theory Ma15 . The on- and off-resonant conditions between the plasmon and SPE were found to change the energy transfer processes. For the off-resonant condition, a coherent Rabi oscillation between the plasmon and SPE was found. You et al. also elucidated the bidirectional energy transfer between the plasmon and SPE based on a model Hamiltonian with the Coulomb interaction You18 , considering excited, injected, and extracted electrons.

In the present work, we focus on the fact that the coupling between the plasmons and SPEs via the T field, as well as the treatment of the microscopic nonlocal response involving the T field, has been missed in the previous frameworks, which is particularly related with the latter case in the above discussion. Recent ab initio calculations successfully implemented coherent coupling or hybridization between plasmon-like and single-particle-like excitations Zhang17 ; Chapkin18 ; Senanayake19 ; Fitzgerald17 ; Ma15 ; You18 ; Schira19 . In those studies, the Coulomb interaction corresponding to part of the induced L field was considered. However, the T field was not precisely considered despite the light-induced excitations. As a phenomenological and semi-classical treatment of the microscopic nonlocal response involving the T field, the hydrodynamic model considering the “pressure” on the conduction electrons was applied to the relation between the polarization (or current) and the electric field Pitarke07 ; McMahan09 ; Raza11 ; Toscano12 ; Moreau13 ; Wang13 ; Mortensen14 ; Christensen14 . In such studies, aspects such as the plasmon peak shift, the presence of additional resonance above the plasma frequency, and the size effect were well discussed Raza11 ; Toscano12 ; Mortensen14 ; Christensen14 . The importance of the nonlocal effect for a dimer nanostructure with a gap much narrower than the light wavelength was noted Toscano12 ; Mortensen14 . However, depending on the sample structure or the situation, the hydrodynamic model requires an additional boundary condition or additional treatments from outside the microscopic model Svendsen20 . Importantly, coherent coupling between the plasmons and SPEs has not yet been fully discussed with this phenomenological model as a basis. For a small nanostructure, the plasmon spectrum becomes indistinguishable from the SPE spectrum Inaoka96 . Therefore, for the nonlocal effect in nanostructures with radiative fields, a microscopic formulation must be developed.

Motivated by the above-mentioned situation, in this study, we develop a quantum mechanical framework for mesoscopic plasmonics, focusing on the L and T components of EM fields and induced electronic responses in nanostructures that describe light-induced plasmons (collective excitations), SPEs (individual excitations), and their coherent coupling via both the L and T fields. An EM field induces charge polarization, charge current, spin fluctuation, and so on, and such electronic responses are described by a constitutive equation with a nonlocal susceptibility. The responses generate both L and T field components, and both field components must be comprehensively considered in a self-consistent treatment of the Maxwell’s and constitutive equations book:cho1 . Our formulation can be applied to arbitrary models of electronic states, such as those based on the Drude model Besteiro17 , jellium model with random phase approximation (RPA) Prodan01 , first-principles calculation Ma15 , or specific exotic materials, e.g., graphene Wang13 ; Koppens11 ; Lu19 . Such applications are expected to facilitate theoretical understanding of the interplay between the individual and collective excitations in nanoscale systems. Very recently, to understand the quantum effect in the light-matter interaction, a combination of the density functional theory and the macroscopic Maxwell’s equation has been discussed Svendsen21 . However, a discussion on the coupling of excitations caused by the L and T field components has not been developed yet. In this study, we develop an understanding of the coherent coupling between the individual and collective excitations in terms of the L and T fields and the charge and current densities. To examine the feasibility of our formulation, we numerically demonstrate the excitation spectrum for a rectangular nanorod. Then, we determine a non-negligible contribution of the T field component for a coherent coupling between the individual and collective excitations and energy dissipation. Further, our formulation provides a platform for the study of unrevealed physics concerning not only plasmons, but also excitons and/or other quasiparticles within nanostructures.

The remainder of this article is structured as follows. In Section II, we describe the basic quantities used in our framework and the nanostructure Hamiltonian. The Hamiltonian is separated into a nonperturbative term and an interaction with the field-induced potentials. This separation determines the treatment of the optical response and the electronic states of the nanostructures. In Section III, we introduce the four-vector representation and deduce the nonlocal susceptibility. In Section IV, we present a formal solution of the Maxwell’s equations with vector and scalar potentials, which is given by the Green’s function describing the propagated fields from the excited current and charge densities. In Section V, we explain the self-consistent treatment of the constitutive and Maxwell’s equations demonstrated in Secs. III and IV. In Section VI, we show the result of the numerical calculation of individual and collective excitations for a rectangular nanorod to demonstrate practicality of our formulation. In Section VII, we discuss possible applications and advantaged developments of our LT formulation for plasmonics and photonics. Finally, we present the summary in Section VIII.

In the present formulation, the fields are described by the vector and scalar potentials, $\bm{A}(\bm{r},t)$ and $\phi(\bm{r},t)$, respectively. In the quantum mechanics, they are essential rather than the electric and magnetic fields, $\bm{E}$ and $\bm{B}$. It has been experimentally proven by the Aharonov–Bohm experiment AB59 ; Tonomura86 . To divide the fields into the L and T components, we apply the Coulomb gauge, ${\rm div}\bm{A}(\bm{r},t)=0$, where the vector potential describes the T field only. In the constitutive equation, the electronic responses corresponding to the plasmon and SPE are described by the charge and current densities, $\rho(\bm{r},t)$ and $\bm{j}(\bm{r},t)$, respectively, which are directly coupled with $\phi$ and $\bm{A}$ in the interaction Hamiltonian. The nonlocal susceptibility is deduced from the Hamiltonian in accordance with quantum linear response theory Kubo57 ; book:cho1 , where the electron eigenstates determine the susceptibility. Therefore, the nanostructure sizes, shapes, and internal structures strongly affect the optical response via the electron states.

To clearly define the light–matter interaction, we first discuss the matter and interaction Hamiltonians in terms of the potentials. Within the nanostructures, mixing of the L and T field components and the densities occurs spontaneously. To formulate the L and T components explicitly, we introduce a four-vector representation for the fields $\bm{\mathcal{A}}=(\bm{A},-\phi/c)$ and densities $\bm{\mathcal{J}}=(\bm{j},c\rho)$. Hence, the nonlocal susceptibility is expressed as a $4\times 4$ matrix. This formulation enables us to exhibit the roles of the LT components of the fields and densities in enhancing the optical response and the interplay between the plasmons and SPEs.

The Hamiltonian for the electrons in the EM fields is

$$\hat{H}=\sum_{j}\left[\frac{\left\{\hat{\bm{p}}_{j}-e\int d\bm{r}\bm{A}(\bm{r},t)\delta(\bm{r}-\hat{\bm{r}}_{j})\right\}^{2}}{2m^{*}}-\varepsilon_{\rm F}\right]+e\sum_{j}\int d\bm{r}\phi(\bm{r},t)\delta(\bm{r}-\hat{\bm{r}}_{j})+\frac{1}{2}\sum_{i\neq j}\frac{e^{2}}{4\pi\varepsilon_{0}|\hat{\bm{r}}_{i}-\hat{\bm{r}}_{j}|},$$

(1)

where $e$ is the charge of electron and $\varepsilon_{0}$ is the dielectric constant in vacuum. $\varepsilon_{\rm F}$ is the Fermi energy and $m$ is the (effective) mass for the conduction electrons. The $\hat{\ }$ symbol indicates the electron operator. In the case of the Coulomb gauge, ${\rm div}\bm{A}=0$, and $\phi(\bm{r},t)=\phi_{\rm ncl}(\bm{r},t)+\phi_{\rm ext}(\bm{r},t)$ in the second term describes the L field due to the nuclei and external origins. The third term is the electron–electron (Coulomb) interaction, which is also longitudinal.

The first term is expanded as

	$\displaystyle\frac{\left\{\hat{\bm{p}}_{j}-e\int d\bm{r}\bm{A}(\bm{r},t)\delta(\bm{r}-\hat{\bm{r}}_{j})\right\}^{2}}{2m^{*}}$	$\displaystyle=$	$\displaystyle\frac{\hat{\bm{p}}_{j}^{2}}{2m^{*}}-\frac{e}{2m^{*}}\int d\bm{r}\left(\hat{\bm{p}}_{j}\delta(\bm{r}-\hat{\bm{r}}_{j})+\delta(\bm{r}-\hat{\bm{r}}_{j})\hat{\bm{p}}_{j}\right)\cdot\bm{A}(\bm{r},t)$		(2)
			$\displaystyle\hskip 85.35826pt+\frac{e^{2}}{2m^{*}}\int d\bm{r}\delta(\bm{r}-\hat{\bm{r}}_{j})\bm{A}(\hat{\bm{r}}_{j},t)\cdot\bm{A}(\bm{r},t).$

With the introduction of a current operator

$$\hat{\bm{I}}(\bm{r})\equiv\frac{e}{2m^{*}}\sum_{j}\left(\hat{\bm{p}}_{j}\delta(\bm{r}-\hat{\bm{r}}_{j})+\delta(\bm{r}-\hat{\bm{r}}_{j})\hat{\bm{p}}_{j}\right),$$

(3)

the second term on the right-hand side (r.h.s.) of Eq. (2) gives the light–matter interaction

$$\hat{H}_{IA}=-\int d\bm{r}\hat{\bm{I}}(\bm{r})\cdot\bm{A}(\bm{r},t).$$

(4)

Then, the Hamiltonian can be separated as follows:

$$\hat{H}=\hat{H}_{0}+\left(\hat{H}_{IA}+\hat{H}_{A^{2}}+\hat{H}_{\rm ext}\right),$$

(5)

with

$$\hat{H}_{0}=\sum_{j}\left[\frac{\hat{\bm{p}}_{j}^{2}}{2m^{*}}-\varepsilon_{\rm F}\right]+e\sum_{j}\int d\bm{r}\phi_{\rm ncl}(\bm{r},t)\delta(\bm{r}-\hat{\bm{r}}_{j})+\frac{1}{2}\frac{e^{2}}{4\pi\varepsilon_{0}}\sum_{i\neq j}\frac{1}{|\hat{\bm{r}}_{i}-\hat{\bm{r}}_{j}|}.$$

(6)

Note that the summations over $i,j$ include the electron spin degrees of freedom.

The second quantization of these Hamiltonians is

	$\displaystyle\hat{H}_{0}$	$\displaystyle=$	$\displaystyle\int d\bm{x}\hat{\Psi}^{\dagger}(\bm{x})\left(-\frac{\hbar^{2}\bm{\nabla}_{x}^{2}}{2m^{*}}-\varepsilon_{\rm F}\right)\hat{\Psi}(\bm{x})+e\int d\bm{x}\hat{\Psi}^{\dagger}(\bm{x})\int d\bm{r}\phi_{\rm ncl}(\bm{r},t)\delta(\bm{r}-\bm{x})\hat{\Psi}(\bm{x})$		(7)
			$\displaystyle\hskip 56.9055pt+\frac{1}{2}\frac{e^{2}}{4\pi\varepsilon_{0}}\int d\bm{x}\int d\bm{x}^{\prime}\hat{\Psi}^{\dagger}(\bm{x})\hat{\Psi}^{\dagger}(\bm{x}^{\prime})\frac{1}{|\bm{x}-\bm{x}^{\prime}|}\hat{\Psi}(\bm{x}^{\prime})\hat{\Psi}(\bm{x})$
		$\displaystyle=$	$\displaystyle\sum_{n,n^{\prime}}\int d\bm{x}\psi_{n^{\prime}}^{*}(\bm{x})\left(-\frac{\hbar^{2}\bm{\nabla}_{x}^{2}}{2m^{*}}-\varepsilon_{\rm F}\right)\psi_{n}(\bm{x})\hat{a}_{n^{\prime}}^{\dagger}\hat{a}_{n}+\sum_{n,n^{\prime}}\int d\bm{x}\rho_{n^{\prime}n}(\bm{x})\phi_{\rm ncl}(\bm{x},t)\hat{a}_{n^{\prime}}^{\dagger}\hat{a}_{n}$
			$\displaystyle\hskip 113.81102pt+\frac{1}{2}\sum_{n,n^{\prime},m,m^{\prime}}\int d\bm{x}\int d\bm{x}^{\prime}\frac{\rho_{n^{\prime}n}(\bm{x})\rho_{m^{\prime}m}(\bm{x}^{\prime})}{4\pi\varepsilon_{0}|\bm{x}-\bm{x}^{\prime}|}\hat{a}_{n^{\prime}}^{\dagger}\hat{a}_{m^{\prime}}^{\dagger}\hat{a}_{m}\hat{a}_{n},$
	$\displaystyle\hat{H}_{IA}$	$\displaystyle=$	$\displaystyle-\sum_{n,n^{\prime}}\int d\bm{x}\bm{I}_{n^{\prime}n}(\bm{x})\cdot\bm{A}(\bm{x},t)\hat{a}_{n^{\prime}}^{\dagger}\hat{a}_{n},$		(8)
	$\displaystyle\hat{H}_{A^{2}}$	$\displaystyle=$	$\displaystyle\frac{e^{2}}{2m^{*}}\sum_{n,n^{\prime}}\int d\bm{x}\psi_{n^{\prime}}^{*}(\bm{x})\bm{A}(\bm{x},t)\cdot\bm{A}(\bm{x},t)\psi_{n}(\bm{x})\hat{a}_{n^{\prime}}^{\dagger}\hat{a}_{n},$		(9)
	$\displaystyle{\hat{H}_{\rm ext}}$	$\displaystyle=$	$\displaystyle{\sum_{n,n^{\prime}}\int d\bm{x}\rho_{n^{\prime}n}(\bm{x})\phi_{\rm ext}(\bm{x},t)\hat{a}_{n^{\prime}}^{\dagger}\hat{a}_{n},}$		(10)

where $\hat{\Psi}(\bm{x})=\sum_{n}\psi_{n}(\bm{x})\hat{a}_{n}$ is the field operator. Note that the $n$ state includes the spin degrees of freedom. Here,

$$\rho_{n^{\prime}n}(\bm{x})=e\psi_{n^{\prime}}^{*}(\bm{x})\psi_{n}(\bm{x})$$

(11)

indicates an element of charge density and

$$\bm{I}_{n^{\prime}n}(\bm{x})=-\frac{e\hbar}{2im^{*}}\Big{[}\left\{\bm{\nabla}_{x}\psi_{n^{\prime}}^{*}(\bm{x})\right\}\psi_{n}(\bm{x})-\psi_{n^{\prime}}^{*}(\bm{x})\left\{\bm{\nabla}_{x}\psi_{n}(\bm{x})\right\}\Big{]}$$

(12)

is derived from the current operator in Eq. (3). Note that $\rho_{n^{\prime}n}(\bm{x})$ has to follow the continuous relation with the charge current density. Further, the single-electron wave function follows the time-dependent Schrödinger equation,

$$i\hbar\frac{\partial\psi(\bm{x},t)}{\partial t}=\left[-\frac{\hbar^{2}}{2m^{*}}\left(\bm{\nabla}_{x}-\frac{ie}{\hbar}\bm{A}(\bm{x},t)\right)^{2}+e\phi_{\rm ncl}(\bm{x},t)\right]\psi(\bm{x},t).$$

(13)

Then, $\rho_{n^{\prime}n}(\bm{x})$ satisfies

$$\frac{\partial\rho_{n^{\prime}n}(\bm{x},t)}{\partial t}=-\bm{\nabla}_{x}\cdot\bm{j}_{n^{\prime}n}(\bm{x},t),$$

(14)

where

$$\bm{j}_{n^{\prime}n}(\bm{x},t)=\bm{I}_{n^{\prime}n}(\bm{x})-\frac{e^{2}}{m^{*}}\psi_{n^{\prime}}^{*}(\bm{x})\bm{A}(\bm{x},t)\psi_{n}(\bm{x})$$

(15)

is the modified current density in the EM field. The second term contributes to the $O(\bm{A}^{2})$ term. The light–matter interaction component with the vector potential in Eq. (5) is approximately

	$\displaystyle\hat{H}_{IA}+\hat{H}_{A^{2}}$	$\displaystyle\simeq$	$\displaystyle-\sum_{n,n^{\prime}}\int d\bm{x}\bm{j}_{n^{\prime}n}(\bm{x},t)\cdot\bm{A}(\bm{x},t)\hat{a}_{n^{\prime}}^{\dagger}\hat{a}_{n}$		(16)
		$\displaystyle=$	$\displaystyle-\int d\bm{x}\hat{\bm{j}}(\bm{x},t)\cdot\bm{A}(\bm{x},t)\equiv\hat{H}_{\rm int}.$

The second and third terms on the r.h.s. of Eq. (7) can be interpreted as the L component of the electric fields due to internal sources. In a following formulation of the constitutive equation, the T and L fields are explicitly distinguished by the vector and scalar potentials, respectively, with the Coulomb gauge. To describe the susceptibility for the vector and scalar potentials, we treat the many-electron system with $\hat{H}_{0}$ on a one-electron basis and rewrite the second and third terms as

	$\displaystyle\int d\bm{x}\sum_{n,n^{\prime}}\hat{a}_{n^{\prime}}^{\dagger}\rho_{n^{\prime}n}(\bm{x},t)\left[\phi_{\rm ncl}(\bm{x},t)+\frac{e}{4\pi\varepsilon_{0}}\int d\bm{x}^{\prime}\frac{\hat{\Psi}^{\dagger}(\bm{x}^{\prime},t)\hat{\Psi}(\bm{x}^{\prime},t)}{|\bm{x}-\bm{x}^{\prime}|}\right]\hat{a}_{n}$
	$\displaystyle=\int d\bm{x}\sum_{n,n^{\prime}}\hat{a}_{n^{\prime}}^{\dagger}\rho_{n^{\prime}n}(\bm{x},t)\left[\hat{\phi}_{\rm{mat}}(\bm{x},t)\right]\hat{a}_{n},$		(17)

with an inherent scalar potential operator for matter,

	$\displaystyle\hat{\phi}_{\rm{mat}}(\bm{x},t)$	$\displaystyle=$	$\displaystyle\phi_{\rm ncl}(\bm{x},t)+\hat{\phi}_{\rm e-e}(\bm{x},t),$		(18)
	$\displaystyle{\hat{\phi}_{\rm e-e}(\bm{x},t)}$	$\displaystyle\equiv$	$\displaystyle{\frac{e}{4\pi\varepsilon_{0}}\int d\bm{x}^{\prime}\frac{\hat{\Psi}^{\dagger}(\bm{x}^{\prime},t)\hat{\Psi}(\bm{x}^{\prime},t)}{|\bm{x}-\bm{x}^{\prime}|}.}$		(19)

Note that, in the absence of $\bm{A}$ and $\phi_{\rm ext}$, a material is electrically neutral, and $\left\langle\hat{\phi}_{\rm{mat}}(\bm{x},t)\right\rangle_{0}=0$. Here, we take the basis for the static state without external fields and consider $\hat{H}_{0}$ in Eq. (7) as a non-perturbative Hamiltonian. When one applies the electromagnetic field of $\bm{A}(\bm{x},t)$ and $\phi_{\rm ext}(\bm{x},t)$ to the material, a polarized charge is induced. It gives an additional interaction between the induced polarized charges,

	$\displaystyle\hat{H}_{\rm p-p}$	$\displaystyle\simeq$	$\displaystyle\int d\bm{x}\frac{\hat{\rho}(\bm{x},t)-\rho_{0}(\bm{x})}{4\pi\varepsilon_{0}}\int d\bm{x}^{\prime}\frac{\hat{\rho}(\bm{x}^{\prime},t)-\rho_{0}(\bm{x}^{\prime})}{|\bm{x}-\bm{x}^{\prime}|}$		(20)
		$\displaystyle=$	$\displaystyle{\int d\bm{x}\delta\hat{\rho}(\bm{x},t)\hat{\phi}_{\rm pol}(\bm{x},t).}$

Here, $\rho_{0}(\bm{x})=\left\langle\hat{\rho}(\bm{x},t)\right\rangle_{0}$ is the average charge density from the electrons [$\hat{\rho}(\bm{x},t)\equiv e\hat{\Psi}^{\dagger}(\bm{x},t)\hat{\Psi}(\bm{x},t)$]. Thus, $\delta\hat{\rho}(\bm{x},t)=\hat{\rho}(\bm{x},t)-\rho_{0}(\bm{x})$ indicates a deviation from the static distribution. The scalar potential due to the induced polarized charge is

$${\hat{\phi}_{\rm pol}(\bm{x},t)=\int d\bm{x}^{\prime}\frac{\delta\hat{\rho}(\bm{x}^{\prime},t)}{4\pi\varepsilon_{0}|\bm{x}-\bm{x}^{\prime}|}.}$$

(21)

Therefore, the L field consists of three elements: $\hat{\phi}_{\rm{mat}}$, $\hat{\phi}_{\rm{pol}}$, and $\phi_{\rm{ext}}$, which are associated with the electron–electron interaction in the static case, the induced charge densities, and the external field, respectively.

In this study, we divide the total Hamiltonian as $\hat{H}=\hat{H}_{0}+\hat{H}_{\rm ind}$ with

$$\hat{H}_{\rm ind}=\hat{H}_{\rm int}+\hat{H}_{\rm ext}+\hat{H}_{\rm p-p}.$$

(22)

$\hat{H}_{0}$ describes the static system without any field irradiation, whereas all induced effects by the external electric and magnetic fields are in $\hat{H}_{\rm ind}$. Then, the incident and induced components are included in the perturbative term. For the scalar potential of the induced polarized charge, we apply the mean field approximation, $\hat{\Psi}^{\dagger}(\bm{x},t)\hat{\phi}_{\rm pol}(\bm{x},t)\hat{\Psi}(\bm{x},t)\simeq\hat{\Psi}^{\dagger}(\bm{x},t)\hat{\Psi}(\bm{x},t)\phi_{\rm pol}(\bm{x},t)$. For the external potential, $\int d\bm{x}\rho_{0}(\bm{x})\phi_{\rm ext}(\bm{x},t)$ corresponds to the system energy shifts without excitation, as it is constant. Subtracting this term from the Hamiltonian, we obtain

	$\displaystyle\hat{H}_{\rm ind}^{\prime}$	$\displaystyle=$	$\displaystyle\hat{H}_{\rm ind}-\left(\int d\bm{x}\rho_{0}(\bm{x})\phi_{\rm ext}(\bm{x},t)\right)$		(23)
		$\displaystyle\simeq$	$\displaystyle-\int d\bm{x}\left[\hat{\bm{j}}(\bm{x},t)\cdot\bm{A}(\bm{x},t)-\delta\hat{\rho}(\bm{x},t)\Big{(}\phi_{\rm ext}(\bm{x},t)+\phi_{\rm pol}(\bm{x},t)\Big{)}\right],$

which describes the interaction between the matter and applied and induced EM fields for both the L and T components. The L field is included in both the constitutive and Maxwell’s equations. However, we do not need to take care of a double count of the electron–electron interaction (see Appendix A).

The susceptibility is obtained from the linear response theory Kubo57 . The susceptibility, field-induced charge density, and induced current density formulations depend on the treatment of the non-perturbative and perturbative Hamiltonians. Note that $\hat{H}_{0}$ includes neither a net charge nor current. Therefore, for the formulation of the susceptibility based on $\hat{H}=\hat{H}_{0}+\hat{H}_{\rm ind}^{\prime}$ is evaluated by the static states $|\mu\rangle$ for $\hat{H}_{0}$. As this susceptibility is irrelevant to the induced fields, we can analyze the T and L components of the field-induced charge and current in terms of $\bm{A}$ (T field) and $\phi$ (L field). A self-consistent relation between the induced charge and $\phi_{\rm{pol}}$ can describe a plasmon (collective excitation) within the RPA Prodan01 , as discussed in the following section.

We can omit $\delta$ of $\delta\hat{\rho}(\bm{x})$ in Eq. (23) for simplicity. In the following discussion, the charge density $\hat{\rho}(\bm{x})$ indicates the induced charge only. This abbreviation does not change the continuous relation, as $\partial\rho_{0}/\partial t=0$.

The light–matter interaction $\hat{H}_{\rm ind}^{\prime}(t)$ in Eq. (23) gives a susceptibility to the external fields, $\bm{A}(\bm{x},t)$ and $\phi_{\rm ind}(\bm{x},t)=\phi_{\rm ext}(\bm{x},t)+\phi_{\rm pol}(\bm{x},t)$ in linear response theory Kubo57 . Here, $\phi_{\rm ind}$ indicates a scalar potential caused by the field irradiation. In the following, let us omit the subscript “${\rm ind}$” on $\hat{H}_{\rm ind}^{\prime}$ and $\phi_{\rm ind}$ for simplicity. For the interaction Hamiltonian $\hat{H}^{\prime}(t)$, let us introduce the following four-vector representation:

	$\displaystyle{\bm{\mathcal{A}}(\bm{x},t)}$	$\displaystyle=$	$\displaystyle{\left(\begin{matrix}\bm{A}(\bm{x},t)\\ -\frac{1}{c}\phi(\bm{x},t)\end{matrix}\right),}$		(24)
	$\displaystyle\hat{\bm{\mathcal{J}}}(\bm{x},t)$	$\displaystyle=$	$\displaystyle\left(\begin{matrix}\hat{\bm{j}}(\bm{x},t)\\ c\hat{\rho}(\bm{x},t)\end{matrix}\right),$		(25)
	$\displaystyle\hat{H}^{\prime}(t)$	$\displaystyle=$	$\displaystyle-\int d\bm{x}\hat{\bm{\mathcal{J}}}(\bm{x},t)\cdot\bm{\mathcal{A}}(\bm{x},t).$		(26)

The time-dependence of $\hat{H}^{\prime}(t)$ in Eq. (26) is attributed to the external field $\bm{\mathcal{A}}(\bm{x},t)$ only, and $c$ is the light velocity in vacuum. Note that we define the sign for the scalar potential component as negative in Eq. (24) for a simpler formulation.

We assume a monochromatic field, $\bm{\mathcal{A}}(\bm{x},t)=\bm{\mathcal{A}}(\bm{x};\omega)e^{-i\omega t}$. The statistical average of the four-vector current $\hat{\bm{\mathcal{J}}}(\bm{x})$ gives the susceptibility by the field $\bm{\mathcal{A}}(\bm{x};\omega)$:

	$\displaystyle\left\langle\hat{\bm{\mathcal{J}}}(\bm{x})\right\rangle(t)$	$\displaystyle=$	$\displaystyle{\bm{\mathcal{J}}(\bm{x};\omega)e^{-i\omega t}}$		(27)
		$\displaystyle=$	$\displaystyle{\bm{\mathcal{J}}_{0}(\bm{x};\omega)e^{-i\omega t}}+\int d\bm{x}^{\prime}\bar{\mathcal{X}}(\bm{x},\bm{x}^{\prime};\omega)\cdot\bm{\mathcal{A}}(\bm{x}^{\prime};\omega)e^{-i\omega t},$

with a nonlocal susceptibility book:cho1

	$\displaystyle\bar{\mathcal{X}}(\bm{x},\bm{x}^{\prime};\omega)$	$\displaystyle\equiv$	$\displaystyle\sum_{\mu,\nu}\left[f_{\nu\mu}\bm{\mathcal{J}}_{\mu\nu}(\bm{x})\left(\bm{\mathcal{J}}_{\nu\mu}(\bm{x}^{\prime})\right)^{\rm t}+h_{\nu\mu}\bm{\mathcal{J}}_{\nu\mu}(\bm{x})\left(\bm{\mathcal{J}}_{\mu\nu}(\bm{x}^{\prime})\right)^{\rm t}\right]$		(28)
		$\displaystyle=$	$\displaystyle\left(\begin{matrix}\bar{\chi}_{jj}(\bm{x},\bm{x}^{\prime};\omega)&\bm{\chi}_{j\rho}(\bm{x},\bm{x}^{\prime};\omega)\\ \bm{\chi}_{\rho j}^{\rm\ t}(\bm{x},\bm{x}^{\prime};\omega)&\chi_{\rho\rho}(\bm{x},\bm{x}^{\prime};\omega)\end{matrix}\right),$

and

	$\displaystyle\bm{\mathcal{J}}_{\nu\mu}(\bm{x})$	$\displaystyle=$	$\displaystyle\left\langle\nu\left|\hat{\bm{\mathcal{J}}}(\bm{x})\right|\mu\right\rangle,$		(29)
	$\displaystyle f_{\nu\mu}$	$\displaystyle=$	$\displaystyle\frac{\rho_{0,\mu}}{\hbar\omega_{\nu\mu}-\hbar\omega-i\gamma},$		(30)
	$\displaystyle h_{\nu\mu}$	$\displaystyle=$	$\displaystyle\frac{\rho_{0,\mu}}{\hbar\omega_{\nu\mu}+\hbar\omega+i\gamma}.$		(31)

The nonlocal susceptibility is a $4\times 4$ matrix in the four-vector space and $|\mu\rangle$ represents the electronic eigenstates of $\hat{H}_{0}$. Thus, the susceptibility is affected by the nanostructure geometry via $|\mu\rangle$. If the electronic system has translational symmetry, e.g., in the case of a bulk metal, the nonlocality in the susceptibility is given only by a relative position, $\bar{\mathcal{X}}(\bm{x}-\bm{x}^{\prime};\omega)$. Further, $\hbar\omega_{\nu\mu}=\varepsilon_{\nu}-\varepsilon_{\mu}$ is the energy difference between the eigenenergies for $\hat{H}_{0}$, $\gamma$ is an infinitesimal positive value for the causality, and $\rho_{0,\mu}$ in the numerators is an element of the density matrix at equilibrium, where

$$\rho_{0,\mu}=\langle\mu|\hat{\rho}_{0}|\mu\rangle=\frac{1}{Z_{0}}\langle\mu|e^{-\beta\hat{H}_{0}}|\mu\rangle$$

(32)

with the partition function $Z_{0}={\rm Tr}\{e^{-\beta\hat{H}_{0}}\}$. The first term in the r.h.s. of Eq. (27) is

$$\bm{\mathcal{J}}_{0}(\bm{x};\omega)=\left(\begin{matrix}\langle\hat{\bm{j}}(\bm{x};\omega)\rangle_{0}\\ c\left(\langle\hat{\rho}(\bm{x})\rangle_{0}-\rho_{0}(\bm{x})\right)\end{matrix}\right)=\left(\begin{matrix}-\frac{e}{m^{*}}\rho_{0}(\bm{x})\bm{A}(\bm{x};\omega)\\ 0\end{matrix}\right)$$

(33)

at equilibrium. The elements of susceptibility in Eq. (28) are related to each other and satisfy the continuous relation $\bm{\nabla}_{x}\cdot\left\langle\hat{\bm{j}}(\bm{x})\right\rangle-i\omega\left\langle\hat{\rho}(\bm{x})\right\rangle=0$, where

	$\displaystyle\bm{\nabla}_{x}\cdot\bar{\chi}_{jj}(\bm{x},\bm{x}^{\prime};\omega)$	$\displaystyle=$	$\displaystyle i(\omega/c)\bm{\chi}_{\rho j}^{\rm\ t}(\bm{x},\bm{x}^{\prime};\omega),$		(34)
	$\displaystyle\bm{\nabla}_{x}\cdot\bm{\chi}_{j\rho}(\bm{x},\bm{x}^{\prime};\omega)$	$\displaystyle=$	$\displaystyle i(\omega/c)\chi_{\rho\rho}(\bm{x},\bm{x}^{\prime};\omega).$		(35)

In the linear response, the $\bm{\mathcal{A}}$-dependence of the susceptibility $\bar{\mathcal{X}}(\bm{x},\bm{x}^{\prime};\omega)$ cannot be discussed. To evaluate the effect of the $\bm{A}(\bm{x},t)$ term in Eq. (28), the higher-order terms must be considered.

Taking the Coulomb gauge, the vector and scalar potentials describe only the T and L components, respectively. For the Maxwell’s equations, we consider the potentials $\bm{\mathcal{A}}(\bm{x},t)$ induced by the densities $\bm{\mathcal{J}}(\bm{x},t)$ via the Green’s function.

The Maxwell’s equations for the potentials, $\phi(\bm{x},t)$ and $\bm{A}(\bm{x},t)$, are expressed as

	$\displaystyle\bm{\nabla}^{2}\bm{A}(\bm{x},t)-\frac{1}{c^{2}}\frac{\partial^{2}\bm{A}(\bm{x},t)}{\partial t^{2}}-\frac{1}{c^{2}}\frac{\partial}{\partial t}\bm{\nabla}\phi(\bm{x},t)$	$\displaystyle=$	$\displaystyle-\mu_{0}\bm{j}(\bm{x},t),$		(36)
	$\displaystyle\bm{\nabla}^{2}\phi(\bm{x},t)$	$\displaystyle=$	$\displaystyle-\frac{\rho(\bm{x},t)}{\varepsilon_{0}}.$		(37)

Note that the L component of the field is relevant only to $\rho(\bm{x},t)$, whereas $\bm{j}(\bm{x},t)$ generates both of the L and T components. Another description may also be considered (see Appendix B). Equations (36) and (37) are unified in the four-vector representation

$$\left(\begin{matrix}\bm{\nabla}^{2}-\frac{\partial^{2}}{\partial(ct)^{2}}&{\frac{\partial}{\partial(ct)}\bm{\nabla}}\\ 0&{-\bm{\nabla}^{2}}\end{matrix}\right)\bm{\mathcal{A}}(\bm{x},t)=-\mu_{0}\bm{\mathcal{J}}(\bm{x},t),$$

(38)

with the statistical average of the densities, $\bm{\mathcal{J}}(\bm{x},t)=\left\langle\hat{\bm{\mathcal{J}}}(\bm{x})\right\rangle(t)$. For the Fourier component, $\bm{\mathcal{O}}(\bm{x},t)=\bm{\mathcal{O}}(\bm{x};\omega)e^{-i\omega t}$ ($\bm{\mathcal{O}}=\bm{\mathcal{A}},\bm{\mathcal{J}}$), the Maxwell’s equation is written as

$$\bar{\mathcal{D}}(\bm{x};\omega)\bm{\mathcal{A}}(\bm{x};\omega)=-\mu_{0}\bm{\mathcal{J}}(\bm{x};\omega),$$

(39)

with the $4\times 4$ differential matrix

$$\bar{\mathcal{D}}(\bm{x};\omega)\equiv\left(\begin{matrix}\bm{\nabla}^{2}+\frac{\omega^{2}}{c^{2}}&{-i\frac{\omega}{c}\bm{\nabla}}\\ 0&{-\bm{\nabla}^{2}}\end{matrix}\right).$$

(40)

From Eqs. (27) and (39), we can construct a self-consistent relation between $\bm{\mathcal{J}}$ and $\bm{\mathcal{A}}$. For the self-consistent equation, let us obtain the formal solution of the Maxwell’s equation (39) using the Green’s function

$$\bm{\mathcal{A}}(\bm{x};\omega)=\bm{\mathcal{A}}_{0}(\bm{x};\omega)-\mu_{0}\int d\bm{x}^{\prime}\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)\bm{\mathcal{J}}(\bm{x}^{\prime};\omega).$$

(41)

The first term in Eq. (41) corresponds to an “incident” field satisfying

$$\bar{\mathcal{D}}(\bm{x};\omega)\bm{\mathcal{A}}_{0}(\bm{x};\omega)=0.$$

(42)

The Green’s function is given as

$$\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)=-\frac{1}{4\pi}\left(\begin{matrix}\frac{\ e^{i\frac{\omega}{c}|\bm{x}-\bm{x}^{\prime}|}}{|\bm{x}-\bm{x}^{\prime}|}\bar{1}\ &\frac{-i\omega}{4\pi c}\int d\bm{x}^{\prime\prime}\frac{e^{i\frac{\omega}{c}|\bm{x}-\bm{x}^{\prime\prime}|}}{|\bm{x}-\bm{x}^{\prime\prime}|}\frac{\bm{x}^{\prime\prime}-\bm{x}^{\prime}}{|\bm{x}^{\prime\prime}-\bm{x}^{\prime}|^{3}}\\ 0&{-\frac{1}{|\bm{x}-\bm{x}^{\prime}|}}\end{matrix}\right).$$

(43)

A detailed derivation of $\bar{\mathcal{G}}$ in Eq. (43) is presented in Appendix C. This Green’s function provides a scheme for field (potential) generation due to the electronic excitations (current and charge densities).

In previous sections, we derived the constitutive equation (27) with the nonlocal susceptibility (28) and the solution of the Maxwell’s equations (41) with Green’s function (43) in a four-vector representation. These expressions form a self-consistent relation.