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

Tomohiro Yokoyama<EMAIL_ADDRESS>Department of
Material Engineering Science, Graduate School of Engineering Science, Osaka
University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan  Masayuki Iio
Department of Material Engineering Science, Graduate School of Engineering
Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan
Department of Physics and Electronics, Graduate School of Engineering, Osaka
Prefecture University, 1-1 Gakuen-cho, Naka-ku, Sakai, Osaka 599-8531, Japan
Takashi Kinoshita  Department of Physics and Electronics, Graduate School of
Engineering, Osaka Prefecture University, 1-1 Gakuen-cho, Naka-ku, Sakai,
Osaka 599-8531, Japan  Takeshi Inaoka  Department of Physics and Earth
Sciences, Faculty of Science, University of the Ryukyus, 1 Senbaru, Nishihara-
cho, Okinawa 903-0213, Japan  Hajime Ishihara  Department of Material
Engineering Science, Graduate School of Engineering Science, Osaka University,
1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan  Department of Physics and
Electronics, Graduate School of Engineering, Osaka Prefecture University, 1-1
Gakuen-cho, Naka-ku, Sakai, Osaka 599-8531, Japan  Center for Quantum
Information and Quantum Biology, Osaka University, 1-3 Machikaneyama,
Toyonaka, Osaka 560-8531, Japan

###### Abstract

A plasmon is a collective excitation of electrons due to the Coulomb
interaction. Both plasmons and single-particle excitations (SPEs) are
eigenstates of bulk metallic systems and they are orthogonal to each other.
However, in non-translationally symmetric systems such as nanostructures,
plasmons and SPEs coherently interact. It has been well discussed that the
plasmons and SPEs, respectively, can couple with transverse (T) electric field
in such systems, and also that they are coupled with each other via
longitudinal (L) field. However, there has been a missing link in the previous
studies: the coherent coupling between the plasmons and SPEs mediated by the T
field. Herein, we develop a theoretical framework to describe the self-
consistent relationship between plasmons and SPEs through both the L and T
fields. The excitations are described in terms of the charge and current
densities in a constitutive equation with a nonlocal susceptibility, where the
densities include the L and T components. The electromagnetic fields
originating from the densities are described in terms of the Green’s function
in the Maxwell equations. The T field is generated from both densities,
whereas the L component is attributed to the charge density only. We introduce
a four-vector representation incorporating the vector and scalar potentials in
the Coulomb gauge, in which the T and L fields are separated explicitly. The
eigenvalues of the matrix for the self-consistent equations appear as the
poles of the system excitations. The developed formulation enables to approach
unknown mechanisms for enhancement of the coherent coupling between plasmons
and the hot carriers generated by radiative fields.

##  I INTRODUCTION

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 plasmons and individual SPEs. Recently, Ma et al. explicitly discussed
the interplay between 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 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 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, SPEs, 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 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 collective and
individual excitations in nanoscale systems. 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 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 equations demonstrated in
Secs. III and IV. In Section VI, we discuss possible applications of our LT
formulation for plasmonics and photonics. Finally, we present the summary in
Section VII.

##  II Basic quantities and Hamiltonian of electrons with electromagnetic
field

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}\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}+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 $m$ and $e$ are the electron mass and charge, respectively, and
$\varepsilon_{0}$ is the dielectric constant of the vacuum. 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 latter is set to have a negative sign for convenience of
the following formulation. Further, 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}\frac{\hat{\bm{p}}_{j}^{2}}{2m}+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}\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}\right)\psi_{n}(\bm{x})\hat{a}_{n}^{\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})=\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})\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})\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.

The total Hamiltonian can be read according to two descriptions, as

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\hat{H}_{0}+\hat{H}_{\rm
p-p}+\hat{H}_{\rm int}+\hat{H}_{\rm ext}=\hat{H}_{0}^{\prime}+\hat{H}_{\rm
int}^{\prime},$ (22) $\displaystyle=$ $\displaystyle\hat{H}_{0}+\hat{H}_{\rm
int}+\hat{H}_{\rm ext}+\hat{H}_{\rm p-p}=\hat{H}_{0}+\hat{H}_{\rm ind},$ (23)

depending on the L-field treatment. For the upper description,

$\displaystyle\hat{H}_{\rm int}^{\prime}$ $\displaystyle\equiv$
$\displaystyle\hat{H}_{\rm int}+\hat{H}_{\rm ext}$ (24) $\displaystyle=$
$\displaystyle-\int
d\bm{x}\left[\hat{\bm{j}}(\bm{x})\cdot\bm{A}(\bm{x},t)+\hat{\rho}(\bm{x})\phi_{\rm
ext}(\bm{x},t)\right].$

Here, the induced charge density interacts with only an external (applied) L
field, whereas both the incident and induced T fields are coupled with the
charge current density. In this case, the electron–electron interaction is
fully included in the non-perturbative Hamiltonian,
$\hat{H}_{0}^{\prime}\equiv\hat{H}_{0}+\hat{H}_{\rm p-p}$. Therefore, one must
prepare a quite complex and time-dependent basis. This approach is reasonable
if one does not consider the self-consistent relation between the constitutive
and Maxwell equations discussed in the following sections. In the description
given by Eqs. (22) and (24), the constitutive equation includes
$\hat{\phi}_{\rm{mat}}$ and $\hat{\phi}_{\rm{pol}}$, and the interaction
between the induced polarization charges is also considered in the Maxwell
equations. Double counting of the electron–electron interaction in the self-
consistent treatment must also be handled. Therefore, this approach is not
useful because of the basis complexity and the double count problem.

For the latter description of Eq. (23), all effects are due to the presence of
external electric and magnetic fields Therefore, 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 shift without excitation, as it is simply constant. Subtracting this
term from the Hamiltonian, we obtain

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

which describes the interaction between the matter and the applied and induced
EM fields for both the L and T components. Note again that, although
$\phi_{\rm pol}$ is attributed to the Coulomb interaction between the internal
charges of the matter, it is present only when the external fields are
applied. Moreover, because of the electric neutrality of the matter,
$\hat{\phi}_{\rm{mat}}$ is not incorporated in the Maxwell equations.
Therefore, no treatment of the double counting of the electron–electron
interaction is required, as the induced Coulomb interaction is irrelevant to
the susceptibility in the constitutive equation and is considered only in the
Maxwell equations.

non-perturb. H | perturb. H | basis $|\mu\rangle$ | Coulomb int. | L field | T field | susceptibility
---|---|---|---|---|---|---
$\hat{H}_{0}+\hat{H}_{\rm{p-p}}$ | $\hat{H}_{\rm{int}}^{\prime}$ | time-dependent | fully in $|\mu\rangle$ | $\phi_{\rm{ext}}$ | $\bm{A}$ | time-dep. via $|\mu\rangle$
$\hat{H}_{0}$ | $\hat{H}_{\rm{ind}}$ | time-independent | partially in $|\mu\rangle$ | $\phi_{\rm{pol}}+\phi_{\rm{ext}}$ | $\bm{A}$ | static
Table 1:  Treatment of Hamiltonian in Eqs. (22) and (23). Here,
$\phi_{\rm{ext}}$ describes an external (applied) L field, $\phi_{\rm{pol}}$
is an induced L field, and $\bm{A}$ consists of the T field of both the
applied and induced components.

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 based on Eqs. (23) and (25), the
susceptibility is evaluated based on the static states $|\mu\rangle$ only,
without the irradiated field. 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. On the other hand, for Eqs. (22) and (24),
a part of L field and time-dependent term $\phi_{\rm pol}$ is included in the
non-perturbative Hamiltonian. Thus, the evaluation of susceptibility by the
time-dependent basis becomes complicated. The induced Coulomb interaction
contributes to the induced charge and current densities via the
susceptibility. This approach is not useful for analysis of the relation
between the L and T components of the induced fields and densities.

The above discussion and classification are summarized in TABLE I and Figure
1. Hereafter, we adopt the representation based on Eqs. (23) and (25). Note
that we can omit the $\delta$ of $\delta\hat{\rho}(\bm{x})$ in Eq. (25) for
simplicity. In the following discussion, the charge density
$\hat{\rho}(\bm{x})$ indicates the induced polarized charge only. This
abbreviation does not change the continuous relation, as
$\partial\rho_{0}/\partial t=0$.

Figure 1:  Schematic summary of two treatments of the Hamiltonian. In upper
case, the electronic state of non-perturbative Hamiltonian has induced
polarized charges. The Coulomb interaction of the induced charges is included
in both the constitutive and Maxwell’s equations. Then, one has to take care
of the double count of Coulomb interaction. Lower case is our scheme, where
the Coulomb interaction of the induced charges is taken into account only in
the Maxwell’s equations. Therefore, one does not need to consider the double
count.

##  III Four-vector representation for nonlocal susceptibility

The light–matter interaction $\hat{H}_{\rm ind}(t)$ in Eq. (25) 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
of $\phi_{\rm ind}$. For the interaction Hamiltonian $\hat{H}_{\rm ind}(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),$ (26)
$\displaystyle\hat{\bm{\mathcal{J}}}(\bm{x})$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\hat{\bm{j}}(\bm{x})\\\
c\hat{\rho}(\bm{x})\end{matrix}\right),$ (27) $\displaystyle\hat{H}_{\rm
ind}(t)$ $\displaystyle=$ $\displaystyle-\int
d\bm{x}\hat{\bm{\mathcal{J}}}(\bm{x})\cdot\bm{\mathcal{A}}(\bm{x},t).$ (28)

The time-dependence of $\hat{H}_{\rm ind}(t)$ in Eq. (28) is attributed to the
external field $\bm{\mathcal{A}}(\bm{x},t)$ only, and $c$ is the light
velocity in vacuum. 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)$:

$\left\langle\hat{\bm{\mathcal{J}}}(\bm{x})\right\rangle(t)=\left\langle\hat{\bm{\mathcal{J}}}(\bm{x})\right\rangle_{0}+\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},$ (29)

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]$ (30) $\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,$
(31) $\displaystyle f_{\nu\mu}$ $\displaystyle=$
$\displaystyle\frac{\rho_{0,\mu}}{\hbar\omega_{\nu\mu}-\hbar\omega-i\gamma},$
(32) $\displaystyle h_{\nu\mu}$ $\displaystyle=$
$\displaystyle\frac{\rho_{0,\mu}}{\hbar\omega_{\nu\mu}+\hbar\omega+i\gamma}.$
(33)

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 is indifferent, 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$
(34)

with the partition function $Z_{0}={\rm Tr}\\{e^{-\beta\hat{H}_{0}}\\}$. At
equilibrium in the absence of an external field,
$\left\langle\hat{\bm{\mathcal{J}}}(\bm{x})\right\rangle_{0}=0$. The elements
of susceptibility in Eq. (30) 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),$ (35)
$\displaystyle\bm{\nabla}_{x}\cdot\bm{\chi}_{j\rho}(\bm{x},\bm{x}^{\prime};\omega)$
$\displaystyle=$ $\displaystyle
i(\omega/c)\bar{\chi}_{\rho\rho}(\bm{x},\bm{x}^{\prime};\omega).$ (36)

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. (15),
the higher-order terms must be considered.

##  IV Maxwell equations for Coulomb gauge

Taking the Coulomb gauge, the vector and scalar potentials describe only the T
and L components, respectively. For the Maxwell 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 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),$ (37)
$\displaystyle\bm{\nabla}^{2}\phi(\bm{x},t)$ $\displaystyle=$
$\displaystyle-\frac{\rho(\bm{x},t)}{\varepsilon_{0}}.$ (38)

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 A). Equations (37) and (38)
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),$
(39)

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 equation
is written as

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

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).$ (41)

From Eqs. (29) and (40), 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 equation (40) 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).$
(42)

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

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

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^{ik|\bm{x}-\bm{x}^{\prime}|}}{|\bm{x}-\bm{x}^{\prime}|}\bar{1}\
&\frac{-ik}{4\pi}\int
d\bm{x}^{\prime\prime}\frac{e^{ik|\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).$ (44)

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

##  V Self-consistent equation

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

###  V.1 $\bm{x}$-space representation

Let us first formulate the self-consistent equation in terms of
$\bm{\mathcal{A}}(\bm{x};\omega)$ and $\bm{\mathcal{J}}(\bm{x};\omega)$.
Substitution of the constitutive equation into the Maxwell equation yields

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

Since the nonlocal susceptibility
$\bar{\mathcal{X}}(\bm{x}^{\prime},\bm{x}^{\prime\prime};\omega)$ shown in Eq.
(30) is a separable integral kernel, by multiplying Eq. (45) by
$\left(\bm{\mathcal{J}}_{\nu^{\prime}\mu^{\prime}}(\bm{x})\right)^{\rm t}$
from the left and integrating with respect to $\bm{x}$, we obtain

$\left(\hbar\omega_{\nu^{\prime}\mu^{\prime}}-\hbar\omega-i\gamma\right)X_{\nu^{\prime}\mu^{\prime}}^{(-)}(\omega)+\sum_{\mu,\nu}\left[K_{\nu^{\prime}\mu^{\prime},\mu\nu}(\omega)X_{\nu\mu}^{(-)}(\omega)+L_{\nu^{\prime}\mu^{\prime},\nu\mu}(\omega)X_{\mu\nu}^{(+)}(\omega)\right]=Y^{\rm(0)}_{\nu^{\prime}\mu^{\prime}}(\omega).$
(46)

Here, the factors are defined as

$\displaystyle K_{\nu^{\prime}\mu^{\prime},\mu\nu}(\omega)$ $\displaystyle=$
$\displaystyle\mu_{0}\rho_{0,\mu}\int d\bm{x}\int
d\bm{x}^{\prime}\left(\bm{\mathcal{J}}_{\nu^{\prime}\mu^{\prime}}(\bm{x})\right)^{\rm
t}\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)\bm{\mathcal{J}}_{\mu\nu}(\bm{x}^{\prime}),$
(47) $\displaystyle L_{\nu^{\prime}\mu^{\prime},\nu\mu}(\omega)$
$\displaystyle=$ $\displaystyle\mu_{0}\rho_{0,\mu}\int d\bm{x}\int
d\bm{x}^{\prime}\left(\bm{\mathcal{J}}_{\nu^{\prime}\mu^{\prime}}(\bm{x})\right)^{\rm
t}\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)\bm{\mathcal{J}}_{\nu\mu}(\bm{x}^{\prime}),$
(48) $\displaystyle X_{\nu\mu}^{(-)}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar\omega_{\nu\mu}-\hbar\omega-i\gamma}\int
d\bm{x}\left(\bm{\mathcal{J}}_{\nu\mu}(\bm{x})\right)^{\rm
t}\bm{\mathcal{A}}(\bm{x};\omega),$ (49) $\displaystyle
X_{\mu\nu}^{(+)}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar\omega_{\nu\mu}+\hbar\omega+i\gamma}\int
d\bm{x}\left(\bm{\mathcal{J}}_{\mu\nu}(\bm{x})\right)^{\rm
t}\bm{\mathcal{A}}(\bm{x};\omega),$ (50) $\displaystyle
Y^{\rm(0)}_{\nu^{\prime}\mu^{\prime}}(\omega)$ $\displaystyle=$
$\displaystyle\int
d\bm{x}\left(\bm{\mathcal{J}}_{\nu^{\prime}\mu^{\prime}}(\bm{x})\right)^{\rm
t}\bm{\mathcal{A}}_{0}(\bm{x};\omega),$ (51)

where $\rho_{0,\mu}$ is an element of the density matrix in Eq. (34). Note
that $X_{\mu\nu}^{(+)}(\omega)\neq\left(X_{\nu\mu}^{(-)}(-\omega)\right)^{*}$.
Similarly, by multiplying
$\left(\bm{\mathcal{J}}_{\mu^{\prime}\nu^{\prime}}(\bm{x})\right)^{\rm
t}\left[=\left(\bm{\mathcal{J}}_{\nu^{\prime}\mu^{\prime}}^{\
*}(\bm{x})\right)^{\rm t}\right]$, we find

$\left(\hbar\omega_{\nu^{\prime}\mu^{\prime}}+\hbar\omega+i\gamma\right)X_{\mu^{\prime}\nu^{\prime}}^{(+)}(\omega)+\sum_{\mu,\nu}\left[M_{\mu^{\prime}\nu^{\prime},\mu\nu}(\omega)X_{\nu\mu}^{(-)}(\omega)+N_{\mu^{\prime}\nu^{\prime},\nu\mu}(\omega)X_{\mu\nu}^{(+)}(\omega)\right]=Y^{\rm(0)}_{\mu^{\prime}\nu^{\prime}}(\omega),$
(52)

with

$\displaystyle M_{\mu^{\prime}\nu^{\prime},\mu\nu}(\omega)$ $\displaystyle=$
$\displaystyle\mu_{0}\rho_{0,\mu}\int d\bm{x}\int
d\bm{x}^{\prime}\left(\bm{\mathcal{J}}_{\mu^{\prime}\nu^{\prime}}(\bm{x})\right)^{\rm
t}\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)\bm{\mathcal{J}}_{\mu\nu}(\bm{x}^{\prime}),$
(53) $\displaystyle N_{\mu^{\prime}\nu^{\prime},\nu\mu}(\omega)$
$\displaystyle=$ $\displaystyle\mu_{0}\rho_{0,\mu}\int d\bm{x}\int
d\bm{x}^{\prime}\left(\bm{\mathcal{J}}_{\mu^{\prime}\nu^{\prime}}(\bm{x})\right)^{\rm
t}\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)\bm{\mathcal{J}}_{\nu\mu}(\bm{x}^{\prime}).$
(54)

Note again that the excited states $\mu$ indicating an electron–hole pair
include their spin degrees of freedom. The definitions of
$K_{\nu^{\prime}\mu^{\prime},\mu\nu}$, $L_{\nu^{\prime}\mu^{\prime},\nu\mu}$,
$M_{\mu^{\prime}\nu^{\prime},\mu\nu}$, and
$N_{\mu^{\prime}\nu^{\prime},\nu\mu}$ are mathematically related. In a
physical sense, $K$ and $N$ indicate the resonant and anti-resonant radiative
correlations, respectively, and $L$ and $M$ are their couplings. Thus, we
retain this description.

From Eqs. (46) and (52), we construct the matrix form of the self-consistent
equation:

$\left[\bar{\Xi}(\omega)\right]\left(\begin{matrix}\bm{X}^{(-)}(\omega)\\\
\bm{X}^{(+)}(\omega)\\\
\end{matrix}\right)=\left(\begin{matrix}\bm{Y}^{\rm(0,-)}(\omega)\\\
\bm{Y}^{\rm(0.+)}(\omega)\\\ \end{matrix}\right),$ (55)

with

$\bar{\Xi}(\omega)\equiv\left[\left(\begin{matrix}\hbar\bar{\Omega}-\left(\hbar\omega+i\gamma\right)\bar{1}&0\\\
0&\hbar\bar{\Omega}+\left(\hbar\omega+i\gamma\right)\bar{1}\\\
\end{matrix}\right)+\left(\begin{matrix}\bar{K}(\omega)&\bar{L}(\omega)\\\
\bar{M}(\omega)&\bar{N}(\omega)\\\ \end{matrix}\right)\right].$ (56)

Let us introduce an integer $N$, the cut-off number of the electron states,
and define the vectors of $X_{\nu\mu}^{(\pm)}$ as

$\bm{X}^{(\pm)}=\left(\begin{matrix}\bm{X}^{(\pm)}_{0}\\\
\bm{X}^{(\pm)}_{1}\\\ \vdots\\\ \bm{X}^{(\pm)}_{N}\end{matrix}\right),$ (57)

with

$\bm{X}^{(-)}_{\mu}=\left(\begin{matrix}X^{(-)}_{0\mu}\\\ X^{(-)}_{1\mu}\\\
\vdots\\\ X^{(-)}_{N\mu}\end{matrix}\right),$ (58)

and

$\bm{X}^{(+)}_{\mu}=\left(\begin{matrix}X^{(+)}_{\mu 0}\\\ X^{(+)}_{\mu 1}\\\
\vdots\\\ X^{(+)}_{\mu N}\end{matrix}\right).$ (59)

For $Y^{(0)}_{\nu\mu}$,

$\bm{Y}^{(0,\pm)}=\left(\begin{matrix}\bm{Y}^{(0,\pm)}_{0}\\\
\bm{Y}^{(0,\pm)}_{1}\\\ \vdots\\\ \bm{Y}^{(0,\pm)}_{N}\end{matrix}\right),\ \
\ {\rm with}\ \ \ \bm{Y}^{(0,-)}_{\mu}=\left(\begin{matrix}Y^{(0)}_{0\mu}\\\
Y^{(0)}_{1\mu}\\\ \vdots\\\ Y^{(0)}_{N\mu}\end{matrix}\right)\ \ ,\ \
\bm{Y}^{(0,+)}_{\mu}=\left(\begin{matrix}Y^{(0)}_{\mu 0}\\\ Y^{(0)}_{\mu 1}\\\
\vdots\\\ Y^{(0)}_{\mu N}\end{matrix}\right).$ (60)

Therefore, $\bm{Y}^{(0,-)}$ and $\bm{Y}^{(0,+)}$ interchange through
permutation of the elements. For the energy differences $\omega_{\nu\mu}$,

$\bar{\Omega}=\left(\begin{matrix}\bar{\Omega}_{0}&&&\\\
&\bar{\Omega}_{1}&&\\\ &&\ddots&\\\ &&&\bar{\Omega}_{N}\end{matrix}\right)$
(61)

is a diagonal matrix with $\bar{\Omega}_{\mu}={\rm
diag}(\omega_{0\mu},\omega_{1\mu},\omega_{2\mu}\cdots,\omega_{N\mu})$. For the
matrixes $\bar{K}$, $\bar{L}$, $\bar{M}$, and $\bar{N}$,

$\bar{Q}=\left(\begin{matrix}\bar{Q}_{00}&\bar{Q}_{01}&\cdots&\bar{Q}_{0N}\\\
\bar{Q}_{10}&\bar{Q}_{11}&&\\\ \vdots&&\ddots&\\\
\bar{Q}_{N0}&&&\bar{Q}_{NN}\end{matrix}\right)\ \ (Q=K,L,M,N),$ (62)

with the following different subscript rules:

$\displaystyle\bar{K}_{\mu^{\prime}\mu}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}K_{0\mu^{\prime},\mu 0}&K_{0\mu^{\prime},\mu
1}&\cdots&K_{0\mu^{\prime},\mu N}\\\ K_{1\mu^{\prime},\mu
0}&K_{1\mu^{\prime},\mu 1}&&\\\ \vdots&&\ddots&\\\ K_{N\mu^{\prime},\mu
0}&&&K_{N\mu^{\prime},\mu N}\end{matrix}\right),$ (63)
$\displaystyle\bar{L}_{\mu^{\prime}\mu}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}L_{0\mu^{\prime},0\mu}&L_{0\mu^{\prime},1\mu}&\cdots&L_{0\mu^{\prime},N\mu}\\\
L_{1\mu^{\prime},0\mu}&L_{1\mu^{\prime},1\mu}&&\\\ \vdots&&\ddots&\\\
L_{N\mu^{\prime},0\mu}&&&L_{N\mu^{\prime},N\mu}\end{matrix}\right),$ (64)
$\displaystyle\bar{M}_{\mu^{\prime}\mu}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}M_{\mu^{\prime}0,\mu 0}&M_{\mu^{\prime}0,\mu
1}&\cdots&M_{\mu^{\prime}0,\mu N}\\\ M_{\mu^{\prime}1,\mu
0}&M_{\mu^{\prime}1,\mu 1}&&\\\ \vdots&&\ddots&\\\ M_{\mu^{\prime}N,\mu
0}&&&M_{\mu^{\prime}N,\mu N}\end{matrix}\right),$ (65)
$\displaystyle\bar{N}_{\mu^{\prime}\mu}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}N_{\mu^{\prime}0,0\mu}&N_{\mu^{\prime}0,1\mu}&\cdots&N_{\mu^{\prime}0,N\mu}\\\
N_{\mu^{\prime}1,0\mu}&N_{\mu^{\prime}1,1\mu}&&\\\ \vdots&&\ddots&\\\
N_{\mu^{\prime}N,0\mu}&&&N_{\mu^{\prime}N,N\mu}\end{matrix}\right).$ (66)

Taking an inversion matrix of Eq. (56), we obtain the vector of
$X_{\nu\mu}^{(\pm)}$ as

$\left(\begin{matrix}\bm{X}^{(-)}\\\ \bm{X}^{(+)}\\\
\end{matrix}\right)=\left[\bar{\Xi}(\omega)\right]^{-1}\left(\begin{matrix}\bm{Y}^{\rm(0,-)}\\\
\bm{Y}^{\rm(0.+)}\\\ \end{matrix}\right).$ (67)

The elements of this vector are the coefficients of the induced densities, as
Eq. (29) is rewritten as

$\left\langle\hat{\bm{\mathcal{J}}}(\bm{x})\right\rangle(t)=\sum_{\mu,\nu}\rho_{0,\mu}\left[X^{(-)}_{\nu\mu}(\omega)\bm{\mathcal{J}}_{\mu\nu}(\bm{x})+X^{(+)}_{\mu\nu}(\omega)\bm{\mathcal{J}}_{\nu\mu}(\bm{x})\right]e^{-i\omega
t}.$ (68)

Therefore, the excitations induced in the nanostructure are described by the
eigenvalues of matrix $\bar{\Xi}(\omega)$. The electromagnetic fields are also
given in terms of $X_{\nu\mu}^{(\pm)}$, following Eq. (45):

$\bm{\mathcal{A}}(\bm{x};\omega)=\bm{\mathcal{A}}_{0}(\bm{x};\omega)-\mu_{0}\sum_{\mu,\nu}\rho_{0,\mu}\int
d\bm{x}^{\prime}\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)\left[X^{(-)}_{\nu\mu}(\omega)\bm{\mathcal{J}}_{\mu\nu}(\bm{x^{\prime}})+X^{(+)}_{\mu\nu}(\omega)\bm{\mathcal{J}}_{\nu\mu}(\bm{x^{\prime}})\right].$
(69)

At zero temperature, i.e., $T=0$, $\rho_{0,\mu}=\delta_{0,\mu}$ in
$f_{\nu\mu}$ and $h_{\nu\mu}$ defined in Eqs. (32) and (33), respectively.
Then, the energy differences $\hbar\omega_{\nu 0}$ in the denominators of
$X_{\nu 0}^{(-)}$ and $X_{0\nu}^{(+)}$ are positive only. For the matrix form
of the self-consistent equation (55), only $\bm{\mathcal{J}}_{\nu 0}$ and
$\bm{\mathcal{J}}_{0\nu}$ are considered. Therefore, the matrix elements are
reduced for $K_{\nu^{\prime}0,0\nu}$, $L_{\nu^{\prime}0,\nu 0}$,
$M_{0\nu^{\prime},0\nu}$, and $N_{0\nu^{\prime},\nu 0}$. Moreover, we can
approximate the constitutive equation (29), where a non-resonant term by
$\hbar\omega_{\nu 0}+\hbar\omega$ is essentially negligible rather than a
resonant term by $\hbar\omega_{\nu 0}-\hbar\omega$. If we neglect the non-
resonant term in Eq. (55), the self-consistent equation is reduced to

$\left(\hbar\omega_{\nu^{\prime}0}-\hbar\omega-i\gamma\right)X_{\nu^{\prime}0}^{(-)}(\omega)+\sum_{\nu}K_{\nu^{\prime}0,0\nu}(\omega)X_{\nu
0}^{(-)}(\omega)=Y^{\rm(0)}_{\nu^{\prime}0}(\omega).$ (70)

The matrix form of the self-consistent equation (55) has an important
advantage. The matrix $\bar{\Xi}(\omega)$ with $\bar{K}$, $\bar{L}$,
$\bar{M}$, and $\bar{N}$ consisting of $\bar{\mathcal{G}}$ and
$\left\langle\nu\left|\hat{\bm{\mathcal{J}}}\right|\mu\right\rangle$ gives the
eigenmodes of the system coupled with the radiative field. The L and T fields
are in $\bar{\mathcal{G}}$ and the properties of the nonlocal response are
given by $|\mu\rangle$. From the real and imaginary components of the
eigenvalues of $\bar{\Xi}(\omega)$, the properties of the collective and
individual excitations are evaluated. Then, the current and charge densities
are discussed in terms of the plasmons and SPEs.

###  V.2 $\bm{k}$-space representation

For a metallic structure having some symmetry, the Fourier transformation of
the fields and densities is useful for the evaluation. In that case, we
consider the $\bm{k}$-space representation for the self-consistent equation.
The fields and densities in the $\bm{k}$-representation are given as

$\tilde{\bm{\mathcal{O}}}(\bm{k};\omega)=\frac{1}{(2\pi)^{\frac{3}{2}}}\int
d\bm{x}\bm{\mathcal{O}}(\bm{x};\omega)e^{-i\bm{k}\cdot\bm{x}},\hskip
14.22636pt(\bm{\mathcal{O}}=\bm{\mathcal{A}},\bm{\mathcal{J}}).$ (71)

The Maxwell equations (40) and their solutions (42) are also transferred to
the $\bm{k}$-representation, with the latter being expressed as

$\tilde{\bm{\mathcal{A}}}(\bm{k};\omega)=\tilde{\bm{\mathcal{A}}}_{0}(\bm{k};\omega)-\mu_{0}\bar{\mathcal{G}}_{\bm{k}}(\omega)\tilde{\bm{\mathcal{J}}}(\bm{k};\omega).$
(72)

Here, the Green’s function is given as

$\bar{\mathcal{G}}_{\bm{k}}(\omega)=\left(\begin{matrix}\frac{1}{-\bm{k}^{2}+\omega^{2}/c^{2}}\bar{1}&-\frac{1}{\bm{k}^{2}}\frac{1}{-\bm{k}^{2}+\omega^{2}/c^{2}}\frac{\omega}{c}\bm{k}\\\
0&-\frac{1}{\bm{k}^{2}}\end{matrix}\right).$ (73)

The incident term satisfies

$\left(\begin{matrix}\left(-\bm{k}^{2}+\frac{\omega^{2}}{c^{2}}\right)\bar{1}&-\frac{\omega}{c}\bm{k}\\\
0&-\bm{k}^{2}\end{matrix}\right)\tilde{\bm{\mathcal{A}}}_{0}(\bm{k};\omega)=0.$
(74)

Note that the solution of the Maxwell equations contains no integral.

The constitutive equation becomes

$\tilde{\bm{\mathcal{J}}}(\bm{k};\omega)=\int
d\bm{k}^{\prime}\bar{\mathcal{X}}_{\bm{k},\bm{k}^{\prime}}(\omega)\tilde{\bm{\mathcal{A}}}(\bm{k}^{\prime};\omega),$
(75)

with

$\bar{\mathcal{X}}_{\bm{k},\bm{k}^{\prime}}(\omega)=\sum_{\mu,\nu}\left[f_{\nu\mu}\tilde{\bm{\mathcal{J}}}_{\mu\nu}(\bm{k})\left(\tilde{\bm{\mathcal{J}}}_{\nu\mu}(-\bm{k}^{\prime})\right)^{\rm
t}+h_{\nu\mu}\tilde{\bm{\mathcal{J}}}_{\nu\mu}(\bm{k})\left(\tilde{\bm{\mathcal{J}}}_{\mu\nu}(-\bm{k}^{\prime})\right)^{\rm
t}\right].$ (76)

Equations (72) and (75) form the self-consistent equation. The matrix form of
the self-consistent equation in the $\bm{k}$-representation is obtained using
the approach described in the previous subsection, with

$\displaystyle\left(\hbar\omega_{\nu^{\prime}\mu^{\prime}}-\hbar\omega-i\gamma\right)X_{\nu^{\prime}\mu^{\prime}}^{(-)}(\omega)+\sum_{\mu,\nu}\left[\tilde{K}_{\nu^{\prime}\mu^{\prime},\mu\nu}(\omega)\tilde{X}_{\nu\mu}^{(-)}(\omega)+\tilde{L}_{\nu^{\prime}\mu^{\prime},\nu\mu}(\omega)\tilde{X}_{\mu\nu}^{(+)}(\omega)\right]$
$\displaystyle=$
$\displaystyle\tilde{Y}^{\rm(0)}_{\nu^{\prime}\mu^{\prime}}(\omega),$ (77)
$\displaystyle\left(\hbar\omega_{\nu^{\prime}\mu^{\prime}}+\hbar\omega+i\gamma\right)X_{\mu^{\prime}\nu^{\prime}}^{(+)}(\omega)+\sum_{\mu,\nu}\left[\tilde{M}_{\mu^{\prime}\nu^{\prime},\mu\nu}(\omega)\tilde{X}_{\nu\mu}^{(-)}(\omega)+\tilde{N}_{\mu^{\prime}\nu^{\prime},\nu\mu}(\omega)\tilde{X}_{\mu\nu}^{(+)}(\omega)\right]$
$\displaystyle=$
$\displaystyle\tilde{Y}^{\rm(0)}_{\mu^{\prime}\nu^{\prime}}(\omega),$ (78)

where the matrix elements are defined, for example, as

$\displaystyle\tilde{K}_{\nu^{\prime}\mu^{\prime},\mu\nu}(\omega)$
$\displaystyle=$ $\displaystyle\mu_{0}\rho_{0,\mu}\int
d\bm{k}\left(\tilde{\bm{\mathcal{J}}}_{\nu^{\prime}\mu^{\prime}}(-\bm{k})\right)^{\rm
t}\bar{\mathcal{G}}_{\bm{k}}(\omega)\tilde{\bm{\mathcal{J}}}_{\mu\nu}(\bm{k}),$
(79) $\displaystyle\tilde{X}_{\nu\mu}^{(-)}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar\omega_{\nu\mu}-\hbar\omega-i\gamma}\int
d\bm{k}\left(\tilde{\bm{\mathcal{J}}}_{\nu\mu}(-\bm{k})\right)^{\rm
t}\tilde{\bm{\mathcal{A}}}(\bm{k};\omega).$ (80)

The Fourier transformation demonstrates that these elements are equivalent to
those of the $\bm{x}$-representation,
$\tilde{K}_{\nu^{\prime}\mu^{\prime},\mu\nu}(\omega)=K_{\nu^{\prime}\mu^{\prime},\mu\nu}(\omega)$,
$\tilde{X}_{\nu\mu}^{(-)}(\omega)=X_{\nu\mu}^{(-)}(\omega)$, etc. Therefore,
the matrix element evaluation can be performed in either the $\bm{x}$\- and
$\bm{k}$-representations depending on the case.

It is worth noting that, in
$\left(\tilde{\bm{\mathcal{J}}}_{\nu^{\prime}\mu^{\prime}}(-\bm{k})\right)^{\rm
t}\bar{\mathcal{G}}_{\bm{k}}(\omega)\tilde{\bm{\mathcal{J}}}_{\mu\nu}(\bm{k})$,
a component between the induced charge densities,

$\tilde{\rho}_{\nu^{\prime}\mu^{\prime}}(-\bm{k})\left(-\frac{1}{\bm{k}^{2}}\right)\tilde{\rho}_{\mu\nu}(\bm{k}),$
(81)

describes the Coulomb interaction, which arises when an external field is
applied. If we apply a free electron model to the basis $|\mu\rangle$ of
$\hat{H}_{0}$, the collective excitation (plasmon) in the RPA can be obtained.
In addition to this
$\tilde{\rho}_{\nu^{\prime}\mu^{\prime}}\left[\bar{\mathcal{G}}_{\bm{k}}\right]_{\rho\rho}\tilde{\rho}_{\mu\nu}$
term, $\tilde{K}_{\nu^{\prime}\mu^{\prime},\mu\nu}$ and the other matrix
components have
$\tilde{\bm{j}}_{\nu^{\prime}\mu^{\prime}}\left[\bar{\mathcal{G}}_{\bm{k}}\right]_{j\rho}\tilde{\rho}_{\mu\nu}$
and
$\tilde{\bm{j}}_{\nu^{\prime}\mu^{\prime}}\left[\bar{\mathcal{G}}_{\bm{k}}\right]_{jj}\tilde{\bm{j}}_{\mu\nu}$
terms, which have been neglected from many previous studies. However, the
modifications of the eigenstates and spectrum of the electronic system due to
the coupled radiative field are given by these components, as noted in the
previous subsection. The nanostructure shape determines the spectrum shift via
the $|\mu\rangle$ in $\tilde{\rho}_{\mu\nu}$ and $\tilde{\bm{j}}_{\mu\nu}$.
Therefore, analysis and exploration of an enhancement condition based on
$\tilde{K}$, $\tilde{L}$, $\tilde{M}$, and $\tilde{N}$, with
$\tilde{\rho}_{\nu^{\prime}\mu^{\prime}}\left[\bar{\mathcal{G}}_{\bm{k}}\right]_{\rho\rho}\tilde{\rho}_{\mu\nu}$,
$\tilde{\bm{j}}_{\nu^{\prime}\mu^{\prime}}\left[\bar{\mathcal{G}}_{\bm{k}}\right]_{j\rho}\tilde{\rho}_{\mu\nu}$,
and
$\tilde{\bm{j}}_{\nu^{\prime}\mu^{\prime}}\left[\bar{\mathcal{G}}_{\bm{k}}\right]_{jj}\tilde{\bm{j}}_{\mu\nu}$
are important.

##  VI Discussion

In the present formulation of (45) and its matrix form (55), the optical
response of the nanostructures is described in terms of
$\bm{\mathcal{J}}_{\nu\mu}(\bm{x})$. The electron eigenstates $|\mu\rangle$
contribute to the optical response via $\bm{\mathcal{J}}_{\nu\mu}(\bm{x})$ in
the susceptibility. When one prepares the $|\mu\rangle$ of $\hat{H}_{0}$,
including the Coulomb interaction at equilibrium, the nanostructure optical
responses, e.g., the plasmon–polariton and SPE spectra, are obtained within a
theoretical framework for the prepared eigenstates. Recently, the first-
principles calculation for a nanostructure or a nanoscale cluster of atoms was
developed Kumar19 ; Zhang17 ; Chapkin18 ; Senanayake19 ; Ma15 ; Schira19 ;
Lu19 , where the electronic states with the electron–electron interaction (the
L component of the EM field) are evaluated precisely. In those previous
studies, however, the T field was not self-consistently considered. In our
formulation, both components of the EM fields are determined self-consistently
on the mesoscopic scale. Further, if one employs the equilibrium states
evaluated in the first-principles calculation in our formulation, the T fields
generated by the electronic responses and the electron—electron interactions
are taken into account com1 . Hence, our approach provides an important and
convenient framework for investigations of mesoscopic plasmonics and
photonics.

The present formulation is based on the classical Maxwell equations for the
incident and induced EM fields. As the induced L field describes part of the
electron–electron interaction, the quantum Maxwell equations for the EM fields
are required for the higher-order electron correlations. For the linear
response framework discussed in this paper, however, the classical treatment
of the fields is equivalent to the quantum treatment.

Let us emphasize the difference in the four-vector representation of Eq. (55)
compared to the conventional nonlocal and self-consistent theory book:cho1 .
In our formulation, both the vector and scalar potentials are fully
considered. The Coulomb gauge separates the L and T components in
$\bm{\mathcal{A}}=(\bm{A},\phi/c)$, which directly provides a picture of the L
and T component mixing due to the nanostructure in the self-consistent
relation in Eq. (55). For the densities $\bm{\mathcal{J}}=(\bm{j},c\rho)$ in
the constitutive equation (29), the current is induced not only by the T
field, but also by the L field, even though $\hat{H}_{\rm
ind}\sim\hat{\bm{j}}\cdot\bm{A}+\hat{\rho}\phi$. The charge density is also
attributed to both the L and T fields. These characteristics are due to the
off-diagonal elements in the $4\times 4$ susceptibility $\bar{\mathcal{X}}$.
By $\bar{\mathcal{G}}$, the L field $\phi$ is generated by the charge only,
whereas both the current and charge generate the T field $\bm{A}$. This “cross
generation” of LT components is essential physics for the nanostructure and is
enlarged when the nanostructure enhances its $\bm{\chi}_{j\rho}$ and
$\bm{\chi}_{\rho j}$ due to the nonlocality. Therefore, the present
formulation provides a guideline for plasmon resonance enhancement since the
current and charge density generation is related to the collective and
individual excitations.

In our formulation, the optical responses of electrons are described in terms
of the current and charge densities. This representation is useful for
considering the responses to the L and T field components. In many previous
studies, the optical properties were discussed in terms of the plasmons and
SPEs. In nanostructures featuring nonlocality, the plasmons and SPEs are
coupled with each other. Hence, the relation between the densities
($\bm{j},c\rho$) and the collective and individual excitations is complicated.
Here, the excitation energy spectrum is evaluated from the eigenvalues of
$\bar{\Xi}(\omega)$ in Eq. (56) as the excitations correspond to the poles of
the system with the radiative field. Hence, the densities $\bm{\mathcal{J}}$
and induced fields $\bm{\mathcal{A}}$ are described in terms of the
excitations in Eqs. (68) and (69) with Eq. (67). This is an advantage over
another existing treatment of the nonlocal effect Pitarke07 . Moreover, from
the poles, coherent coupling between the plasmon (collective excitation) and
the carrier generation due to the SPE can be investigated.

The coherent coupling of the plasmon-like and single-particle-like excitations
in the nanostructures gives rise to a bidirectional energy transfer Ma15 . The
plasmon and individual excitations form discrete and continuous spectra,
respectively. Therefore, their coupling may demonstrate the Fano resonance as
evidence of coherent coupling. In such a scenario, the self-consistency of the
L and T total fields and the induced charge and current density are important.
In addition to the energy, the electron wavefunction in the nanostructures may
significantly affect the coupling. Therefore, our microscopic-theory-based
formulation reveals such light–nanostructure interaction and its enhancement.
As a multi-dimensional integral is included in Eqs. (47)–(51), the presence of
numerous electrons, even in submicro-scale metallic structures, would make a
feasibility of numerical calculations difficult. The approach based on the
first-principles calculation Ma15 ; You18  has more limited applications than
the analytical approach based on our formulation. Note also that we can
overcome this difficulty for several nanostructures, e.g., a 2D sheet
(including graphene) and rectangular rod. We will demonstrate these
applications in a subsequent paper Yokoyama21 , in which an analytical
evaluation is performed. In that evaluation, the analytical expression reveals
a clear relation between the collective excitation and induced densities. Such
developments, which are achievable using our microscopic approach, are not
limited to several adaptable nanostructures, but can be applied to arbitrary
nanostructures by considering appropriate approximations to restrict the
bases. Our approach can also be applied to an array of two or more
nanostructure units, in which the nonlocal effect is enlarged with the
nanogap. Hence, our theory is expected to reveal an energy transfer through
cooperation between the collective excitation and radiative field via an
appropriately fabricated plasmonic array.

##  VII Summary

To correctly describe coherent coupling between collective excitations
(plasmons) and single-particle excitations (SPEs) via a transverse (T)
electric field, which is an aspect that has been neglected from previous
studies, we developed a self-consistent and microscopic nonlocal formulation
for the light–matter interaction in plasmonic nanostructures. Our formulation
is based on linear response theory with a nonlocal susceptibility and the
classical Maxwell equations describing the T and longitudinal (L)
electromagnetic fields. The nonlocal susceptibility
$\bar{\mathcal{X}}(\bm{x},\bm{x}^{\prime})$ is obtained from the interaction
Hamiltonian and the eigenstates of the nanostructure. The Coulomb interaction
is included not only in the non-perturbative Hamiltonian $\hat{H}_{0}$, but
also in the L electric fields induced by the collectively excited “plasmons”
and the carriers generated by the “SPEs.” The developed formulation can be
applied to various frameworks by utilizing the electronic states obtained
according to those respective frameworks. In the formulation, the L and T
components of the fields and optically induced responses are described in
terms of the four-vector potential $\bm{\mathcal{A}}(\bm{x},t)$ and density
$\bm{\mathcal{J}}(\bm{x},t)$ with the Green’s function
$\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime})$ in the four-vector space. The
self-consistent equation is rewritten in matrix form with the incident and
induced fields, $Y^{(0,\pm)}_{\mu\nu}$ and $X^{(\pm)}_{\mu\nu}$. From the
poles of this formulated matrix, the excitation spectrum can be examined.

Our formulation could be combined with real-time density functional theory
based on the first-principles calculation, in which the electron–electron
interaction is considered accurately. The four-vector formulation can describe
the effects of T fields on the eigenstates given by the first-principles
calculation at each step in a real-time simulation. In this treatment, the
exchange correlation and quantum fluctuation can be partially considered via
the excited states. Application of the present formulation to various types of
metallic meso- and nanostructures will facilitate theoretical understanding of
the interplay between collective and individual excitations in nanoscale
systems, which will aid the design of metallic systems to control photo-
induced hot carrier generation.

###### Acknowledgements.

T.Y. and H.I. are supported by JSPS KAKENHI JP16H06504 in Scientific Research
on Innovative Areas, “Nano-Material Optical-Manipulation.” T.Y. is supported
by JSPS KAKENHI 18K13484. H.I. is supported by JSPS KAKENHI 18H01151.

##  Appendix A L-T separation of Maxwell’s equations

In the main text, we formulated the self-consistent equation in terms of
$\\{(\bm{A},\phi/c),(\bm{j},c\rho)\\}$ with the Coulomb gauge. Let us consider
another description in terms of
$\\{(\bm{A}^{\rm(T)},\phi^{\rm(L)}/c),(\bm{j}^{\rm(T)},\bm{j}^{\rm(L)})\\}$.
Here, to emphasize the separation of the L and T components in the fields and
densities, we present the superscripts explicitly. Equation (37) can be
expressed as

$\displaystyle\left(\bm{\nabla}^{2}-\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\right)\bm{A}^{\rm(T)}(\bm{x},t)$ $\displaystyle=$
$\displaystyle-\mu_{0}\bm{j}^{\rm(T)}(\bm{x},t),$ (82)
$\displaystyle-\frac{1}{c^{2}}\frac{\partial}{\partial
t}\bm{\nabla}\phi^{\rm(L)}(\bm{x},t)$ $\displaystyle=$
$\displaystyle-\mu_{0}\bm{j}^{\rm(L)}(\bm{x},t),$ (83)

where $\bm{j}(\bm{x},t)=\bm{j}^{\rm(T)}(\bm{x},t)+\bm{j}^{\rm(L)}(\bm{x},t)$.
The scalar potential $\phi^{\rm(L)}$ generates the L component only. Equation
(38) also describes the L only, where

$\bm{\nabla}^{2}\phi^{\rm(L)}(\bm{x},t)=-\frac{\rho^{\rm(L)}(\bm{x},t)}{\varepsilon_{0}}.$
(84)

Equations (83) and (84) satisfy the continuous relation for the L component,

$\bm{\nabla}\cdot\bm{j}^{\rm(L)}(\bm{x},t)+\frac{\partial\rho^{\rm(L)}(\bm{x},t)}{\partial
t}=0.$ (85)

Let us refer to the T component in the scalar potential as $\phi^{\rm(T)}$.
The T component should be
$\bm{\nabla}\cdot\bm{E}=\bm{\nabla}\cdot\left(-\bm{\nabla}\phi^{\rm(T)}\right)=0$.
Thus, it must have a linear dependence up to $\phi^{\rm(T)}(\bm{x},t)\sim
ax+b$. These terms cannot be induced from the internal charge density of the
system. We incorporate the external scalar potential $\phi_{\rm
ext}(\bm{x},t)$ in the interaction Hamiltonian (25). The T component in the
external field, however, is included in the vector potential. Therefore, the
scalar potential induces the L component only, and we find
$\bm{E}^{\rm(L)}=-\bm{\nabla}\phi^{\rm(L)}$ and
$\bm{E}^{\rm(T)}=-\partial_{t}\bm{A}^{\rm(T)}$ for
$\bm{E}=\bm{E}^{\rm(L)}+\bm{E}^{\rm(T)}$.

Such L and T field separation has been also discussed by Cho book:cho2 , where
the scalar potential is reduced from the Hamiltonian (or Lagrangian) and only
the T components, $\bm{A}^{\rm(T)}$ and $\bm{j}^{\rm(T)}$, is described in the
Maxwell equations. The L field, $\bm{E}^{\rm(L)}$, is treated as the external
field only.

The representation of ($\bm{j}^{\rm(T)},\bm{j}^{\rm(L)}$) is a modification of
($\bm{j},c\rho$). It may be useful to separate the L and T components clearly
in the source terms. Moreover, the Green’s function for Eqs. (82) and (83) is
simpler than Eq. (44), as the off-diagonal components are absent. However, for
this representation, the constitutive equation (29) [and the perturbative
Hamiltonian (25)] should be rewritten in terms of
($\bm{j}^{\rm(T)},\bm{j}^{\rm(L)}$). Then, the susceptibility
$\bar{\mathcal{X}}$ must be a $6\times 4$ matrix, which is not desirable in
terms of the mathematics.

##  Appendix B Derivation of Green’s function

The Green’s function for the Maxwell equations is defined according to the
differential operators in the equation, which corresponds to the matrix
$\bar{\mathcal{D}}(\bm{x};\omega)$ in Eq. (40). Therefore, the Green’s
function $\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)$ is the $4\times 4$
matrix,

$\bar{\mathcal{D}}(\bm{x};\omega)\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)=\bar{1}\delta(\bm{x}-\bm{x}^{\prime}).$
(86)

The formal solution of the Maxwell equation is

$\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).$
(87)

The first term in Eq. (87) satisfies

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

To obtain an explicit form, $\bar{\mathcal{G}}$ is divided into the following
matrix elements:

$\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)=\left(\begin{matrix}\bar{g}_{AA}(\bm{x},\bm{x}^{\prime};\omega)&\bm{g}_{A\phi}(\bm{x},\bm{x}^{\prime};\omega)\\\
\bm{g}_{\phi A}^{\ {\rm
t}}(\bm{x},\bm{x}^{\prime};\omega)&g_{\phi\phi}(\bm{x},\bm{x}^{\prime};\omega)\end{matrix}\right),$
(89)

where $\bar{g}_{AA}$ is a $3\times 3$ matrix, $\bar{g}_{\phi\phi}$ is scalar,
and the other elements are vectors.

We have a mathematical preparation. The solution of

$\left(\bm{\nabla}^{2}+k^{2}\right)g_{k}(\bm{x},\bm{x}^{\prime})=\delta(\bm{x}-\bm{x}^{\prime})$
(90)

is given as

$g_{k}(\bm{x},\bm{x}^{\prime})=-\frac{e^{ik|\bm{x}-\bm{x}^{\prime}|}}{4\pi|\bm{x}-\bm{x}^{\prime}|}.$
(91)

Here, we avoid constant and linear terms,
$\alpha_{0}+\bm{\alpha}_{1}\cdot\bm{x}$, in Eq. (91), because of the boundary
condition $g(|\bm{x}|\to\infty)=0$.

The scalar potential component of the Green’s function satisfies

$\bm{\nabla}^{2}g_{\phi\phi}(\bm{x},\bm{x}^{\prime};\omega)=\delta(\bm{x}-\bm{x}^{\prime}),$
(92)

the solution of which is given by Eq. (91) as follows:

$g_{\phi\phi}(\bm{x},\bm{x}^{\prime};\omega)=g_{0}(\bm{x},\bm{x}^{\prime})=-\frac{1}{4\pi|\bm{x}-\bm{x}^{\prime}|}.$
(93)

Further, $g_{\phi\phi}$ gives another matrix element $\bm{g}_{A\phi}$. From
Eq. (86),

$\left(\bm{\nabla}^{2}+k^{2}\right)\bm{g}_{A\phi}(\bm{x},\bm{x}^{\prime};\omega)+ik\bm{\nabla}g_{\phi\phi}(\bm{x},\bm{x}^{\prime};\omega)=0.$
(94)

Here, $k=\omega/c$. By substituting Eq. (93),

$\left(\bm{\nabla}^{2}+k^{2}\right)\bm{g}_{A\phi}(\bm{x},\bm{x}^{\prime};\omega)=-ik\bm{\nabla}g_{\phi\phi}(\bm{x},\bm{x}^{\prime};\omega)=-\frac{ik}{4\pi}\frac{\bm{x}-\bm{x}^{\prime}}{|\bm{x}-\bm{x}^{\prime}|^{3}}.$
(95)

This is an inhomogeneous equation, which is related to Eq. (90). The solution
of $\bm{g}_{A\phi}$ is also obtained using Eq. (91), with

$\bm{g}_{A\phi}(\bm{x},\bm{x}^{\prime};\omega)=\bm{g}_{k}(\bm{x})+\int
d\bm{x}^{\prime\prime}g_{k}(\bm{x},\bm{x}^{\prime\prime})\left(-\frac{ik}{4\pi}\frac{\bm{x}^{\prime\prime}-\bm{x}^{\prime}}{|\bm{x}^{\prime\prime}-\bm{x}^{\prime}|^{3}}\right).$
(96)

The first term satisfies
$\left(\bm{\nabla}^{2}+k^{2}\right)\bm{g}_{k}(\bm{x})=0$. It gives
$\bm{g}_{k}(\bm{x})=A_{0}e^{i\bm{k}\cdot(\bm{x}-\bm{x}_{0})}$ with
$|\bm{k}|=\omega/c$. Here, $A_{0}$ and $\bm{x}_{0}$ are constant. However,
this term should be zero because of the boundary condition at
$|\bm{x}|\to\infty$. Then, we obtain

$\bm{g}_{A\phi}(\bm{x},\bm{x}^{\prime};\omega)=\frac{ik}{(4\pi)^{2}}\int
d\bm{x}^{\prime\prime}\frac{e^{ik|\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}}.$
(97)

Next, we consider $\bar{g}_{AA}$ and $\bm{g}_{\phi A}^{\ {\rm t}}$. From Eq.
(86), $\bm{g}_{\phi A}^{\ {\rm t}}$ satisfies

$\bm{\nabla}^{2}\bm{g}_{\phi A}^{\ {\rm t}}(\bm{x},\bm{x}^{\prime};\omega)=0.$
(98)

Note that $\bm{g}_{\phi A}^{\ {\rm t}}$ is a $1\times 3$ row vector. The
solution is $\bm{g}_{\phi A}^{\ {\rm
t}}(\bm{x},\bm{x}^{\prime};\omega)=\left\\{\bar{\beta}_{1}\bm{x}+\bm{\beta}_{0}\right\\}^{\rm
t}$, where $\bar{\beta}_{1}$ and $\bm{\beta}_{0}$ are a constant matrix and
vector, respectively. For the boundary condition, $\bm{\beta}_{0}=0$ and
$\bar{\beta}_{1}=0$, and

$\bm{g}_{\phi A}^{\ {\rm t}}(\bm{x},\bm{x}^{\prime};\omega)=0.$ (99)

The last element of the Green’s function, $\bar{g}_{AA}$, satisfies

$\left(\bm{\nabla}^{2}+k^{2}\right)\bar{g}_{AA}(\bm{x},\bm{x}^{\prime};\omega)=\bar{1}\delta(\bm{x}-\bm{x}^{\prime}),$
(100)

which is equivalent to Eq. (90), and the solution is

$\bar{g}_{AA}(\bm{x},\bm{x}^{\prime};\omega)=\bar{1}g_{k}(\bm{x},\bm{x}^{\prime})=-\frac{e^{ik|\bm{x}-\bm{x}^{\prime}|}}{4\pi|\bm{x}-\bm{x}^{\prime}|}\bar{1}.$
(101)

Finally, Eqs. (93), (99), (97), and (101) give the Green’s function in matrix
form:

$\bar{\mathcal{G}}(\bm{x},\bm{x}^{\prime};\omega)=-\frac{1}{4\pi}\left(\begin{matrix}\frac{\
e^{ik|\bm{x}-\bm{x}^{\prime}|}}{|\bm{x}-\bm{x}^{\prime}|}\bar{1}\
&\frac{-ik}{4\pi}\int
d\bm{x}^{\prime\prime}\frac{e^{ik|\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).$ (102)

This Green’s function provides a picture of the field (potential) generation
due to the electronic excitations (current and charge densities).

The Green’s function in the $\bm{k}$-space representation (73) is obtained via
a conventional Fourier transformation of Eq. (40). By applying

$\bm{\mathcal{O}}(\bm{x};\omega)=\frac{1}{(2\pi)^{\frac{3}{2}}}\int
d\bm{k}\tilde{\bm{\mathcal{O}}}(\bm{k};\omega)e^{i\bm{k}\cdot\bm{x}}$ (103)

for $\bm{\mathcal{O}}=\bm{\mathcal{A}},\bm{\mathcal{J}}$, we find

$\displaystyle\bar{\mathcal{D}}(\bm{x};\omega)\left\\{\frac{1}{(2\pi)^{\frac{3}{2}}}\int
d\bm{k}\tilde{\bm{\mathcal{A}}}(\bm{k};\omega)e^{i\bm{k}\cdot\bm{x}}\right\\}$
$\displaystyle=$ $\displaystyle\frac{1}{(2\pi)^{\frac{3}{2}}}\int
d\bm{k}\left(\begin{matrix}\left(-\bm{k}^{2}+\frac{\omega^{2}}{c^{2}}\right)\bar{1}&-\frac{\omega}{c}\bm{k}\\\
0&-\bm{k}^{2}\end{matrix}\right)\tilde{\bm{\mathcal{A}}}(\bm{k};\omega)e^{i\bm{k}\cdot\bm{x}}$
(104) $\displaystyle=$ $\displaystyle\frac{1}{(2\pi)^{\frac{3}{2}}}\int
d\bm{k}(-\mu_{0})\tilde{\bm{\mathcal{J}}}(\bm{k};\omega)e^{i\bm{k}\cdot\bm{x}}.$

Hence, the Green’s function is

$\bar{\mathcal{G}}_{\bm{k}}(\omega)=\left(\begin{matrix}\left(-\bm{k}^{2}+\frac{\omega^{2}}{c^{2}}\right)\bar{1}&-\frac{\omega}{c}\bm{k}\\\
0&-\bm{k}^{2}\end{matrix}\right)^{-1}=\left(\begin{matrix}\frac{1}{-\bm{k}^{2}+\omega^{2}/c^{2}}\bar{1}&-\frac{1}{\bm{k}^{2}}\frac{1}{-\bm{k}^{2}+\omega^{2}/c^{2}}\frac{\omega}{c}\bm{k}\\\
0&-\frac{1}{\bm{k}^{2}}\end{matrix}\right).$ (105)

The formal solutions of $\bm{\mathcal{A}}=(\bm{A},\phi/c)$ can also be
obtained from the separated Maxwell equations in Eqs. (37) and (38). The
latter equation gives

$\phi(\bm{x};\omega)=\phi_{0}(\bm{x};\omega)-\frac{1}{\varepsilon_{0}}\int
d\bm{x}^{\prime}g_{\phi}(\bm{x},\bm{x}^{\prime})\rho(\bm{x}^{\prime};\omega).$
(106)

The scalar Green’s function $g_{\phi\phi}(\bm{x},\bm{x}^{\prime})$ is
equivalent to $g_{\phi\phi}(\bm{x},\bm{x}^{\prime})$ in Eq. (93). The first
term satisfies $\bm{\nabla}^{2}\phi_{0}(\bm{x};\omega)=0$. In the Coulomb
gauge, however, the T field is described by $\bm{A}(\bm{x};\omega)$ only and,
hence, $\phi_{0}=0$. The solution of the former equation is

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

Here, the dyadic Green’s function $\bar{G}_{\bm{A}}$ for the vector potential
is expressed as

$\displaystyle\bar{G}_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)$
$\displaystyle=$ $\displaystyle
g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)\bar{1}-\int
d\bm{x}^{\prime\prime}g_{\bm{A}}(\bm{x},\bm{x}^{\prime\prime};\omega)\bm{\nabla}_{\bm{x}^{\prime\prime}}g_{\phi}(\bm{x}^{\prime\prime},\bm{x}^{\prime})\bm{\nabla}_{\bm{x}^{\prime}}\cdot$
(108) $\displaystyle=$ $\displaystyle
g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)\bar{1}-\int
d\bm{x}^{\prime\prime}\left[\bm{\nabla}_{\bm{x}^{\prime\prime}}g_{\bm{A}}(\bm{x},\bm{x}^{\prime\prime};\omega)\right]\left[\bm{\nabla}_{\bm{x}^{\prime}}g_{\phi}(\bm{x}^{\prime\prime},\bm{x}^{\prime})\right]\cdot$

with

$g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)=-\frac{e^{i\frac{\omega}{c}|\bm{x}-\bm{x}^{\prime}|}}{4\pi|\bm{x}-\bm{x}^{\prime}|}.$
(109)

For the derivation, Eq. (106) is substituted and the continuous equation
$\bm{\nabla}\cdot\bm{j}(\bm{x};\omega)-i\omega\rho(\bm{x};\omega)=0$ is used,
as the source term is described by $\bm{j}(\bm{x};\omega)$ only in this
representation. In the dyadic Green’s function
$\bar{G}_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)$ in Eq. (108), $g_{\phi}$ is
included, which implies that the L component is related to the vector
potential. This seems strange at first glance, as the vector potential in Eq.
(107) must satisfy $\bm{\nabla}_{x}\cdot\bm{A}=0$.

To see it, we separate the T and L components of the current density,
$\bm{j}=\bm{j}^{\rm(T)}+\bm{j}^{\rm(L)}$. This is a description in terms of
$\\{(\bm{A}^{\rm(T)},\phi^{\rm(L)}/c),(\bm{j}^{\rm(T)},\bm{j}^{\rm(L)})\\}$,
which is discussed in Appendix A. The formal solution for Eq. (82) is given by
$g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)$,

$\bm{A}^{\rm(T)}(\bm{x};\omega)=\bm{A}_{0}^{\rm(T)}(\bm{x};\omega)-\mu_{0}\int
d\bm{x}^{\prime}g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)\bm{j}^{\rm(T)}(\bm{x}^{\prime};\omega).$
(110)

Note that Eq. (110) differs slightly from Eq. (107) in its second term on the
r.h.s. For the scalar potential, we apply the continuous relation,
$i\omega\rho=\bm{\nabla}\cdot\bm{j}^{\rm(L)}$, to Eq. (106),

$\displaystyle\phi^{\rm(L)}(\bm{x};\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{i\omega\varepsilon_{0}}\int
d\bm{x}^{\prime}g_{\phi}(\bm{x},\bm{x}^{\prime})\bm{\nabla}_{x^{\prime}}\cdot\bm{j}^{\rm(L)}(\bm{x}^{\prime};\omega).$
(111)

The set of Eqs. (110) and (111) give the solutions of the Maxwell’s equation
for the L–T separation. For the source term in Eq. (110), we apply Eqs. (83),
(111), and the continuous relation,

$\displaystyle\bm{j}^{\rm(T)}(\bm{x}^{\prime};\omega)$ $\displaystyle=$
$\displaystyle\bm{j}(\bm{x}^{\prime};\omega)-\bm{j}^{\rm(L)}(\bm{x}^{\prime};\omega)$
(112) $\displaystyle=$
$\displaystyle\bm{j}(\bm{x}^{\prime};\omega)-\frac{-i\omega}{\mu_{0}c^{2}}\bm{\nabla}_{x^{\prime}}\phi^{\rm(L)}(\bm{x}^{\prime};\omega)$
$\displaystyle=$
$\displaystyle\bm{j}(\bm{x}^{\prime};\omega)-\bm{\nabla}_{x^{\prime}}\int
d\bm{x}^{\prime\prime}g_{\phi}(\bm{x}^{\prime},\bm{x}^{\prime\prime})\bm{\nabla}_{x^{\prime\prime}}\cdot\bm{j}^{\rm(L)}(\bm{x}^{\prime\prime};\omega)$
$\displaystyle=$ $\displaystyle\int
d\bm{x}^{\prime\prime}\left\\{\bar{1}\delta(\bm{x}^{\prime}-\bm{x}^{\prime\prime})-\bm{\nabla}_{x^{\prime}}g_{\phi}(\bm{x}^{\prime},\bm{x}^{\prime\prime})\bm{\nabla}_{x^{\prime\prime}}\cdot\right\\}\bm{j}(\bm{x}^{\prime\prime};\omega).$

Then, Eqs. (110) and (107) are equivalent. Note that the last term can be
rewritten as a dyadic function,

$\bm{\nabla}_{x^{\prime}}g_{\phi}(\bm{x}^{\prime},\bm{x}^{\prime\prime})\bm{\nabla}_{x^{\prime\prime}}\cdot=\overleftrightarrow{\left(\bm{\nabla}_{x^{\prime}}\bm{\nabla}_{x^{\prime}}\cdot\right)}g_{\phi}(\bm{x}^{\prime},\bm{x}^{\prime\prime}).$
(113)

As $g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)$ depends only on
$|\bm{x}-\bm{x}^{\prime}|$, we find

$\displaystyle\bm{\nabla}_{x}\cdot\int
d\bm{x}^{\prime}g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)\bm{j}^{\rm(T)}(\bm{x}^{\prime};\omega),$
$\displaystyle=$ $\displaystyle-\int
d\bm{x}^{\prime}\left(\bm{\nabla}_{x^{\prime}}\cdot
g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)\right)\bm{j}^{\rm(T)}(\bm{x}^{\prime};\omega),$
$\displaystyle=$ $\displaystyle\int
d\bm{x}^{\prime}g_{\bm{A}}(\bm{x},\bm{x}^{\prime};\omega)\left(\bm{\nabla}_{x^{\prime}}\cdot\bm{j}^{\rm(T)}(\bm{x}^{\prime};\omega)\right)$
$\displaystyle=$ $\displaystyle 0.$

This follows the transverse vector potential, $\bm{\nabla}_{x}\cdot\bm{A}=0$.
Therefore, the solutions (106) and (107) obtained from the separated Maxwell’s
equations have no inquiries.

These treatments of the source term $\bm{j}$ and the L and T field components
complicate the theoretical framework. In our four-vector representation, such
complexity becomes clear as the two components $(\bm{A},\phi/c)$ of the fields
and the two source components $(\bm{j},c\rho)$ are related to each other in
matrix form.

## References

  * (1) M. S. Tame, K. R. McEnery, Ş. K. Özdemir, J. Lee, S. A. Maier, and M. S. Kim, Nat. Phys. 9, 329 (2013).
  * (2) C. J. Powell and J. B. Swan, Phys. Rev. 116, 81 (1959).
  * (3) L. Novotny and B. Hecht, Principles of Nano-Optics Second Edition (Cambridge University Press, United Kingdom) 2012.
  * (4) A. Otto, Z. Phys. 216, 398 (1968).
  * (5) E. Kretschmann and H. Raether, Z. Naturforsch. 23A, 2135 (1968).
  * (6) P. J. Feibelman, Prog. Surf. Sci. 12, 287 (1982).
  * (7) R. R. Gerhardts and K. Kempa, Phys. Rev. B 30, 5704 (1984).
  * (8) M. A. Garcia, J. Phys. D: Appl. Phys. 45, 389501 (2011).
  * (9) S. Zeng, K.-T. Yong, I. Roy, X.-Q. Dinh, X. Yu, and F. Luan, Plasmonics 6, 491 (2011).
  * (10) M. S. A. Gandhi, S. Chu, K. Senthilnathan, P. R. Babu, K. Nakkeeran, and Q. Li, Appl. Sci. 9, 949 (2019).
  * (11) A. O. Govorov, H. Zhang, H. V. Demir, and Y. K. Gun’ko, Nano Today 9, 85 (2014).
  * (12) M. L. Brongersma, N. J. Halas, and P. Nordlander, Nat. Nanotech. 10, 25 (2015).
  * (13) T. P. White and K. R. Catchpole, Appl. Phys. Lett. 101, 073905 (2012).
  * (14) A. O. Govorov, H. Zhang, and Y. K. Gun’ko, J. Phys. Chem. C 117, 16616 (2013).
  * (15) P. V. Kumar, T. P. Rossi, D. Marti-Dafcik, D. Reichmuth, M. Kuisma, P. Erhart, M. J. Puska, and D. J. Norris, ACS Nano 13, 3188 (2019).
  * (16) K. Ueno, T. Oshikiri, and H. Misawa, Chem. Phys. Chem. 17, 199 (2016).
  * (17) W. Li and J. G. Valentine, Nanophotonics 6, 177 (2017).
  * (18) T. Tatsuma, H. Nishi, and T. Ishida, Chem. Sci. 8, 3325 (2017).
  * (19) C. Clavero, Nat. Photo. 8, 95 (2014).
  * (20) R. Sundararaman, P. Narang, A. S. Jermyn, W. A. Goddard III, and H. A. Atwater, Nat. Commun. 5, 5788 (2014).
  * (21) M. Bernardi, J. Mustafa, J. B. Neaton, and S. G. Louie, Nat. Commun. 6, 7044 (2015).
  * (22) H. Zhang and A. O. Govorov, J. Phys. Chem. C 118, 7606 (2014).
  * (23) M. B. Ross and G. C. Schatz, J. Phys. D: Appl. Phys. 48, 184004 (2015).
  * (24) L. V. Besteiro, X.-T. Kong, Z. Wang, G. Hartland, and A. O. Govorov, ACS Photonics 4, 2759 (2017).
  * (25) A. Manjavacas, J. G. Liu, V. Kulkarni, and P. Nordlander, ACS Nano 8, 7630 (2014).
  * (26) K. D. Chapkin, L. Bursi, G. J. Stec, A. Lauchner, N. J. Hogan, Y. Cui, P. Nordlander, and N. J. Halas, Proc. Natl. Acad. Sci. USA 115, 9134 (2018).
  * (27) J. Ma, Z. Wang, and L.-W. Wang, Nat. Commun. 6, 10107 (2015).
  * (28) X. You, S. Ramakrishana, and T. Seideman, J. Phys. Chem. Lett. 9, 141 (2018).
  * (29) R. Zhang, L. Bursi, J. D. Cox, Y. Cui, C. M. Krauter, A. Alabastri, A. Manjavacas, A. Calzolari, S. Corni, E. Molinari, E. A. Carter, F. J. G. de Abajo, H. Zhang, and P. Nordlander, ACS Nano 11, 7321 (2017).
  * (30) J. M. Fitzgerald, S. Azadi, and V. Giannini, Phys. Rev. B 95, 235414 (2017).
  * (31) R. D. Senanayake, D. B. Lingerfelt, G. U. Kuda-Singappulige, X. Li, and C. M. Aikens, J. Phys. Chem. C 123, 14734 (2019).
  * (32) R. Schira and F. Rabilloud, J. Phys. Chem. C 123, 6205 (2019).
  * (33) J. M. Pitarke, V. M. Silkin, E. V. Chulkov, and P. M. Echenique, Rep. Prog. Phys. 70, 1 (2007).
  * (34) J. M. McMahan, S. K. Gray, andG. C. Schatz, Phys. Rev. Lett. 103, 097403 (2009).
  * (35) S. Raza, G. Toscano, A.-P. Jauho, M. Wubs, and N. A. Mortensen, Phys. Rev. B 84, 121412(R) (2011).
  * (36) G. Toscano, S. Raza, A.-P. Jauho, N. A. Mortensen, and M. Wubs, Opt. Exp. 20, 4176 (2012).
  * (37) A. Moreau, C. Ciracì, and D. R. Smith, Phys. Rev. B 87, 045401 (2013).
  * (38) W. Wang and J. M. Kinaret, Phys. Rev. B 87, 195424 (2013).
  * (39) N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, and S. I. Bozhevolnyi, Nat. Commun. 5, 3809 (2014).
  * (40) T. Christensen, W. Yan, Søren Raza. A.-P. Jauho, N. A. Mortensen, and M. Wubs, ACS Nano 8, 1745 (2014).
  * (41) M. K. Svendsen, C. Wolff, A.-P. Jauho, N. A. Mortensen, and C. Tserkezis, J. Phys.: Condens. Matter 32, 395702 (2020), and related reference therein.
  * (42) T. Inaoka, J. Phys.: Condens. Matter 8, 10241 (1996).
  * (43) K. Cho, Optical Response of Nanostructures: Microscopic Nonlocal Theory Springer Series in Solid-State Sciences (Springer-Verlag, Tokyo) 2003.
  * (44) E. Prodan and P. Nordlander, Chem. Phys. Lett. 349, 153 (2001).
  * (45) F. H. L. Koppens, D. E. Chang, and F. J. G. de Abajo, Nano Lett. 11, 3370 (2011).
  * (46) D. Lu, Z. Hou, X. Liu, B. Da, and Z. J. Ding, J. Phys. Chem. C 123, 25341 (2019).
  * (47) Y. Aharonov and D. Bohm, Phys. Rev. 115, 485 (1959).
  * (48) A. Tonomura, N. Osakabe, T. Matsuda, T. Kawasaki, J. Endo, S. Yano, and H. Yamada, Phys. Rev. Lett. 56, 792(1986).
  * (49) R. Kubo, J. Phys. Soc. Jpn., 12, 570 (1957).
  * (50) Complete inclusion of the correlation and quantum fluctuation for the excited states is a task beyond the simple preparation of the exact eigenstates $|\mu\rangle$ including the Coulomb interaction at equilibrium. This point will be a subject of future research based on the present study.
  * (51) T. Yokoyama, M. Iio, T. Inaoka, H. Ishihara, to be submitted.
  * (52) K. Cho, Reconstruction of Macroscopic Maxwell Equations: A Single Susceptibility Theory Springer Tracts in Modern Physics Volume 237 (Springer-Verlag, Berlin Heidelberg) 2010.